Transition radiation in an infinite one-dimensional structure interacting with a moving oscillator—the Green's function method

Transition radiation in an infinite one-dimensional structure interacting with a moving

oscillator—the Green's function method

Faragau, Andrei B.; Mazilu, Traian; Metrikine, Andrei V.; Lu, Tao; van Dalen, Karel N.

DOI

10.1016/j.jsv.2020.115804

Publication date

2020

Document Version

Final published version

Published in

Journal of Sound and Vibration

Citation (APA)

Faragau, A. B., Mazilu, T., Metrikine, A. V., Lu, T., & van Dalen, K. N. (2020). Transition radiation in an

infinite one-dimensional structure interacting with a moving oscillator—the Green's function method. Journal

of Sound and Vibration, 492, 1-22. [115804]. https://doi.org/10.1016/j.jsv.2020.115804

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

ContentslistsavailableatScienceDirect

Journal

of

Sound

and

Vibration

journalhomepage:www.elsevier.com/locate/jsv

Transition

radiation

in

an

infinite

one-dimensional

structure

interacting

with

a

moving

oscillator—the

Green’s

function

method

Andrei

B.

F

˘ar

˘ag

˘au

a ,∗

_,

_Traian

_Mazilu

b

_,

_Andrei

_V.

_Metrikine

a

_,

_Tao

_Lu

a

_,

Karel

N.

van

Dalen

a

a Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, CN Delft 2628, Netherlands b University Politehnica of Bucharest, Splaiul Independentei nr. 313, sector 6, Bucharest, CP 060042, Romania

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 30 April 2020 Revised 2 October 2020 Accepted 20 October 2020 Available online 20 October 2020 Keywords:

Transition radiation Infinite system Inhomogeneous system Green’s function method Moving oscillator Non-reflective boundaries

a

b

s

t

r

a

c

t

Transition zones in railwaytracks areareas with considerable variation of track prop-erties (e.g.,foundation stiffness) whichmaycausestrong amplification oftheresponse, leading torapiddegradationofthe trackgeometry.Two possible indicatorsof degrada-tioninthesupportingstructureareidentified,namelythewheel-railcontactforceandthe power/energyinputbythevehicle.Thispaperanalysestheinfluenceofaccountingforthe interactionbetweenthevehicleandthesupportingstructureonthecontactforceandon the power/energyinput.To thisend,aone-dimensionalmodel isformulated,consisting ofaninfiniteEuler-BernoullibeamrestingonalocallyinhomogeneousKelvinfoundation, interactingwithamovingloadedoscillatorthathasanonlinearHertzianspring.The solu-tionisobtainedusingtheGreen’s-functionmethod.ToobtaintheGreen’sfunctionofthe inhomogeneousand infinitebeam-foundationsub-system,thefinitedifferencemethodis usedforthespatialdiscretizationandnon-reflectiveboundaryconditionsareapplied. Ac-countingforthe interactionbetweenthemovingoscillatorandthe supportingstructure generallyleadstostrongerwaveradiation,causedbythevariationofthevertical momen-tum ofthemovingmass.Resultsshowthatforrelativelyhighvelocityandsmall transi-tionlength,themaximumcontactforceaswellastheenergyinputexhibitasignificant increasecomparedtothemovingconstantloadcase.Furthermore,forrelativelyhigh ve-locities, themaximumcontact forcealsoincreasessignificantlywithincreasingstiffness dissimilarity,findingswhichsupplementtheexistingliterature.Finally,thetwo degrada-tionindicatorscanbeusedinthepreliminarystagesofdesigntoassesstheperformance ofrailwaytracktransitionzones.

© 2020 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Transition zones in railway tracks are areas with substantial variation of track properties (e.g., foundation stiffness) encountered near rigid structures such as bridges, tunnels or culverts. These zones require 3–8 times more frequent

∗_{Corresponding author.}

E-mail address: A.B.Faragau@tudelft.nl (A.B. F ˘ar ˘ag ˘au).

https://doi.org/10.1016/j.jsv.2020.115804

0022-460X/© 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )

maintenance than the regular partsof therailwaytrack[1,2] ,leadingto highcosts andreducedavailability ofthe track. Asubstantialpartofthemaintenanceperformedintransitionzonesisconcernedwithrestoringtheverticalpositionofthe track, whichchanges overtime dueto soilandballastsettlement. Inordertoimprovethedesignoftransitionzonesand propose mitigation measures, theunderlying mechanismsleadingto soilandballast settlement needto be better under-stood.Tothisend,some studiesfocusontheballast modelling[e.g.,3–8 ],andothersonthevehicle-structure interaction. Thispaperaddressesthelattercase.

Transition radiationis the radiation emitted when a source (without an inherent frequency) moves rectilinearly with constantvelocityandactsonornearaninhomogeneousmedium[9,10] .Itoccurs,forexample,whenatrainapproachesand crossesaninhomogeneityintherailwaytrack(i.e.,atransitionzone).Whenthishappens,thetrainmustspendadditional energytomaintainitsconstantvelocity[11] .Transitionradiationisthustheradiationcausedsolelybytheinhomogeneity, which is emitted into the railway track in the form of waves.Transition radiation [9,10] has been addressed as one of the primary causes of track and foundation settlement in transition zones due to the often strong amplification ofthe stress andstrain fields[12–15] .Initspurest form,transitionradiationis addressedinstructuresactedupon by amoving constantload.Thefirststudyontransitionradiationofelasticwavesinanelasticallysupportedinfinitestringsubjectedtoa constantmovingloadwascarriedoutbyVesnitskiiandMetrikin[16] .Toalsoaccountfortheflexuralrigidity,Vesnitskiiand Metrikin[11] analysedtransitionradiationinasemi-infinitebeam.Lateron,theproblemofabeamrestingona piecewise-homogeneous Winklerfoundation subjectedto a moving constant load wasaddressed usingdifferent solution methods: modalexpansiontechniques[17,18] andthemoving-elementmethod[19] .Tostudywavepropagationinthegrounddueto transitionradiation,2-Dmodels ofapiecewise-homogeneouscontinuum acteduponby amoving load[20–23] aswell as 3-Dmodelsofthecompletetrack/soilsystem[24,25] wereanalysed.

Modellingthevehicleasamovingconstantloadhasobviouslimitations,suchastheneglectionofthevehicleinertia,the absenceofthecorrespondinginteractionbetweenthevehicleandthesupportingstructure,andtheimpossibilityofunstable motions ofthevehicle.Forthat reason,researchers includedtheinteraction betweenthe vehiclewithits owndegreesof freedom andthesupporting structure instudieson transitionradiation.LeiandMao[26] found that thecontactforce is not greatlyinfluenced bythemagnitudeofthestiffnessjump(i.e.,stiffnessdissimilarity),butthatitis mostlyinfluenced by theangleoftherailelevationafterdifferentialsettlementshaveoccurred.AlsoBanimahdetal.[27] concludedthat the so-calledfaultsinthetrack(i.e.,differentialsettlements)influencethecontactforcemuchmorethanthestiffnessvariation does.Moreover,thefindingsofAngandDai[19] areinagreementwiththeconclusionsofLeiandMao[26] andofBanimahd etal.[27] .However,Germonpré etal.[28] foundthatforlargevelocities,thecontactforceinthecaseofanabrupttransition significantly increasescompared to the smooth transition. Transition radiationwas alsostudied using complex multiple-degrees-of-freedomvehicles interactingwithinhomogeneous1-D supportingstructures [29–31] ,andwithinhomogeneous 2-Dand3-Dsupportingstructures[e.g.,32–34 ],confirmingthattheinertiaofthevehiclecanplayanimportantroleinthe transitionradiationphenomenonandinthetrack-geometrydegradation.

Apartfromthetransitionradiationstudies,theproblemofanoscillatorinteractingwithasupportingstructurehasbeen extensively studied fromthe stabilitypoint of view. The wave-induced instability of an object moving on a distributed elastic systemwasforthefirst time describedby Denisov etal.[35] andBogacz etal.[36] .The parametricinstability of a mass moving ona periodicallyinhomogeneous elastic systemwasthen addressedby VesnitskiiandMetrikin [37] .The physical interpretationofthe causeofinstability,namelytheexcitation ofanomalous Doppler wavesinthe super-critical regime,wasfirstgivenbyMetrikin[38] .Theinstabilityofvibrationswasalsostudiedforarandomly-inhomogeneouselastic supportstructure[39] ,foranaxiallycompressedbeamrestingonviscoelasticfoundation[40] ,andstructuresincontactwith an elastichalf-space[41,42] .TheinstabilityofvibrationsconsideringtheHertziancontactrelationbetweenthevehicleand the supporting structure wasinvestigated by Maziluetal. [43–45] .Furthermore, Abe etal.[46] describedthe instability of avehicle movingon adiscretely supported infinitebeam, whilean instability studyincorporatingtheeffect ofa two-parameterfoundationispresentedbyDimitrovová[47] .

Tostudythetransitionradiationphenomenon,detailedmodelsoffiniteandinhomogeneoussystemsthatinteractwith complex vehiclemodelsare availableinthe literature,asexplainedabove.However, fundamentalstudiesusingsimplified modelsoftransitionzonesininfinitesystemsinteractingwithamovingoscillator,similartothestabilitystudies,arescarce. This motivatesthe aimofthe currentwork, whichis toformulate a 1-D modelof an infiniteEuler-Bernoulli beamona smoothly inhomogeneousKelvinfoundation,interactingwithamovingoscillator,inordertostudytheeffectofaccounting forthevehicle’sinternaldegreeoffreedomontransitionradiation.Thesolutioncanbeobtainedbymeansofthe Green’s-functionmethod[48] becausethesupportingstructureisassumedtobehavelinearly.Tothisend,theresponseofthebeam to the moving vehicle, which ismodelled asa loaded oscillator, isexpressed through a convolution integral in terms of the unknown contact force andthe known time-domain Green’s function ofthe beam-foundation sub-system.Then, the contact relation (i.e.,integral equation) is solved iteratively forthe contactforce. Moreover, the finitedifference method is usedforthe spatialdiscretizationto accommodatethe smoothly inhomogeneousfoundation. The infiniteextentofthe systemisaccountedforthroughasetofnon-reflectiveboundaryconditions,andthroughtheinitialconditionsbasedonthe steady-state responseof abeamwithhomogeneous foundationsubjecttoa moving load(the steadystate isassumedto havedevelopedpriortoreachingthetransition).

The novelty ofthe current work lies in the studyinto the influence of the degree of freedom of the vehicle on the transition radiation in an infinite andsmoothly inhomogeneous structure by comparing it to the situation in whichthe vehicleactionisapproximatedbyitsdeadweight.Furthermore,twopossibleindicatorsofthetransition-zoneperformance

are identified. The current paperhas some similaritiesto the worksof Lei and Mao [26] and of Ang andDai [19] , but is distinct inthe following ways.Firstly, the vehiclemodels inboth studies are more complexthan the oneused in the currentpaper,making theobservations difficulttorelate/compare tothe caseofa singlemovingconstant load. Secondly, both worksarelimitedto piecewise-homogeneousfoundations,while thecurrentpaperaddressesasmooth stiffnessand damping variation ofthe supporting structure.Thirdly,the infiniteextent ofthe systemis properlyaccounted forinthe presentwork.Finally,theextentoftheparametricstudypresentedinbothprevious worksisconsiderablyenlargedinthe currentstudy,andthus,moregeneralobservationsofthesystem’sbehaviourcanbemade,relevantforengineeringpractice. Modellingtherailwaytrackasa1-Dsystem(withfrequency-independentparameters)hasitsshortcomings.Forexample, thewave propagationintheverticaldirectionaswellassurfacewavesarenotcapturedbythismodel,aspects whichcan influence thedynamicresponse.However, thequalitative observationsandthe phenomenaidentifiedwiththissimplified modelwillbepresentinamorecomprehensivemodel.Therefore,forthepurposeofstudyingtheinfluenceofthevehicle’s own degree-of-freedomon thetransitionradiationphenomenon, the1-D modelis sufficient.Itmust beemphasized that this paper focuses on transition radiation asthe process of energy emission, which is characterized by the energy flux though aboundarythatencompassestheemitter,withnoemphasisonthetypeofwavesthataregenerated(propagating and/orevanescent).

2. Modelandsolution

2.1. Problemstatement

Inthissection,a1-Dmodelisformulated,consistingofaninfiniteEuler-Bernoullibeamrestingonasmoothly inhomo-geneousKelvinfoundation,interactingwithamovingone-massloadedoscillator(Fig. 1 ).Theinfinitesystemisdividedinto threedomains:thecomputationaldomain(i.e.,x∈[0,L]ascanbeseeninFig. 1 )thatcontainstheentireinhomogeneityof the supportingstructureanditsvicinity, andtwohomogeneous semi-infinitedomains,oneateach boundaryofthe com-putationaldomain.Theinitialposition(i.e.,att₌0)ofthemovingoscillatorisatx₌0anditsoscillationsarerestrictedto thecomputationaldomain.Theequationofmotionforthethreedomainsreads

˜

∂

4_w ˜

∂

x4 +

ρ

˜

∂

2_w ˜

∂

t2 +cd

(

x

)

˜

∂

w ˜

∂

t +kd

(

x

)

w=− 1 EIQ

(

t

)

δ

(

x−

v

t

)

,

∀

x,

∀

t, (1) w

(

x,t

)

=



_w l

(

x,t

)

, x≤ 0, wc

(

x,t

)

, 0≤ x≤ L, wr

(

x,t

)

, x≥ L, Q

(

t

)

=



_Q 0, t ≤ 0, Qc

(

t

)

, 0≤ t ≤ L_v, Q0, t ≥ L_v, kd

(

x

)

=



_k d,l, x≤ 0, kd,c

(

x

)

, 0≤ x≤ L, kd,r, x≥ L, cd

(

x

)

=



_c d,l, x≤ 0, cd,c

(

x

)

, 0≤ x≤ L, cd,r, x≥ L, where ∂_˜˜ ∂x and ˜ ∂ ˜

∂t represent generalizedderivatives withrespect tospace andtime, respectively. Comparedto the

classi-cal derivatives,thegeneralizedonescanhavediscontinuities,andare necessaryinEq. (1) duetothe Diracdeltafunction on the right-hand side, whichis a generalizedfunction [49] .Furthermore,all parameters in Eq. (1) have beenscaled by the beam’s bending stiffness EI;

ρ

representsthe scaled mass per unit length of the beam;kd,l and cd,l are the scaled (homogeneous) foundationstiffnessanddampingoftheleftsemi-infinite domain,respectively;k_d,r and c_d,r representthe same quantities of the rightsemi-infinite domain; kd,c

(

x

)

andcd,c

(

x

)

are the scaled foundationstiffness andfoundation

Fig. 1. Schematization of the model: infinite Euler-Bernoulli beam continuously supported by a smoothly inhomogeneous Kelvin foundation, interacting with a moving loaded oscillator.

dampingofthecomputationaldomain,respectively;Qc

(

t

)

,Q0and

v

representthetime-dependentcontactforceinthe com-putational domain,the steady-statecontactforce (due to thedead loadofthe carbody) andthevelocity of themoving oscillator,respectively;

δ

(

...

)

denotestheDiracdeltafunction.Finally,wl,wrandwcaretheunknowndisplacementsofthe leftandrightsemi-infinitedomains,andofthecomputationaldomain,respectively.Thespaceandtimedependencyofthe unknown displacementsisomittedfrommostexpressionsforbrevity.Itmustbeemphasizedthatthismodeldoesnot in-corporatetheinfluenceofthediscretenatureoftherailsupport.Moreover,byadopting≤ and≥ signsinthedefinitionof

w

(

x_,t

)

_,Q

(

t

)

_,kd

(

x

)

_,andc_d

(

x

)

foralldomains itisemphasizedthat thereiscontinuityinthesequantitiesattheinterfaces betweenthethreedomains.

Inthisformulation,thevehicleismodelledasaloadedone-massoscillator(i.e.,loadedmass-springsystem),wherethe constantdeadloadQ0 (whichincludesthedeadweightofthewheel)correspondstothevehicleweightonthewheel, the mass M represents the wheelofthe vehicleand thespring representsthe contactelasticitybetweenthe wheelandthe rail. Thecar bodyandthe bogiescould be modelledby additionalsuspended masses. However,the frequenciesat which thesesuspended massesoscillate aremuch lower thantheoscillation frequencyofthewheelin theinteraction withthe rail.Therefore,thedegreesoffreedomofthesuspendedmassesofthecarbodyandbogiesareneglected.Intheremainder ofthepaper,thetermoscillatorisusedtodenotetheloadedone-massoscillator.Theequationofmotionfortheoscillator reads

Mu¨=Q

(

t

)

− Q0, (2)

whereuisthedisplacementofthemassandoverdotsdenoteclassicalderivativeswithrespecttot.Theinteractionbetween theoscillatorandthebeam-foundationsub-systemisdescribedbytheHertziancontactlaw,wherethepossibilityofcontact lossisincorporatedaspresentedinthefollowingequation:



_Q

₍

_t

₎

CH



2 3

=

(

w

(

vt,t

)

− u

(

t

)

)

H

(

w

(

vt,t

)

− u

(

t

)

)

, (3)

whereCH isHertz’sconstantandH

(

...

)

denotestheHeavisidefunction.Itmustbenotedthatforlargetransitionlengths, the inertiaofthesuspendedmassesmayhavea non-negligibleeffectonthebehaviour dueto thelowfrequencies ofthe interaction[32] .

From Eq. (1) it can be deduced that the shear force is continuous everywhere except under the moving load where the force is discontinuous. Thisdiscontinuity isaccounted forinthe equation of motion(Eq. (1) ), for all time moments, dueto thepresenceofthegeneralizedderivatives.Consequently,asinterface conditionsbetweenthedomains, continuity in displacementandslope, aswell asinshearforce andbending momentis imposed.The set ofboundaryconditions is completed byimposingthat,duetothepresenceofdamping,thedisplacementsoftheleftandrightdomainsarezeroat infinitedistancefromthemovingoscillator:

wl

(

0,t

)

=wc

(

0,t

)

, wc

(

L,t

)

=wr

(

L,t

)

, (4) w_l

(

0,t

)

=wc

(

0,t

)

, wc

(

L,t

)

=wr

(

L,t

)

, (5) w_l

(

0,t

)

=wc

(

0,t

)

, wc

(

L,t

)

=wr

(

L,t

)

, (6) w_l

(

0,t

)

=wc

(

0,t

)

, wc

(

L,t

)

=wr

(

L,t

)

, (7)



limx→−∞wl

(

x,t

)

=0, t>−∞, lim₍x−vt)→−∞wl

(

x,t

)

=0, t→−∞.



limx→∞wr

(

x,t

)

=0, t<∞, lim₍x−vt)→∞wr

(

x,t

)

=0, t→∞. (8)

whereprimesdenoteclassicalderivativeswithrespecttox.

The computational domaincould be madeverylarge such thatthe influenceofthe initialconditions ontheresponse dies out before theoscillatorreachesthe transitionzone.However, forcomputational efficiency,x₌0should be asclose aspossibletotheinhomogeneity(i.e.,thecomputationaldomainshouldbeassmallaspossible),meaningthatthechoice of initial conditions affects theresponse in the transitionzone. Prior to reaching a transition zone, the response caused by atrain cangenerallybeconsideredto beinthesteadystate.Consequently, theinitialconditionsmustbe chosensuch that thesystemreachesthesteady-stateregime instantaneouslyatthestart ofthe simulation,andare thusbasedon the steady-state fieldofthe approachingload, alsoreferred toastheeigenfield we_{.Fortheconsidered problem, inthesteady} state, themovingoscillatorisequivalenttoamoving constantloadwherethemagnitudeoftheforce isequaltoQ0,and theeigenfieldforsuchasystemreads[50,51]

we

₍

_x,t

₎

₌



A1e−ik e 2(x−vt)+B₁e−ik e 3(x−vt), x<

v

t, A2e−ik e 1(x−vt)+B₂e−ik4e(x−vt), x≥

v

t, (9)

whereA1,A2,B1 andB2 representcomplex-valuedamplitudesandk1e,ke2,ke3 andke4 representthecomplexwavenumbersof theeigenfield.TheirexpressionsaregiveninAppendix A .Itmustbenotedthatwe_{inEq. (9) isreal-valued.}

The initialconditionsbasedontheeigenfieldare correctonlyiftheinitial fields(i.e.,displacementandvelocityfields) are not disturbedbythe inhomogeneitybecause theeigenfieldrelatesto ahomogeneous Kelvin foundation.Consequently,

x=0 should be positioned such that the initial fields basedon the eigenfield are negligiblysmall atthe inhomogeneity (theycannotbepreciselyzeroduetotheeigenfield’sexponentialdecayinspace).Theinitialconditionsoftheleft domain areimposedastheeigenfieldpartbehindtheoscillator,theonesofthecomputationaldomainareimposedastheeigenfield partinfrontoftheoscillator,andtheinitialconditionsoftherightdomainaretrivial:

wl

(

x,t=0

)

=we

(

x,t=0

)

, w˙l

(

x,t=0

)

=w˙e

(

x,t=0

)

, −∞<x<0, (10) wc

(

x,t=0

)

=we

(

x,t=0

)

, w˙c

(

x,t=0

)

=w˙e

(

x,t=0

)

, 0≤ x≤ L, (11) wr

(

x,t=0

)

≈ 0, w˙r

(

x,t=0

)

≈ 0, L<x<∞. (12) The initial conditionsof the oscillatorshould be chosen such that it alsoreaches the steadystate instantaneously at the start ofthe simulation.Consideringthat the contactforce isconstant andequal toQ0 inthe steadystate, theinitial conditionsfortheoscillatoraregivenasfollows:

u

(

t=0

)

=−



Q0 CH



2 3

+we

₍

_0,0

₎

_{, u˙}

₍

_t=0

₎

_{=0. (13)}

Eqs. (1) to(13) constituteacompletedescriptionoftheproblem.Next,thesolutionmethodispresented.

2.2. Solutionmethod

The system describedin Section 2.1 is nonlinear dueto the nonlinear contactrelation (Eq. (3) ). However, the beam-foundationsub-systembehaveslinearly.Thisallowstoobtainthesolutionusingthetime-domainGreen’s-functionmethod [48] whichisexplainednext.

2.2.1. Generalprocedure

The goalof theGreen’s-function methodis toexpressthe displacement ofthebeamat thecontactpoint interms of theunknown contactforceQ

(

t

)

.The obtaineddisplacementcan besubstitutedinthecontactrelation(Eq. (3) ).Then, the displacementofthemasscanbeexpressedintermsofthecontactforce,andtheresultingexpressioncanbesubstitutedin theequationofmotionofthemass.Theresultingintegro-differentialequationcanthenbesolvedfortheunknowncontact force numerically.Thisapproachissuitedforanytypeofvehiclemodel.However,takingadvantageofthelinearity ofthe vehiclemodelusedinthispaper(exceptforthecontactspring),thedisplacementofthevehicleisexpressedintermsofthe unknowncontactforcethroughtheGreen’s-functionmethodtoo.Then,thecontactrelationreducestoanintegralequation forthecontactforce,andissolvediteratively.

Toexpressthe displacementofthe beamintermsofthe contactforce,the forwardLaplace transformisappliedover timetotheequationofmotionofthecomputationaldomain,resultinginthefollowingexpression:

ˆ

wc +[

ρ

s2+cd,c

(

x

)

s+kd,c

(

x

)

]wˆ=FˆcML+FˆcIC, (14) wherewˆcrepresentstheunknownLaplace-domaindisplacementofthecomputationaldomain,ands=

σ

+i

ω

istheLaplace variable,with

σ

a smallandpositive realnumberand

ω

theangularfrequency.Furthermore,FˆML

c representsthe Laplace-domainmoving-loadforcingofthecomputationaldomain,andFˆIC

c representstheLaplace-domaininitial-conditionsforcing ofthecomputationaldomain,whicharegivenby

ˆ FML c

(

x,s

)

=−1vQˆc



x v



e−sxv, (15) ˆ FIC c

(

x,s

)

=



ρ

s+cd,c

(

x

)



wc

(

x,t=0

)

+

ρ

w˙c

(

x,t=0

)

, (16)

where Qˆc represents the Laplace-domain contact force in the computational domain. It must be noted that the contact force Qˆc isyetunknown,while theinitial-conditionsforcing iscompletelydetermined (Eq. (11) ).It isthereforechosento determine theresponsecausedby themoving loadandtheresponsegeneratedby theinitial conditionsseparately,using different approaches. The displacementwˆIC

c generated by theinitial conditionscan be determined inthe Laplacedomain andtheprocedureiselaboratedinSection 2.2.2 .Thedisplacementcausedbytheunknowncontactforcecanbeobtainedin the timedomainby meansoftheGreen’s-functionmethod.Tothisend,theLaplace-domain displacementwˆML

c causedby theunknowncontactforceisexpressedintermsoftheGreen’sfunctionasfollows:

ˆ wMLc

(

x,s

)

= L 0 ˆ g

(

x,

ξ

,s

)

FˆcML

(

ξ

,s

)

d

ξ

, (17)

where

ξ

representsthespatialvariableofintegrationandgˆ

(

x_,

ξ

_,s

)

representstheLaplace-domainGreen’sfunction ofthe inhomogeneousandinfinitebeam-foundationsystemwhichisdeterminedinSection 2.2.2 .Toobtainthetime-domain dis-placementwML

c ,Eq. (17) isexpressedinthetimedomainusingaconvolutionintegral: wML c

(

x,t

)

=  L 0 t 0 g

(

x,

ξ

,t−

τ

)

Qc

(

τ

)

δ

(

ξ

−

v

τ

)

d

τ

d

ξ

= t 0 g

(

x,

v

τ

,t−

τ

)

Qc

(

τ

)

H

(

L−

v

τ

)

d

τ

, t≥ 0, (18) where g representsthe time-domain Green’sfunction andisobtained byapplying theinverseLaplace transformto gˆ_, as showninSection 2.2.2 .

The totaldisplacement ofthe beaminthetime domainis asuperposition ofthedisplacement causedby themoving loadandthedisplacementgeneratedbytheinitialconditions:

wc

(

x,t

)

=

t

0

g

(

x,

v

τ

,t−

τ

)

Qc

(

τ

)

H

(

L−

v

τ

)

d

τ

+wICc

(

x,t

)

+wLc

(

x,t

)

, (19) where wL

c representsthe displacementofthe computationaldomain causedby theoscillatorcontinuingits movement in the rightdomain,andwL

c∼ H

(

v

t− L

)

meaning thatwcL isnon-zero onlyaftertheoscillatorhasentered therightdomain. Becausetheoscillatorisassumedtobeinthesteadystatewhenitenterstherightdomain,wL

cdependsonthesteady-state contactforceQ0andnotontheunknowncontactforceQc

(

t

)

.Therefore,wLc isobtainedusingthesameprocedureemployed forobtainingwIC

c,whichiselaboratedinSection 2.2.2 .

Thedisplacementofthebeamhasnowbeenexpressedintermsoftheunknowncontactforce,asseeninEq. (19) .Taking advantageofthelinearityofthevehiclemodel(exceptforthecontactspring),thesameprocedureisappliedforobtaining thedisplacementofthemassasafunctionoftheunknowncontactforce.TheLaplace-domaindisplacementofthemassuˆ reads ˆ u

(

s

)

=−Q0 Ms3 + ˆ Qc Ms2+ u

(

t=0

)

s + ˙ u

(

t=0

)

s2 . (20)

Thetime-domaindisplacementuofthemassisobtainedbyevaluatingtheinverseLaplacetransformofEq. (20) analytically, anditsexpressionreads

u

(

t

)

=−Q0 2Mt 2₊ 1 M t 0

(

t−

τ

)

Qc

(

τ

)

d

τ

+u

(

t=0

)

+u˙

(

t=0

)

t, t≥ 0. (21) Boththedisplacementofthemassandthatofthebeamarenowexpressedintermsoftheunknowncontactforce.The contactforce Qc

(

t

)

issolvedforiterativelyateach timemomentsuchthat thecontactrelationissatisfied.Theprocedure forobtainingthecontactforceisthoroughlyexplainedin[43] ,andisbriefly summarizedhere.Outsidethecomputational domain,thecontactforceisassumedtobe constantandequaltoQ0.Therefore,thecontactforce needstobedetermined onlyinthetimeintervalwhentheoscillatorisinsidethecomputationaldomain.Uponevaluatingthedisplacementofthe beam at x=

v

t, both displacements (Eqs. (19) and (21) ) are substituted inthe contact relation (Eq. (3) ), resulting inan integralequation:



_Qc

₍

_t

₎

CH



2 3 =R

(

t

)

H

(

R

(

t

)

)

, 0≤ t≤L

v

, R

(

t

)

= t 0 g

(

vt,

v

τ

,t−

τ

)

Qc

(

τ

)

d

τ

+wICc

(

vt,t

)

−

−Q0 2Mt 2₊ 1 M  t 0

(

t−

τ

)

Qc

(

τ

)

d

τ

+u

(

t=0

)

+u˙

(

t=0

)

t

, (22) where H



R

(

t

)



denotes the Heaviside function with argument R

(

t

)

. Considering the time discretization t0,t1,..., tn, with

t0=0,tn=t (thepresenttimemoment)and

t=ti− ti−1wherei=

{

1,...,n

}

,Eq. (22) becomes



_Q c

(

tn

)

CH



2 3 =RnH

(

Rn

)

, 0≤ tn≤L

v

, Rn= n  i=1  ti ti−1 g

(

v

tn−

v

τ

,tn−

τ

)

Qc

(

τ

)

d

τ

+wICc

(

v

tn,tn

)

−



−Q0 2Mt 2 n+ n  i=1 1 M ti ti−1

(

tn−

τ

)

Qc

(

τ

)

d

τ

+u

(

t=0

)

+u˙

(

t=0

)

tn

. (23) It mustbeemphasizedthat althoughEqs. (22) and(23) arevalidonlyfor0_{≤ t}n≤_vL,thesimulationisnotlimitedtothis

timeintervalprovidedthattheoscillationsofthemasshavestoppedbeforetheoscillatorreachestherightboundaryx=L. Assuming that boththe Green’sfunction ofthe beam-foundationsystemandthecontactforce havea linearvariation duringthe timeinterval

τ

_∈[t_i−1,t_i]_, theintegralsinEq. (23) are evaluatedanalytically, resultingina nonlinearalgebraic equation which can be solved iteratively fortheunknown contactforce at presenttimetn.Oncethe time historyof the

contact force has been obtained, the response of the mass and that of the beam can be computed by substituting the determinedcontactforceinEqs. (19) and(21) .

Atthispoint,theGreen’sfunctionoftheinhomogeneousandinfinitebeam-foundationsub-systemg,thedisplacement

wIC

c causedby theinitialconditions,andthedisplacement wLc generatedbythe oscillatorcontinuingits movementinthe rightdomainhavenotyetbeendetermined.Thisisdoneinthenextsub-section.

2.2.2. Thebeam-foundationGreen’sfunction,theresponsetoinitialconditionsandtheresponsecausedbytheoscillator continuingitsmovementintherightdomain

Duetothespatialvariationofthefoundationstiffnessanddamping(asseeninEq. (14) ),theGreen’sfunctiongˆofthe beam-foundation sub-system, the displacement wˆIC

c caused by theinitial conditions,and the displacement wˆLc generated by theoscillatorcontinuingits movementin therightdomaincannot besolved analyticallyforall stiffnessanddamping profiles.Therefore,thefinitedifferencemethod(centraldifferenceschemeofO

(

x2

₎

_{)isusedtospatiallydiscretizeEq.(14) .} Upon discretizationofthecomputationaldomain,theLaplace-domainequationofmotiondescribingthebeam-foundation sub-systemsubjectedtoinitialconditionsreads

[Kˆi j+

(

ρ

s2+cd,c,is+kd,c,i

)

Ii j]wˆICc, j=FˆcIC,i+FˆcD,i, (24)

whereKˆi j representsthebending-stiffnessmatrixofthebeam,FˆcD,irepresentsa forceexertedbythesemi-infinitedomains

attheboundariesderivedinSection 2.2.3 ,andIi j istheidentitymatrix.Thetime-domaindisplacementisobtainedbyleft

multiplyingEq. (24) withtheinverseofthedynamicstiffnessmatrix(expressioninthesquarebrackets)andbyevaluating theinverseLaplacetransformnumerically.However,theLaplace-domainspectrumofwˆIC

c, jexhibitsapoordecayduetothe

appliednon-zeroinitialconditions.Consequently,thenumericalintegrationmustinprinciplebeperformeduptoveryhigh frequenciesleadingtoasignificantcomputationaleffort.Amethodofincorporatingthehighfrequencieswithoutincreasing thecomputationaleffortispresentedin[51] andisalsoappliedheretoaccuratelyobtainwIC

c, j.

ThesameprocedureusedforobtainingwIC

c isusedtoobtainwLc.TheLaplace-domainequationofmotionreads

[Kˆi j+

(

ρ

s2+cd,c,is+kd,c,i

)

Ii j]wˆcL, j=FˆcL,i, (25)

whereFˆL

c,irepresentsaboundary-forcingvectorexertedattherightboundarybytheoscillatorasitcontinuesitsmovement

intherightdomain,andisderivedinSection 2.2.3 .

ToobtaintheLaplace-domainGreen’sfunctiongˆneededinEq. (17) ,theright-handsideofEq. (14) isreplacedbyaDirac deltafunctioninspaceEI−1

δ

(

x−

ξ

)

.Upon spatialdiscretization,theLaplace-domainequationofmotionthatdescribesthe impulseresponsereads

[Kˆi j+

(

ρ

s2+cd,c,is+kd,c,i

)

Ii j]gˆj=Fˆcδ,i. (26)

ThediscretizedDiracdeltafunctionisavectorwithasinglenon-zeroentryatx₌

ξ

_,anditsexpressionisgivenasfollows:

ˆ

Fδ_c =



0,0,· · · , 1

EI

x,· · · ,0,0



T

. (27)

It mustbe notedthat Eq. (26) hasto besolved individually foreach positionoftheimpulse

ξ

_∈[0_,L].Consequently, the Green’s functions are storedin a matrixgˆi j wherei indicates the discretized excitation variable

ξ

i and j thediscretized

observationvariablexj.Inaddition,thelengthofthediscreteelementattheboundariesis

x/2,whichhastobeaccounted

forinEq. (27) when

ξ

_i=0or

ξ

_i=L.

SpecialattentionshouldbegiventothenumericalevaluationoftheinverseLaplacetransformoftheGreen’sfunctions. Firstly,inordertoevaluatetheGreen’sfunctionsatxj=

v

tjand

ξ

i=

v

τ

i(Eq. (22) )thesamplingrequirements

x=

v

t and

ξ

=

v

τ

needtobesatisfied.Furthermore,dependingontherequiredmaximumfrequencyofintegration

ω

max,thetime stepcan becomeverysmallleadingto anunrealisticallysmallspatialstep.Tofulfilltheabove-mentionedrequirements,a spatialinterpolationschemeisusedforgˆi j.Secondly,forthetotalresponsewc toinstantaneouslyreachthesteadystateat

t₌0_,theGreen’s functionsmustsatisfy trivialinitial conditionsatt₌

τ

_,especiallyclosetot₌0.Tothisend, thecausal character oftheresponseisemployed andtheimaginary-partapproachisselectedsinceitautomaticallysatisfiesthetrivial initial conditions even if the inverse Laplace integral is truncated. Consequently, the time-domain Green’s functions are computedas gi j=− 2eσ (tj−τi)

π

 _ωmax 0 Im

(

gˆi j

)

sin



ω

(

tj−

τ

i

)



d

ω

, i=1,...,j; j=1,...,n. (28) Atthispoint,theonlyunknownsaretheboundaryconditionsofthecomputationaldomain.Uponapplicationofthefinite difference scheme,theboundary conditionsare incorporatedinthestiffnessmatrixofthebeamKˆi j andinthe

boundary-forcingvectorsFˆD

c,iandFˆcL,i.Thenon-reflectiveboundaryconditionsaredeterminedinthenextsub-section.

2.2.3. Non-reflectiveboundaryconditions

The procedureinthissub-sectionis basedontheformulationdescribedthoroughlyin[51] ,andissummarizedin the following. As boundaryconditions forthe computational domain,the reaction forcesof thesemi-infinite domains atthe interfaceswiththecomputationaldomainareimposed.Thegoalistoexpresstheinterfacereactionforces(bendingmoment andshearforce)ofthe leftandrightdomainsasfunctionsoftheunknown displacementandslopeofthecomputational domainatthecorrespondinginterfaces.Tothisend,theforwardLaplacetransformisappliedovertimetotheequationsof motionofthetwosemi-infinitedomains,(Eq. (1) forx∈

(

−∞,0]andx∈[L,∞

)

),leadingto

ˆ

w_l − k4

wherewˆ_landwˆrrepresenttheunknownLaplace-domaindisplacementsoftheleft andrightdomains,respectively;kland

krrepresentthewavenumbersofthetwosemi-infinitedomainsandreadkh= 4



−

ρ

s2_{− c}

d,hs− kd,h withh=

{

l,r

}

,where

the branchesofthecomplex-valuedwavenumbersare chosen such thatIm

(

kh

)

<0 andRe

(

kh

)

>0.Furthermore,FˆlIC and ˆ

FML

r representtheLaplace-domaininitial-conditionsforcingoftheleftdomainandtheLaplace-domainmoving-loadforcing actingon therightdomain,respectively,andtheir expressionsare thesameasEqs. (15) and(16) ,butwithQ0 insteadof

Qc



x

v



andwiththeinitialconditionsoftheleftdomain(Eq. (10) )insteadofthoseofthecomputationaldomain.Asforthe boundary conditions,thezero-displacements condition(the Laplace-domain counterpartofEq. (8) ) isimposed atinfinity, whiletheunknownLaplace-domaindisplacementandslopeofthecomputationaldomain(theLaplace-domaincounterparts ofEqs. (4) and(5) )areprescribedattheinterfaces.

TheLaplace-domaindisplacementswˆlandwˆrcanbeobtainedanalyticallybysolvingEqs. (29) withtheabove-discussed boundaryconditions.Theinterfacereactionforcesofthetwosemi-infinitedomainsarethenexpressedasfunctionsofthe prescribeddisplacementsandslopes,andtheyread

ˆ w_l

(

0,s

)

ˆ w_l

(

0,s

)

=

kl,Vυ kl,Vφ kl,Mυ kl,Mφ

ˆ wc

(

0,s

)

ˆ w_c

(

0,s

)

− Dl, (30)

ˆ w_r

(

L,s

)

ˆ wr

(

L,s

)

=

kr,Vυ kr,Vφ kr,Mυ kr,Mφ

ˆ wc

(

L,s

)

ˆ wc

(

L,s

)

−

ˆ VL

₍

_s

₎

ˆ ML

₍

_s

₎

, (31)

wheretheentriesofthematricesrepresentthedynamicstiffnesscoefficientsassociatedwithboundaryforcesproportional totheunknowndisplacementandslopeattheinterfaces;subscriptVstandsforshearforce,Mforbendingmoment,_{υ for} translation,and

φ

forrotation.ThecoefficientsaregivenbyEqs.(46)and(47)of[51] .Inaddition,vectorDlincorporatesthe influence oftheinitialconditions onthereactionforces,givingrise toboundaryforcesindependent oftheunknown dis-placementandslopeofthecomputationaldomain,whichareaccountedforinEq. (24) throughFˆD

c,i.Duetothefactthatin

Section 2.2.1 theresponsewˆIC

c causedbytheinitialconditions(Eq. (24) )andtheresponsewˆMLc causedbythemoving oscil-latoraretreatedseparately,vectorDlisonlyincorporatedintheboundaryconditionforthedeterminationoftheresponse causedbytheinitialconditions.Thevectorisgivenbythefollowingexpression:

Dl=

ˆ w_l,p

(

0,s

)

ˆ w_l,p

(

0,s

)

+

kl,Vυ kl,Vφ kl,Mυ kl,Mφ

ˆ wl,p

(

0,s

)

ˆ w_l,p

(

0,s

)

, (32)

wherewˆl,p

(

0,s

)

istheparticularsolutionthataccountsforFˆlICinEq. (29) .Itcanbeobtainedasfollows:

ˆ wl,p

(

0,s

)

= 0 −∞gˆl

(

x−

ξ

,s

)





x=0 ˆ FIC l

(

ξ

,s

)

d

ξ

, (33)

wheregˆlistheLaplace-domainGreen’sfunctionofahomogeneousandinfinitebeam-foundationsystemwiththeproperties oftheleftdomain,anditsexpressionisgivenbyEq.(54)of[51] .TheintegralinEq. (33) canbecomputedanalytically,but thesolutionisnotpresentedhereforbrevity.Furthermore,VˆL_andMˆL _{inEq. (31) aretheshearforceandbendingmoment,} respectively,exertedbythemovingoscillatorontherightboundaryafterithasenteredtherightsemi-infinitedomain,and theirexpressionsread

ˆ VL

₍

_s

₎

₌iQ0



s+

(

1+i

)

kr

v



(

kr

v

+s

)(

kr

v

− is

)

e−svL, MˆL

(

s

)

= −iQ0

v

(

kr

v

+s

)(

kr

v

− is

)

e−sLv. (34)

It must benoted that Eq. (34) iscorrect onlyifthe oscillations ofthevehicle havestopped before reaching x=L(as assumedinSection 2.1 andexplainedinSection 2.2.1 ),suchthat thecontactforce isagainconstantandequaltoQ0;this represents the criterion for choosing L. Also, Eq. (34) is used in the boundary conditions only fordetermining the dis-placement wˆL

c generatedbytheoscillatorcontinuingits movementinthe rightdomain,andnot fordeterminingwˆICc, nor fordetermining gˆ.Itmust alsobe mentionedthat forthe resultspresentedinSection 3 ,Eq. (34) isnot usedbecause the simulations areperformedonlyuntilthemovingoscillatorreachestherightboundary.However, accordingtotheproblem statementinSection 2.1 ,ifacertainanalysisisconcernedwiththeresponseofthecomputationaldomainafterthemoving oscillatorhasenteredtherightdomain(i.e.,freevibration),Eq. (34) mustbeincorporated.Moreover,thepresentedsolution methodcanhandleunstablevibrationsoftheoscillatoronsetbytheinhomogeneity(e.g.,breakingthewave-velocitybarrier whengoingfromthestiff domaintothesoftone).However,inthiscasethemaximumtimeofthesimulationislimitedto themomenttheloadreachestherightboundaryofthecomputationaldomainbecauseEq. (34) isonlyvalidforaconstant contactforce.

Eqs. (30) and(31) representthenon-reflectiveboundaryconditionswhichare prescribedtothecomputationaldomain through theLaplace-domain imageofEqs. (6) and(7) .ThetermDl isaccountedforthroughthe boundary-forcingvector

ˆ

FD

c,i while ˆVL and ˆML are accounted for through the boundary-forcing vector ˆFcL,i, because theyare not proportional to the

unknowndisplacementandarethusconsideredasexternalforces.TheremainingpartsofEqs. (30) and(31) areaccounted forinthebeam’sstiffnessmatrixKˆi j.Theapplicationoftheseboundaryconditionsensuresthatthebehaviourofthefinite

computationaldomainisnotinfluencedbyartificialboundariesandthattheinfiniteextentofthemodelisrespectedinan exactmanner.

3. Resultsanddiscussion

In this section, the proposed solutionis firstly validated by considering two limit cases and comparing the obtained resultstoabenchmarksolution.Then,thetime-domainresponseofthewheelandoftherailunderthewheelarepresented fortwospecificcases,andtheinfluenceoftheoscillatorontheresponseatthecontactpointishighlighted.Subsequently, the influenceofthe oscillatoronthe powerinput intothebeam-foundation systemisdiscussed.From theseresults,two indicatorsofpotentialdamagecausedtothesupportingstructureareidentified,namelythemaximumcontactforceandthe energy/powerinput.Finally,theinfluenceofthetransitionlength,oscillatorvelocity,andstiffnessratioonthetwoindicators isassessedthroughaparametricstudy.

Toinvestigatetheinfluenceofthe transitionlength lt on themaximumcontactforceandtheenergy/powerinput,the spatialprofileofthefoundationstiffnessischosenbasedonasinesquared(depictedinFig. 2 )asfollows:

kd,c

(

x

)

=

⎧

⎪

⎨

⎪

⎩

kd,l, 0≤ x<xtc−l2t, kd,l



1+sin





x−xtc lt + 1 2



π 2



2

(

p− 1

)



, xtc−₂lt ≤ x≤ xtc+l₂t, kd,r, xtc+l2t <x≤ L, (35)

where xtc representsthe locationof thecentreofthe transitionandp isthe stiffnessratiobetweenthestiff and soft domains.Thelocationofthetransitionzoneinsidethecomputationaldomainisdictatedbythespatialextentoftheinitial conditions,asdiscussed inSection 2.1 ,while thevalue ofLischosen suchthat theoscillations ofthemoving masshave stoppedbeforeitreachesthepositionx=L,asdiscussedinSection 2.2.3 .Inaddition,thespatialvariationofthefoundation dampingis chosensuch that thedampingratio

ζ

_,definedsimilarly tothat ofasingle-degree-of-freedomsystem, iskept constant.Consequently,thefoundationdampingprofilereads

cd,c

(

x

)

=2

ζ



ρ

kd,c

(

x

)

. (36)

3.1. Parametervalues

Thechoiceoftwoparameters,namelythestiffnessratiop andthevelocity

v

oftheoscillator,requiresspecialattention. Inthiswork,thestiffnessratiop=kd,stiff

k_d,soft isthestiffnessratiobetweenthestiff andsoftdomainsofthesupportingstructure

(i.e.,notincludingthebeam’sstiffness).Thisisdifferentfromthedissimilarityqoftheverticalstiffnessofthetrack,which isthedifferenceinstiffnessexperiencedbythewheel(i.e.,includingthebeam’sstiffness).Thetwoarerelated(instaticor quasi-staticconditions)throughthefollowingexpression[52] :

q≈ p3

4. (37)

Severalresearchershavemeasuredthevertical trackstiffnessinareaswithsoftsoilsandhavereporteddifferencesinthe vertical trackstiffnessattransitionzonesbetweenq≈ 1.5–3[12,14,53,54] ,whichcorresponds to p≈ 1.5–4.5. Formostof theresultsintheremainderofthepaper, p=5ischosen.Althoughitisonthehighend,itisnotunrealisticallyhigh. The value waschosentoaccentuatethetransitionradiationmechanismssuchthat theyare easierto identify.InSection 3.3.1 ,

p=10(i.e.,q=5.6)ischosen toillustrate thewheel-railseparation.Moreover,inSection 3.5.3 , p=1–20(i.e.,q=1–9.5); althoughfromp=5onwardsthevaluesarelessrealistic(orevennotatallrealisticcloseto p=20),theyarepresentedin ordertoobservetheoveralltrend.

As forthevelocity oftheoscillator,its value ispresentedin thispaperasa fractionofthecritical velocity(i.e., mini-mum phasevelocity cph,min forthis1-D system).TheKelvin foundationisknowntosignificantly overestimatethe critical velocity ofthe railwaytrack;forexample,fortheparameters used inthispaper,c_{ph,min≈ 410}m/s(i.e.,1475 km/h). This

Fig. 2. The spatial profile of the foundation stiffness given by Eq. (35) for different transition lengths l t ; the centre of the transition x tc = 5 m and the

Table 1

System parameters (the overbar represents that the coefficient has not been scaled by the beam’s bending stiffness).

Parameter Symbol Value Unit

Bending stiffness EI 6.42 ·10 6 _Nm 2

Mass per unit length ρ 268.33 kg/m

Dead weight Q 0 80 ·10 3 N

Foundation stiffness k d,l 83.33 ·10 6 N/m 2

Foundation damping ratio ζ 0.25

Wheel mass M 750 kg

Hertz’s constant C H 1.1864 ·10 11 N/m 3/2

is muchhigherthanthemeasured criticalvelocityofrailwaytracksinareaswithsoftsoils,wherethe so-calledRayleigh wave velocity (which usually is thecritical velocity)can be aslow as40–50m/s(i.e., 150–180km/h)[55,56] . Therefore, choosingarealisticvaluefortheoscillatorvelocity(e.g.,100–200km/h)togetherwiththeKelvinmodelthatoverestimates thecriticalvelocitywillleadtoaquasi-staticresponseequivalenttothetrainvelocityof20km/honatrackwithrealistic minimumcriticalvelocityof180km/h.Consequently,inthispaper,theoscillatorvelocityischosenrelativetothec_ph,minto representrealisticbehaviour.Aviablealternativeistoincorporatetherail,sleepers,andballastintothebeamelementwith equivalentmass,bendingstiffness,andinternaldamping[57] .Doingso,thecriticalvelocitywillhavemorerealisticvalues; however,thismodelhasitsownshortcomings(e.g.,obtainingtheequivalentbendingstiffnessandinternaldampingisnot straightforward). Nonetheless, thiswillnot change thequalitative observationssince, asexplainedabove, theratioofthe train’svelocitytothecriticaloneisimportant,andnotsomuchtheabsolutevaluesofthesequantities.

The maximum designvelocity of trainsin areaswithsoftsoils isaround 60–70% ofthe minimumcritical velocity in thetrack.However, inextremecases,trainshaveevenexceeded thecriticalvelocity(e.g.,theSwedishX-2000high-speed train whichrunsalong theWestCoastLine betweenGothenburgandMalmö[56] ).Inthispaper,forSections 3.3 and3.4 , the oscillatorvelocity is chosen as90–95% ofcph,min; althoughthese velocitiesare highfor regular trainoperation, they are chosen tooncemore accentuatethetransitionradiationmechanismssuch that theyare easiertoidentify.It mustbe emphasized that theidentified mechanismsare presentalsoforlower velocities,butthey arejustharderto observe.For Section 3.5 ,theoscillatorvelocity isvariedbetween13–110% ofcph,min.Again,although super-criticalvelocitiesare rarely observedinpractice,thewiderangeofvelocitiesisgiventoanalyzetheoveralltrendinbehaviour.

The parameters which are kept constant throughoutthe presented resultsare givenin Table 1 .The bending stiffness corresponds to the UIC 60 rail andthe mass per unit length includes the rail mass per meter and the contribution of the (concrete)sleepers uniformlydistributed along therail [58] ,such that themodel capturesthe trackresponse atlow frequencies(includingitsfirstnaturalfrequency).

3.2. Validationandconvergence

TovalidatethesolutionderivedinSection 2 ,twolimitcasesareconsideredandarecomparedtoabenchmarksolution. Firstly, tovalidate theGreen’sfunctionsofthebeam-foundationstructureandthesuppressionofinitialtransientsthrough the initialconditionsapplied,the limitcaseofa movingconstant loadisconsidered,wherethe magnitudeofthe loadis equaltoQ0.Inthislimitcase,thedegreeoffreedomofthemassandthecontactequationvanish,whilethetime-domain displacementofthebeam(Eq. (19) )simplifiesto

wA lim

(

v

t,t

)

=Q0 t 0 g

(

v

t,

v

τ

,t−

τ

)

d

τ

+wIC c

(

v

t,t

)

. (38)

Secondly,toincorporatetheiterativeschemeofsolvingthecontactequationinthevalidation,thesamelimitcaseis consid-ered(movingconstantload)bytakingthelimitofthewheel’smassgoingtozero.However,thewheel’smasscannotbeset tozerobecauseitleadstoasingularityintheexpressionofthewheel’sdisplacement,ascanbeseeninEq. (21) .Therefore, thewheel’smassisgradually decreasedtostudytheconvergenceofthesecond limit-casesolutionwB

limasthemassgoes tozero.Thefoundationisconsideredaspiecewisehomogeneousforbothcases.

Thelimit-caseresponsesarecomparedtothesemi-analyticalsolutionforthepiecewisehomogeneoussystemsubjected to a moving constant load. Thisbenchmark solution is obtained by solving a system of two semi-infinite domains with differentfoundationstiffnessbyusingtheFourierTransformovertime.IntheFourierdomain,thedisplacementsofthetwo domains canbe obtainedanalytically byimposinginterface conditionsandtheconditionofzerodisplacementsatinfinity [e.g.,22,50 ,51 ].Toobtainthetime-domainsolution,theinverseFourierintegralisevaluatednumerically.

The errorseA _andeB _{presentedinFig. 3 ,corresponding tothe firstandsecondlimit cases,respectively,aredefinedas} follows:

eh

₍

_x=

_v

_t

₎

₌

|

wbench

(

x=

v

t,t

)

− whlim

(

x=

v

t,t

)

|

|

wbench

(

x=

v

t,t

)

|

, h=

{

A,B

}

, (39)

wherewbench isthebenchmarksolution.Tostudytheconvergenceofthepresentedsolution,thetruncationfrequencyfor the inverse Laplacetransform inthe first limit caseisvaried, while a constant truncationfrequency of fmax=16kHz is

Fig. 3. Error under the moving load for different maximum frequencies (panel (a)) and for different values of the wheel’s mass (right panel); the abrupt transition is at x tc = 10 m , v = 0 . 95 c l,ph,min , and p = 5 ; f max = 8 kHz (panel (b)).

usedforthebenchmarksolution.Notethatbychangingthemaximumfrequencyforthelimitcase,accordingtotheNyquist samplingrule,thetimesteppingalsochanges.

Panel(a)ofFig. 3 showsthatthefirstlimit-casesolutionconvergestothesemi-analyticaloneasthemaximumfrequency isincreased.Theerrorinthesoftdomainissmallerthantheoneinthestiff domainbecausethehigherfoundationstiffness leadsto smallerdisplacementsand, thus,to ahigherrelative error.TheerroreA _{ismainlycausedby two factors,namely} thetruncation(infrequency)oftheinverseLaplaceintegral(Eq. (28) )andthespatialdiscretizationintroducedbythefinite differencemethod.Also,tosatisfy

x₌

v

t (discussedinSection 2.2.2 ),thespatialstepsize

xcanbecomeunrealistically small depending on the time step

t. To this end, a spatial step size which leads to accurate results is chosen (

x= 0.05m),andto satisfy

x=

v

t aninterpolationschemeisused.However, theerrorcausedby theinterpolationscheme is negligible.Itcanbe observedinFig. 3 thatfortruncationfrequency fmax=8 kHz,theerrorislessthan 0.5%,whichis more than satisfactory.Ahigher maximumfrequencyleadsto smaller error;however,the computational effortincreases significantly.Forallotherresultspresentedinthissection,themaximumfrequencywasthereforechosenas8kHz.

In panel (b) of Fig. 3 , the second limit-casesolution wB

lim is observed to convergeto the benchmarksolution for the wheel’s masstending to zero.Moreover, theerrors inthe steady-stateregimes arethe sameforall mass valuesandare equaltotheerrorsinthefirstlimitcaseeA_{.Thissuggeststhattheerrorintroducedbyiterativesolverisnegligible.Itmust} be notedthat error eB _{inthe area of the transitionzone is mainly due to theincapability of considering the truelimit} case(M₌0), anddoesnot implyanythingabouttheerrorinthecases studiedfurther,wheretheinertia ofthewheelis significant.

To investigateifthe whole solutionis validand not justatthe contactpoint, thefull displacement field in thelimit caseA(i.e.,Eq. (38) ,butnotlimitedto x₌

v

t)is presentedfordifferenttimemoments inFig. 4 ,andiscompared tothe benchmarksolutionw_bench

(

x,t

)

.ItcanbeseenfromFig. 4 thattheresultsagreeverywell fromthefactthatthetwolines are indistinguishable.Inpanel(a) itcanbe seenthat theinitialtransientshavebeensuccessfullysuppressedthroughthe implementationofthecorrectinitialconditions(i.e.,theresponseisinthesteady-stateregime);hereitcanalsobeclearly observed whytheinitial conditionsbasedontheeigenfieldmustdecaybeforereaching theinhomogeneity (xtc=10 min this case).Moreover, inpanel (e), the fact that thetwo solutions agree verywell demonstratesthat all waves generated insidethecomputationaldomainpropagateawayfromthetransitionwithnoreflectionattheboundaries.Also,inpanel(f) itcan beseenthattheeigenfield passestotherightdomainwithoutanydisturbance.Therefore,itcanbeconcludedthat thesolutionderivedinSection 2 iscorrect.

Finally,avalidationofthevehiclemodelemployed inthispaperispresentedinFig. 5 .Morespecifically,thesuitability oftheloadedone-massoscillatorforstudyingthewheel-railinteractionisinvestigated.Thenormalizedmaximumcontact force obtained withtheloaded one-mass oscillatoriscompared tothe one obtainedwitha morecomprehensivevehicle model,namelythethree-massoscillatorusedbyAngandDai[19] .Forthehighervelocity(panel(a)),theresultsagreevery well,whileforthelowervelocity(panel(b)),theresultspresentsomedifferences,butnotmajorones.Thiscanbeexplained by thefactthat thelowerthe velocity,thelower thefrequencies ofinteraction,which leadstothe sprungmasseshaving a moresignificantinfluence.Overall, theresultsagreewell renderingtheloaded one-massoscillatorsuitable forstudying wheel-railinteraction.Nonetheless,itmustbenotedthatavehiclemodelwithmultiplecontactpointscaninfluencethe re-sultsmoresignificantly.Also,forlargetransitionlengths,theinertiaofthesuspendedmassesmayhaveamorepronounced influenceduetothelowfrequenciesofinteraction.

3.3. Time-domainresponses

Tostudytheinteractionbetweenthemovingoscillatorandthebeam-foundationstructure,thedisplacementofthemass and that of the beam atthe contactpoint, as well asthe contactforce are presented in Fig. 6 . Two specific cases are considered:theoscillatormovingfromthesoftdomaintothestiff oneandthatmoving fromthestiff domain tothesoft one. Forboth cases,soft-to-stiff andstiff-to-soft,the velocityischosen as90% oftheminimumphasevelocity inthesoft domain (cl,ph,min forthe soft-to-stiff case andcr,ph,min forthestiff-to-soft case),and thetransitionlength lt is chosen as

Fig. 4. The displacement field for the moving constant load obtained in the limit case A (black-continuous line) and with the semi-analytical approach (grey-dashed line) for different time moments: (a) t = 0 s, (b) t = 0 . 0187 s, (c) t = 0 . 0246 s, (d) t = 0 . 0299 s, (e) t = 0 . 0368 s, (f) t = 0 . 0618 s; the abrupt transition is at x tc = 10 m , v = 0 . 95 c l,ph,min , and p = 5 ; f max = 8 kHz.

Fig. 5. The normalized maximum contact force obtained with the loaded one-mass oscillator (asterisks) compared to the results obtained by Ang and Dai

[19] with the three-mass oscillator (circles) forv = 90 m/s (panel (a)) and v = 70 m/s (panel(b)); the parameters used are the ones from Ang and Dai [19] .

0.1m,whichisclosetothepiecewise-homogeneouscase.Toensurethattheinitialdisplacementandvelocityfieldsdonot interactwiththeinhomogeneity,thecentreofthetransitionzonextcispositionedatx=19mforthesoft-to-stiff caseand at x₌11m forthe stiff-to-softcase. Itis clearfromFig. 6 that theresponse reachesthesteady stateinstantaneously as the oscillatorenters the computational domainand thus,that the initial transientshavebeen suppressed.Note that this initial plateauisonlydepictedinFig. 6 todemonstratethattheinitialtransientshavebeensuppressed,andisomittedin thefollowingfiguressinceitdoesnotprovideanyinformationabouttheinteraction.

In thesoft-to-stiff case(left panelsinFig. 6 ), the displacementofboth themass andthebeam(at thecontactpoint) exhibit significant oscillations (panel(a) inFig. 6 ). The oscillationsare causedby the interactionofthe moving oscillator withtheinhomogeneouselasticstructure,wherethreephenomenaplayarole:

1. Bending wavesareexcited astheloadapproachesandpassesthe transitiondueto theeigenfield interactingwiththe inhomogeneity.

Fig. 6. Displacement of the mass (dashed line) and of the beam at the contact point (solid line) (panels (a) and (b)), and the normalized contact force (panels (c) and (d)); l t = 0 . 1 m, x tc = 19 m and v = 0 . 90 c l,ph,min for the soft-to-stiff case (left panels); l t = 0 . 1 m, x tc = 11 m and v = 0 . 90 c r,ph,min for the

stiff-to-soft case (right panels); p = 5 .

2. Avibration ofthemass(i.e.,variation oftheverticalmomentum)induced bybending wavesthatkinematicallyexcite themovingoscillatoratthecontactpoint.

3. Avibrationofthemassinducedbyaparametricvariationofthesystempropertiesasexperiencedbythemoving oscil-lator(i.e.,varyingfoundationstiffness).

The threephenomenadescribedherearecoupledasfollows.Thefirst andsecondphenomenaare completely interde-pendent; strongerwave radiationintroduces more oscillationsatthe contactpoint,leadingto astronger vibrationofthe mass, whichinturn influencesthe wave radiation.Inaddition, althoughstronger vibrationofthe massleadsto stronger wave radiation, implying that the third phenomenon influences the first one, stronger wave radiation does not lead to stronger parametricvariation,meaningthat thefirstphenomenon doesnot influencethethirdone.Furthermore,thefirst phenomenon isconsideredtobe theprimaryone,asitalsotakesplaceinthecaseofa movingconstantload[51] ,while thesecondandthirdphenomenaarecharacteristicofthemovingoscillatorcase.

The strongerradiationcausedbytheoscillationofthemass hasa clearinfluenceonthevariationofthecontactforce (panel(c)inFig. 6 ).Intheconsideredcase,themaximumcontactforceisalmostthreetimeslargerthaninthesteadystate. Moreover,themaximumpeakinthecontactforceisfollowedbyadecreaseincontactforcewhichinthiscasereachesthe value ofzeroforashorttimeinterval.Thismeansthat thewheelloses contactwiththerail.Thelossofcontactbetween thewheelandtherailisdiscussedmorethoroughlyinSection 3.3.1 .Thelargevariationinthecontactforce couldbeone ofthefactorsleadingtodegradationofthefoundationinthevicinityoftransitionzones.Therefore,themaximumcontact force isselected inSection 3.5 asan indicatorofthe damageinthe supportingstructure andthus oftheperformance of thetransitionzone.

Inthestiff-to-softcase(rightpanelsinFig. 6 ),thedisplacementofboththemassandthebeam(atthecontactpoint) exhibitsmalleroscillationscomparedtothesoft-to-stiff case,whichisinagreementwiththeliterature[19,51] .Thiscanbe explainedbasedonthephenomenadescribedpreviously:

a)Thesmallerindisplacementandlessbroadeigenfieldturnsouttoyieldweakerbending-waveexcitationatthe inho-mogeneityandtheexcitedlarge-amplitudewavespropagateintothesoftdomain(relatedtothefirstphenomenon). b) Smalleroscillationsatthecontactpointleadtosmallervibrationsofthemass(relatedtothesecond phenomenon),

whichinturnintroduceweakeradditionalwavesinthebeam(relatedtothefirstphenomenon).

c) Astheoscillatorpassesthetransition,thedisplacementatthecontactpointincreasesduetolowerfoundation stiff-ness;however,theinertiaofthemasskeepsit(initially)atthesamelevel.Consequently,inthestiff-to-softcase,the inertia force pointsupwardswhile thedead loadpointsdownwards,partially cancellingeach other (related tothe thirdphenomenon).

Theeffectofitemsb)andc)areclearlyvisibleinthelowervariationofthecontactforceinthestiff-to-softcase(panel (d)inFig. 6 ).

Itisimportanttonotethatthestiff-to-softcaseanalyzedinthispaperassumesthatthesysteminitiallyisinthesteady state. In reality, the systemreaches the steady state in the stiff domain (e.g., a bridge) only if the stiff domain is long enough.Forashortstiff domainembeddedinasoftdomain(e.g.,aculvert[30] ),theresponsemightnotreachthesteady state,leadingtodifferentbehaviourinthestiff-to-softtransition.

Tohighlighttheinfluenceofthemovingoscillatoronthevibrationsatthecontactpoint,thedisplacementofthebeam atthecontactpointiscomparedtothedisplacementofthebeamunderamovingconstantloadinFig. 7 .Thesteady-state