A lattice model for transition zones in ballasted railway tracks

A lattice model for transition zones in ballasted railway tracks

de Oliveira Barbosa, João Manuel; Faragau, Andrei B.; van Dalen, Karel N.

DOI

10.1016/j.jsv.2020.115840

Publication date

2021

Document Version

Final published version

Published in

Journal of Sound and Vibration

Citation (APA)

de Oliveira Barbosa, J. M., Faragau, A. B., & van Dalen, K. N. (2021). A lattice model for transition zones in

ballasted railway tracks. Journal of Sound and Vibration, 494, 1-22. [115840].

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

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

A

lattice

model

for

transition

zones

in

ballasted

railway

tracks

João

Manuel

de

Oliveira

Barbosa

∗

,

Andrei

B.

F

˘ar

˘ag

˘au,

Karel

N.

van

Dalen

a TU Delft, Faculty of Civil Engineering and Geosciences (CiTG), Department of Engineering Structures, Section of Dynamics of Solids and Structures. Stevinweg 1, 2628 CN Delft, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 20 May 2020 Revised 4 November 2020 Accepted 9 November 2020 Available online 16 November 2020

Keywords:

Moving load Infinite lattice

Support stiffness variation

a

b

s

t

r

a

c

t

The procedurepresented inthiswork aimsatprovidingaframework for studying

set-tlementofballastinzoneswithstiffnessvariationoftherailwaytracksupport.The

pro-posed procedureresults fromexpandingan existinginfiniteperiodicmodelofarailway

tracktoaccountforvariationsinthestiffnessofthefoundation.Ballastissimulated via

alinearlattice,whosedynamicresponsediffersfromthatofacontinuum.Theexpanded

model iscomposedofthreeregions:aleftregion,whichissemi-infiniteandperiodic; a

mid-region,offinitelengthandwherethepropertiesofthefoundationcanchange;anda

rightregion,whichisalsosemi-infiniteandperiodicandwhosepropertiescandifferfrom

thoseoftheleftregion.Theequationsofmotionofthemid-regionarewrittendirectlyin

thetimedomain,withtherailbeingdescribedbyfiniteelements.Ontheotherhand,the

leftandrightsemi-infiniteregionsaretreatedsemi-analyticallyinthefrequencydomain,

and afterwards theirresponsesareconverted tothe timedomain,resulting in

convolu-tionintegralsprescribedattheboundariesofthemid-regionthatsimulatenon-reflective

boundaries. Thefinalmodel onlycontainsthedegrees offreedomcorresponding tothe

mid-region(whichcanbeasshortastheregionwherethestiffnessvariationispresent),

and thatleadstofastercalculationsthaniftheboundarieswereplacedfurtherawayto

dissipateundesired reflections.Themethodiscastinthetimedomain,andallelements

areassumedtobehavelinearly.Inthefuture,themodelwillbeexpandedtoincorporate

non-linearbehaviouroftheballast.Thepresentedmethodisvalidatedbymeansofsimple

examples, andafterwardsappliedtoarealscenarioinwhichaculvertcrossesarailway

track. Aspresented,themethodcanbeusedtostudy thelineardynamicsoftransitions

zones, study mitigationmeasures,and infer about indicators likeforce transmitted and

energydissipated,whichmightbeusefultoassessthedevelopmentofsettlementofthe

ballast.

© 2020TheAuthor(s).PublishedbyElsevierLtd.

ThisisanopenaccessarticleundertheCCBYlicense

(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Previousstudiesandmeasurements indicatethat differentialballastsettlementoccursatafasterpaceinregions where railwaytracksexperiencevariationsintheirsupportstiffnessthaninregionswithhomogeneoussupport[1,2].The acceler-ateddifferentialsettlementdemands formorefrequentmaintenance operations,whichtranslatesintoadditionalcosts and

∗_{Corresponding author.}

E-mail address: j.m.deoliveirabarbosa@tudelft.nl (J.M. de Oliveira Barbosa).

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

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

disruptionsofthetrainschedules.Theseregions ofstiffnessvariationarethereforeofinteresttorailwayresearchers,since understandingthephenomenabehindtheaccelerateddegradationhelpsdesigningsolutionstomitigatetheissue.

Studiesfocusingonsettlementofballastcanbecategorizedintofield measurementsandconditionmonitoring[3,4],in which thepositionof railsandthe geometryofballast aremeasured andmonitoredinthe periodbetweenmaintenance operations;labtests[5–7],inwhichballast layersarecyclicallyloaded andtheevolutionofitsgeometryisrecorded;and numericalsimulations,inwhichmoreorlesscomplexmathematicalmodelsareformulated topredicttheresponseofthe track. Whilstthe first category is very useful fordetecting main features of the settlement, it isless viable to optimize mitigation measures,since anyalterationwouldrequirenewconstructions, andthat wouldfirstrequirean approvalfrom therelevantsafetyauthorities,whichareusuallyveryconservativeandstrict.Ontheotherhand,labtestsallowforabetter control of the loading process and in thisway are useful to derive constitutive laws for the ballast (which can include mitigationmeasures suchasgeogrids[8] orbindingpolymers[9]),butforstudyingrailwaytracksupportstiffnesschanges thisapproachbecomesunfeasibleduetothespacerequiredtoaccommodatesuchfeatures.Numericalsimulationshavethe disadvantageofbelongingtoafantasyworld,inwhicheverythingiswelldefinedandcontrolled,butareveryversatileand thus usefulforassessing theefficiencyofmitigation measures.The threecategoriesofstudiescomplementeach other.In thiswork,thefocusliesonnumericalsimulations.

Awidevarietyofmodelsfocussingonlongitudinalvariationsofthetracksupporthavebeenpresentedintheliterature. Theyincludesimple 1Dmodels ofrails restingonviscoelastic Winklerfoundationswithlongitudinallyvarying properties [10–13],simple2Dmodelsofhalf-spaceswithstepvariationsintheirproperties[14,15],andverydetailed2Dor3Dmodels of track, foundationandpotential mitigation measures [16–19]. Forthe latter,the mostcommonsolution strategy isthe use ofthe finite element method(FEM) todiscretize the domainand its differentcomponents, whicheventually can be describedwithnon-linearconstitutivelawsthatdefinethesettlementbehaviourofballast[20].However,becauseFEMcan only handlefinitedomains, inorderto avoidboundaryeffects, largevolumesof domainhaveto bemodelled, leadingto hugecomputationaltimesandrenderingtheapproachnotsuitableforlargeparametricstudies.

Inwhatconcernsballastdescription,onecanassumecontinuummodelsordiscretemodels.Thefirstassumptionis read-ilyavailableinanyFEMsoftware,andcanbecombinedwithnon-linearconstitutivelawsthatdefinetheballastbehaviour [21].However, representingballast insuch awayfailstocapturecertainfeatures ofits dynamicresponse,namelythe ex-istenceofathresholdfrequencyabove whichwavescannot propagate,aswell asthespecificdispersivebehaviourathigh frequencies;thismayverywellaffectthelocaldynamicbehaviouroftheballastandthusthedevelopmentofpermanent deformations.Discretemodels, ontheother hand,intendtodescribe each particleofballastindependently, andthuscan provide moreaccuratedescriptions ofits behaviour.Thediscrete elementmethod(DEM) [22–26] is themethodusedthe most, sinceitcandescribethecontactbetweenparticles quitewellandcapturedifferentfailure modesofballast,butdue toitshighcomputationalcost,itcannotbeusedtosimulatethelengthsofthetrackneededtoconsiderchangesinthe stiff-nessofthesupport.AlternativelytoDEM,latticemodels[27,28]alsointendtodescribeeachballastparticleindividually,but theyassumethatallparticlesareequalandregularlydistributedandthatcontactsbetweenparticlesdonotchange,inthis wayavoidingthecontactcalculationsandreducing thecomputational time.Bymakingthelattice connectionsnon-linear, settlementbehaviourcanbesimulated.

The procedurepresentedinthisworkaims atprovidingaframework forstudyingtrain-passage inducedsettlement of ballast inzoneswithstiffnessvariationsof thesupport.The purposeofthe framework isto allowparametric studieson mitigation solutions,andthereforethecalculationsshouldnot betoocomputationallyintensive.Therefore,atime-domain 2D modelisproposed,inwhichballastandfoundationare describedbya linearlattice(infuturework,non-linearlattice connections willbe addedto lookatsettlementbehaviour),andtheinfinitecharacter ofthetrackinthe longitudinal di-rectionisachievedwithnon-reflectiveboundaries. Theproposedapproachconsistsindividingtheballastedtrackinthree regions:onebeforethestiffnessvariationthatissemi-infiniteinthelongitudinaldirection;oneafterthestiffnessvariation thatisalsosemi-infinite;andamiddleregionoffinitelengthwherethestiffnessvariationisdescribed.Becausethedomains beforeandafterthestiffnessvariationarereplacedbytransmittingboundaries,thefinalmodelonlycontainsthedegreesof freedom correspondingtothemiddleregion(whichcanbeasshortastheregionwherethestiffnessvariationispresent), andthatleadstofastercalculationsthaniftheboundarieswereplacedfurtherawaytodissipateundesiredreflections.The resultingapproachisanexpansionofthemodelpresentedinapreviouswork[29],inwhichasemi-analyticalsolutionfor thetrackwithoutstiffnessvariationispresented.

Thismanuscriptisorganizedasfollows.Section2 introducesthemodelintermsofthethreeregionsmentionedinthe previous paragraphandformulatestheequationsofmotionofeachregioninthetimedomain.InSection3,theequations derivedinapreviouswork[29] areworkedoutinordertoobtainthetransmittingboundariesthatreplacethesemi-infinite regions. InSection 4,themethodisverifiedby comparingeachintermediate stepwithequivalentmodels,andthe three-regionmodelwiththeanalyticalexpressionsfrom[29].InSection5,theusageoftheproposedapproachisillustratedbased ona realstretchoftheDutchrailwayline.Inthelast section,Section6,theworkissummarizedandfinalconsiderations areproffered.

2. Threeregionmodel

The procedurepresented inthisarticlebuildson themodelintroduced ina previous work[29],inwhich therailway trackisconsideredtobeperiodicandissolvedsemi-analyticallyformovingandnon-movingloads.There,railsaremodelled

Fig. 1. Schematization of the periodic model (adapted from [29] ).

Fig. 2. Left, middle and right regions.

Fig. 3. Mid-region and interaction forces.

asan Euler-Bernoullibeam,sleepers asrigidmasses,andballastandfoundationasahorizontallylayeredlattice, whichis a regular networkofequallyspaced massesconnectedviaviscoelastic connections,asdefinedinthe workby Suikerand collaborators[27,28].Theweightofthelattice massesandthepropertiesoftheviscoelasticconnections maychangefrom layer tolayer inorder todistinguish thedifferentmaterial propertiesofballast andfoundation. Fortheballast layer, the propertiescan be inferredfromlabcompression tests(seeSection 4), whileforsoil/foundationlayers, thepropertiescan be inferredfromequivalentcontinuousmodels[27] (beawarethatlatticeslimit theapproximationofcontinuumlayersto Poisson’sratiobelow0.25,seeSection5.1).Fig.1 schematizestheperiodicmodel.

The expansiondescribedinthiswork consistsinaddingamid-regioninwhichthepropertiesofthefoundation(orof anyothercomponent, thoughheretheinterest isinthefoundation)canvary.Inthisway,themodelisdividedintothree parts: aleft region,whichisperiodicandsemi-infiniteandformulated semi-analyticallyinthefrequencydomain;a mid-region,withinwhichmaterialpropertiescanchangeandthatisformulatedinthetimedomain;arightregion,whichlike the left region issemi-infiniteand formulated semi-analyticallyin thefrequencydomain. The left andrightregions may havedistinctballastorfoundationproperties.Fig.2 schematizesthedivisionofthesub-regions.

In thefollowingsubsections eachofthe sub-regionsisdescribed, andthe procedureto couplethemisexplained.The mid-regionisdescribedfirst.

2.1. Mid-region

Fig. 3 shows the mid-region isolated from the rest ofthe domain. In the figure, the mid-region is represented by 3 regularcellsofthesametypeasthoseoftheleftandrightregions,butinprinciplethemid-regioncanconsistofanytype ofstructure,providedthatthegeometriesatitsleftandrightboundariesarecompatiblewiththesideregions.Inaddition, thematerial propertiesofthemid-regioncanchangefromlocationtolocation,anditslength canassumewhatevervalue necessary. Thisregionis acteduponby external forces(fext,not representedinFig.3), whichinmostcaseswill bea set

ofloads moving ontopofthe railfromone edgetothe other,andbytheinteraction forcesattheleft (f_l) andright(fr)

boundaries.Theseinteractionforcesrepresenttheeffectofthesemi-infiniteregionsonthemiddledomain,andwillensure thatthesimulationisnotaffectedbyunwantedreflectionsattheseedges.Theinteractionforcesoneachsidearecomposed byashearforceandamomentattherail,andacoupleofhorizontalandverticalforcesateachmassofthelattice.

Fig. 4. Left and right regions and interaction forces.

Due to its simplicity, the Finite Elementmethod is used to describe the behaviour of the rail,which in this work is dividedintoelementsthatarehalfthelengthofthelatticecharacteristicdistanced (inprinciple,anydiscretizationlength canbeused,butforconvenienceassociatedwithmeshing,d/2isused).Thus,afterdiscretizingtherailandassemblingall railelementsandviscoelasticconnectionsofthelatticeintothestiffnessanddampingmatrices(KandC),andalltheinertial elementsintothemassmatrixM,theequationsofmotionofthemid-regioncanbewrittenas



_M ll Mlm Mml Mmm Mmr Mrm Mrr









M



_u¨ l ¨ um ¨ ur









¨ u +



_C ll Clm Cml Cmm Cmr Crm Crr









C



_u˙ l ˙ um ˙ ur









˙ u +



_K ll Klm Kml Kmm Kmr Krm Krr









K



_u l um ur









u =fext+



_f l 0 fr



(1)

Forconvenience,thedegreesoffreedomhavebeendividedintothoseattheleftinterface(l),thoseattherightinterface (r),andtheremainingones(m). Thesub-indicesinthematricesK,CandM,andintimedependantvectorsuandf rep-resentthesegroupsofdegreesoffreedom.InEq.(1),thetime dependantvectorucontainstheresponseofeachnode(for nodesoftherailsandsleepers,theresponsesaretheverticaltranslationandtherotation;forthelatticenodes,theyarethe horizontalandtheverticaltranslations).ThemassmatricesM_ll andMrr arezeroeverywhereexceptatentries

correspond-ing tothedegreesoffreedomassociated withtherail (notethatatthelateralboundariesofthemid-regionthereare no latticeparticles;theseareconsideredintheleftandrightregions;seeFig.4 forthedefinitionoftheleftandrightregions). Also, themassmatricesM_lm=MT

ml andMrm=M T

mr arenon-zerobecause theapplicationofthefiniteelement methodto

describetherailleadstoinertialtermscouplingtheleftandrightnodesofeachrailelement.

Eq.(1) issolveddirectlyinthe timedomainusingtheNewmarkmethod[30].Applicationofthismethodleadsto the followingtime-steppingscheme

⎧

⎪

⎪

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎪

⎪

⎩

¯ Kuj_=fj ext+

⎡

⎣

f j l 0 frj

⎤

⎦

+M



a0uj−1+a2u˙j−1+a3u¨j−1



+C



a1uj−1+a4u˙j−1+a5u¨j−1



˙ uj_=a 1



uj_{− u}j−1



_{− a} 4u˙j−1− a5u¨j−1 ¨ uj_=a 0



uj_{− u}j−1



_{− a} 2u˙j−1− a3u¨j−1 (2)

wheretheupperindex j referstothetimeinstantt=j



t (



t istheconstanttime step;thepresenceoftheupperindex meansthevariableistimedependant),K¯ istheequivalentdynamicstiffnessmatrixgivenby

¯

K=a0M+a1C+K (3)

andtheconstantsa0 toa5aregivenbythefollowingexpressions[31]:

a0= 1

β



t2, a1=

γ

β



t, a2= 1

β



t a3= 1 2

β

− 1, a4=

γ

β

− 1, a5=



t

 γ

₂

_β

− 1



(4) The constants

β

and

γ

aretheNewmarkparameters,andaresetat

β

₌0_.25and

γ

₌0_.5,whichleadstounconditionally stablesolvers,atleastinthecaseoflinearandfinitesystems.Inthecurrentcase,theexistenceoftheboundaryforcesfland

fr,whichsimulatesemi-infinitedomains, violatesthefinitesystemcriteria, meaningthatthemethoddoesnot necessarily

havetobe unconditionallystable. ThechoiceoftheseNewmarkparametersassumesthatwithin thetimeinterval



t the

accelerationsremainconstantandequaltotheaveragebetweenthepreviousandthenexttime-step. Tosolvethetime-steppingalgorithminEq.(2) foruj_,u˙j_andu¨j_{,theinteractionforcesf}

landfrmustfirstbedetermined.

Iftheseforceswereresultingfromsomeknownexcitation,ofwhichoneknewtheirtimeevolution,thenthesystemcould be solvedasitis.However,theseforcesaredependanton theresponseoftheboundaryinsucha waythat the displace-ments ofthethreeregions mustbecompatibleattheinterfaces,andthus thesetermsmust befurtherexpanded.Thatis addressedinthefollowingsubsections.

2.2. Leftandrightregions

Fig. 4 showsthe left andrightregions isolated fromthe mid-region.Theyare actedupon by the interactionforces fl

andfronthecorrespondingsides(now actingintheoppositedirection:action-reactionpair),andby otherexternalforces

fext,landfext,r(notrepresentedinFig.4),withforcesmovingalongtherailbeingthemostcommonexternalexcitation.The

responseoftheboundariesofthesedomainscanbewrittenas

u_α

(

t

)

=− t

−∞Hα

(

t−

τ

)

fα

(

τ

)

d

τ

+uα,ext

(

t

)

(

α

=l,r

)

(5)

where H_α

(

t

)

is the impulse response matrix ofthe left (

α

=l) and right (

α

=r) boundaries(vertical displacement and rotationoftherailandhorizontalandverticaltranslationsoflatticemassesattheboundariesduetosuddenunit-impulse forcesappliedatt₌0andatthesameboundaries),andu_α,ext

(

t

)

istheresponseoftheseboundariesduetotheexternally

appliedforcesfext,landfext,r.Theminussigninfrontoftheintegralisduetotheaction-reactionpairattheinterface.Also,

u_α

(

t

)

correspondstoulandurinEq.(1),sincetheremustbecompatibilityofdisplacementsattheinterface.

Duetothediscretenatureofthetimesolverusedforthemid-region,theforcesf_α

(

t

)

canonlybeknownattheintervals

t=j



t:f_αj =f_α

(

j



t

)

.Thus,inordertoevaluatetheintegralinEq.(5),atimeevolutionhastobeassumedfortheseforces betweenthediscretetimeinstants.Here,theapproachproposedbyBodeetal.[32] isfollowed.Thatapproachassumesthat withinatimeintervaltheforceremainsconstantandequaltothemeanoftheforceatthebeginningandattheendofthe timeinterval,i.e.,

f_α

(

t

)

≈f

j−1 α +f_αj

2 ,

(

j− 1

)



t<t<j



t (6)

Whencomparedtootherassumptionsfortheforceevolution,suchaslinearvariation,constantandequaltoforceat be-ginning,orconstantandequaltoforceatend,theapproximationproposedin[32] istheonethatshowedmorerobustness intermsofstabilityofthesolver,andthusistheoneadoptedinthiswork.Nonetheless,thisapproximation(oranyother) isnot consistentwiththeassumptions madewiththeNewmarkmethodusedtosolvethemid-region(inthemid-region, assumptionsaremadeabouttheaccelerationevolution,whileatthesideregions,assumptionsaremadeabouttheevolution oftheinteractionforces),andthatcanleadtounwantedreflectionsattheboundaries.Thesereflectionscanbeminimized byreducingthetimestep,asdiscussedlaterinSection4.

BasedontheapproximationinEq.(6) andtheexpressionsinEq.(5),theboundaryresponsesatthetimeinstantt₌j



t

aregivenby u_αj =u_α

(

j



t

)

=− j  i=−∞  it (i−1)tHα

(

j



t−

τ

)

d

τ

fi−1 α +fiα 2 +uα,ext

(

j



t

)

=−1 2 _t 0 H_α

(

τ

)

d

τ

f_αj − j−1  i=−∞ 1 2  _t −tHα

(

(

j− i

)



t+

τ

)

d

τ

f i α+uα,jext (7)

whichcanberewrittenas

u_αj =−F0 αfαj− j−1  i=−∞ F_αj−ifi α+uα,jext (8)

wherematricesF_αj arecalculatedwith

F_αj =1 2

_t −t

H_α

(

j



t+

τ

)

d

τ

(9)

Eq.(8) splits theinterface displacements inthree components:thedisplacements induced by theforcesactingonthe sametimeinterval;theaccumulateddisplacementsduetoforcesactingontheprevioustime intervals(historyterm);and thedisplacementsinducedbytheexternallyactingforces.Thatequationprovidestheframeworktostudytheleftandright regions, butexpressions forthedisplacementsu_α,j_ext andfortheimpulse responsematricesH_α

(

t

)

andF_αj still needtobe defined. Duetothesemi-infinitenatureofthesidedomains,thesequantitiescanbe easiercalculatedsemi-analytically in thefrequencydomainandlateronconvertedtothetimedomain.SuchisexplainedinSection3,butfirstthecouplingof thesideregionstothemid-regionisexplained,sothatthewholeproblemcanbesolved.

2.3. Couplingthesideregionstomid-region

Eq.(8) canberearrangedtoisolatef_αj ontheleft-handside,leadingto

f_αj =−



F0 α



−1



u_αj+ j−1  i=−∞ F_αj−ifi α− uα,jext



(10)

If Eq. (10) is now inserted into the time-stepping algorithm (2), after some rearrangement, namely moving the factor

(

F0

α

)

−1uαj totheleft-handside,oneobtainsanexpandedtime-steppingalgorithmoftheform

⎧

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎩

ˆ¯ Kuj_=fj ext+

⎡

⎣

¯f j l 0 ¯fj r

⎤

⎦

+M



a0uj−1+a2u˙j−1+a3u¨j−1



+C



a1uj−1+a4u˙j−1+a5u¨j−1



˙ uj_=a 1



uj_{− u}j−1



_{− a} 4u˙j−1− a5u¨j−1 ¨ uj_=a 0



uj_{− u}j−1



_{− a} 2u˙j−1− a3u¨j−1 f_αj =−



F0 α



−1uαj +¯fαj,

α

=l,r (11) wherevectors¯fj l and¯f j r aredefinedby ¯fj α=−



F0_α



−1



_j −1  i=−∞ F_αj−ifi_α− uj α,ext



(12)

andwherethenewequivalentstiffnessmatrixKˆ¯ isgivenby

ˆ¯ K=a0M+a1C+

⎡

⎣

Kll+



F0 l



−1 Klm Kml Kmm Kmr Krm Krr+



F0 r



−1

⎤

⎦

(13)

In short,time-steppingalgorithms (2)and(11) arevery similar, withthedifference that some oftheterms haveslightly differentexpressions, namelytheequivalentstiffnessmatrix,whichinEq.(11) containsstiffnessterms

(

F0

l

)

−1 and

(

F0r

)

−1

associated withthe leftandrightregions,andtheinterface force vectors,whichforthealgorithm inEq.(11) contain the history ofwhat happenedat theboundary beforet=j



t, andare givenby Eq.(12). The extraexpression in Eq.(11) is neededtocalculatetheinteractionforcesf_ljandfrj,whichcanonlybeevaluatedaftertheinterfacedisplacementsuljandu

j r

havebeencomputed.

Note that, as presented,the mid-region is acted upon by external forces f_extj and incident displacements u_α,j_extat the left and rightinterfaces; both thesequantities are prescribed by the user. Whereasu_α,j _ext can be calculated taking into consideration the possible interaction withan oscillator moving on the rail of the side regions, the external forces f_extj cannotconsidersuchinteractionwithoutfirstaugmentingthesystemofEqs.(11) withtheequationsgoverningthemotion oftheoscillator.Thisinteractionwithanoscillatorwillbeconsideredinafuturework;inthiswork,f_extj andu_α,j _extrepresent theactionofamovingconstantload.

3. Semi-analyticalevaluationofinterfacematricesandvectors

Section 2 explainsthemainframeworkofthemodelproposed inthiswork. Itreliesontheknowledgeoftheimpulse responsematricesH_α

(

t

)

anddisplacementsu_α,ext

(

t

)

,bothreferringtotheinterfacesbetweensideregionsandmid-region.

Duetothesemi-infinitecharacteristicoftheformer,thecalculationofthementionedquantitiesisnotstraightforward,and itiseasierdonefirstsemi-analyticallyinthefrequencydomain,andthenconvertedtothetimedomain.Thepresentsection explainshowthisisachieved,startingwiththematricesH_α

(

t

)

.

3.1. Calculationofimpulseresponsematrices

The impulseresponsematricesH_α

(

t

)

give theresponseofall degreesoffreedom atthe interfaceoftheleft (

α

=l) or right(

α

=r)region,duetounit-impulseforcesappliedatanydegreeoffreedomofthesameinterface(forthecalculationof thesematrices,theleft/rightregionstandsalone;apartfromtheappliedimpulseforces,thatinterfaceisforce/stress-free). Theirfrequency-domaincounterparts,H˜_α

(ω

)

,aretermeddynamicflexibilitymatricesandcontaintheharmonicresponseof alldegreesoffreedomduetounitharmonicforcesappliedatanydegreeoffreedom.Therowsofthematricesrefertothe degreeoffreedomatwhichtheresponseisobserved,whilethecolumnsrefertothedegreeoffreedomwheretheforceis applied.MatricesH_α

(

t

)

canbeobtainedfromH˜_α

(ω)

throughaninverseFouriertransform

H_α

(

t

)

= 1 2

π

_+∞ −∞ ˜ H_α

(

ω

)

eiωt_d

_ω

₍₁₄₎

andmatricesF_αj,asdefinedinEq.(9),canbecalculatedvia

F_αj = 1 2

π

 _+∞ −∞ ˜ H_α

(

ω

)

sin

(

ω



t

)

ω

eiωjtd

ω

(15)

The dynamicflexibility matrices H˜_α can be calculated semi-analytically based on the expressions derived in [29].In that referencework,thedisplacements,internal shearandinternal momentoftherail,displacementsofthesleepers,and displacementsofthelatticearegivenforarbitraryforcesactingontherailoronthelattice,assuminginfiniteperiodicity.To gofromtheseinfinite-domainsolutionsandobtainthoseofasemi-infinitedomainofthesametype(towhichthedynamic flexibility matricesH˜_α refer), one must impose that there is no transmission offorces betweenthe left and right semi-infinitedomains that composethe infinitedomain.Thisconditionisequivalenttocreatinga semi-infinitedomainwitha freeboundary(asdesired),andcanbeachievedbyaddingfictitiousforcesofappropriateamplitudeatthepartoftheinfinite domainthatistoberemoved(i.e.,forthecalculationofH˜lthefictitiousforcesareappliedtotherightofthefreeboundary;

forthecalculationofH˜r,totheleft).ThisprocedureisexplainednextforH˜l(forH˜r,theprocedureisverysimilaranddoes

notrequirefurtherexplanation).

Consider a periodically infinite domain asrepresented in Fig. 1. At each vertical alignment going through the lattice masses, it ispossible to apply unit forces Frail _{andmoments M}rail _{at the rail, andhorizontal F}latj

x and verticalFzlatj unit forces at the lattice masses, and for each of these loads, the expressions fromwork [29] can be used to calculate the displacements v rail_,rotations

_θ

rail_{,internalshears}rail_andmomentmrail _{oftherailatanyposition,andalsothehorizontal}

xlatj_andverticalzlatj_{displacementsofthelatticemassesatanyverticalalignment(jstandsforthedepthindexofthelattice} mass). Applythementioned setofforcesatalignmenti indicatedinFig.1 (the forcesattherailshould beappliedatan infinitesimaldistancetotheleftofthealignment),andcollectthedisplacementsofthesamealignmentinmatrixUi,iand the displacementsofalignment ii(alsorepresentedinFig. 1)inmatrixUii,i (firstsubscript refers tothealignment where theresponseiscalculated;secondsubscriptreferstowheretheloadisapplied).Forexample,matrixU_i,iisoftheform(the partofthesuperscriptafterthecommaineachofthematrixentriesrepresentstheloadlocationanddirection)

Ui,i=

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

v

rail,Frail i,i

v

rail,M rail i,i

v

rail,Flat1 x i,i

v

rail,Flat1 z i,i · · ·

v

rail,FlatN x i,i

v

rail,FlatN z i,i

θ

rail,Frail i,i

θ

rail,M rail i,i

θ

rail,Flat1 x i,i

θ

rail,Flat1 z i,i · · ·

θ

rail,FlatN x i,i

θ

rail,FlatN z i,i xlat1,Frail i,i xlat1,M rail i,i x lat1,Flat1 x i,i x lat1,Flat1 z i,i · · · x lat1,FlatN x i,i x lat1,FlatN z i,i

z_ilat_,i1,Frail z_ilat_,i1,Mrail zlat1,Fxlat1

i,i z lat1,Flat1 z i,i · · · z lat1,FlatN x i,i z lat1,FlatN z i,i . . . ... ... ... ... ... ... xlatN,Frail i,i xlatN,M rail i,i x latN,Flat1 x i,i x latN,Flat1 z i,i · · · x latN,FlatN x i,i x latN,FlatN z i,i

z_i,ilatN,Frail z_i,ilatN,Mrail zlatN,Fxlat1

i,i z latN,Flat1 z i,i · · · z latN,FlatN x i,i z latN,FlatN z i,i

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

(16)

Matrix Ui,i doesnot correspond to the dynamicflexibility matrix H˜l because the expressions from work [29] do not

satisfy the free boundary conditionneeded to simulate a semi-infinite domain,i.e., atalignment i there is transmission offorcesfromthe leftto therightside.Theseforces,heretermedconnecting forces,correspond totheinternal shearand moment of the rail at alignment i, and to the resultant forces at the lattice masses dueto the viscoelastic connections betweenalignmentsiandii(seeFig.1).CollectalltheseconnectingforcesinamatrixS_i,ioftheform

Si,i=

⎡

⎢

⎣

srail,Frail i,i srail,M rail i,i s rail,Flat1 x i,i s rail,Flat1 y i,i · · · s rail,FlatN x i,i s rail,FlatN y i,i

m_irail_,i ,Frail m_irail_,i ,Mrail mrail,Fxlat1

i,i m rail,Flat1 y i,i · · · m rail,FlatN x i,i m rail,FlatN y i,i Fi,i

⎤

⎥

⎦

(17)

(the firstsubscript referstothealignment wherethe connectingforcesare calculated,whilethesecondsubscript refersto the alignment where the forces are applied) wheresubmatrix Fi,i, which containsthe connecting forces at the lattice, is calculatedwith



F_i,i F_ii,i



=K˜_i,ii



Ulat i,i Ulat ii,i



(18) InEq.(18),matrixK˜i,iiisthedynamicstiffnessmatrixwhereallthehorizontalanddiagonalconnectionsbetweenalignments

iandiiareassembled,andmatricesUlat

i,i andUlatii,iarethesubmatricesofUi,iandUii,i,respectively,thatcollectonlythelattice displacements,i.e.,afterremovingthefirsttworows.

In orderto compensatefor thenon-zero connectingforces Si,i, asecond set offorces(termedfictitious forces) mustbe applied at alignment iisuch that the combinationof the connecting forces causedby the loads atalignment i andthose causedbythefictitiousforcesiszero,i.e.,(Oisamatrixofzeros)

Si,i+Si,iiAfic=O (19)

MatrixAfic_{,whichcontainstheamplitudesofthefictitiousforcesneededtosatisfythefree-boundaryconditionofthe}

semi-infinitedomain,isasquarematrixofwhicheachcolumncorrespondstothefictitiousforceamplitudesforeachloadscenario atalignmenti.ItisobtainedbyrearrangingEq.(19),anditreads:

Afic₌₋

₍

_S

The dynamic flexibility matrix H˜l of the left domain relatingforces anddisplacements atthe free surface ofthe left

semi-infinitedomainisnowobtainedbyaddingtheinfinite-domaindisplacementsinduced bythefictitiousforcestothose inducedbytheloadsappliedatalignmenti,i.e.,

˜

Hl=Ui,i+Ui,iiAfic

=Ui,i− Ui,ii

(

Si,ii

)

−1Si,i

(21) ThecalculationofmatrixH˜rfollowsthesamestepsasformatrixH˜l,withthealignmentsiandiiswapped,anditsfinal

expressionis

˜

Hr=Uii,ii− Uii,i

(

Sii,i

)

−1Sii,ii (22) Notethattocalculatetheleft andrightflexibilitymatrices,itmaybenecessarytoconsiderdifferentmaterialpropertiesif the leftandrightsemi-infinitedomainsaredifferent.Theexpressions fortheircalculation arestillEqs.(21) and(22),but foreach,thequantitiesmustbecalculatedwiththecorrespondingmaterialproperties.

3.2. AlternativeapproachforcalculatingH˜_α

Thoughthemethodexplainedintheprevioussubsectionworkswell andissolelybasedonthesemi-analytical expres-sionsfromwork[29],thefactthattheseexpressionsrequirethenumericalevaluationofintegralsoverwavenumbersleads to timeconsumingcalculations. Alternatively,afastermethod,which takesadvantageoftheperiodiccharacteristicofthe structure,canbeapplied.Itreliesonthediscretizationoftherailandrecursivecloningofcells,leadingtosomenumerical approximation(andpropagationoferrors).Nevertheless,sincethediscretizationappliedisquitefinecomparedtothe char-acteristicwavelengthintherail(thefocusisonlowfrequencies,below5kHz),theintroducederrorsareminimaland,thus, this alternative approach leadsto virtually the same resultsas themethod presented inthe previous subsection (which introducesnodiscretizationerrors).

ConsidertheportionoftheperiodicallyinfinitedomainrepresentedinFig.1 thatislimitedbythealignmentsiiandiii. After discretizingtherailintosmallbeamelements(say,withlengthl=d/2,whered isthecharacteristicdistanceofthe lattice – see Fig.1), thefrequency-domainequations ofmotionofthe resultingsystem(latticemasses, sleeper massand discretizednodesofthebeam)arewrittenas

⎡

⎣

_K˜K˜_mlll _K˜K˜_mmlm _K˜_mr ˜ Krm K˜rr

⎤

⎦



_u˜ l ˜ um ˜ ur



=

⎡

⎣

˜_f˜f_extext_,,_ml ˜ fext,r

⎤

⎦

(23)

wherematricesK˜_αβ(

α

_,

β

₌l_,m_,r)arethedynamicstiffnessmatrices(K˜_αβ=−

ω

2_M

αβ+i

ω

Cαβ+Kαβ)andwhere˜fext,α are theexternal forcesappliedateachnode. Forconvenience, thedegreesoffreedomhavebeendividedintothosebelonging totheleftedge(alignmentii,denotedwithl),thosebelongingtoalignmentrightedge(alignmentiii,denotedwithr),and thosebelongingtothemiddle(denotedwithm).Takingintoconsiderationthatthemiddlenodesarefreeofexternalforces, i.e.,˜fext,m=0,thenthemiddledegreesoffreedomcanberemovedfromthesystemofequations,obtaininganewrelation

ofthetype



˜ K0 ll K˜ 0 lr ˜ K0 rl K˜0rr



˜ u1 l ˜ u1 r



=



_˜ f1 ext,l ˜_f1 ext,r



(24) wherethesubmatricesK˜0

αβ (

α

,

β

=l,r)areobtainedthroughtheexpression

˜

K0_αβ=

δ

αβK˜_αβ− ˜Kαm



K˜mm



−1K˜mβ (25)

InEq.(25),

δ

_αβistheKroneckerdelta(

δ

_αβ=1if

α

=

β

;

δ

_αβ=0if

α

=

β

),thesuperscript0inmatrixK˜ referstostepzero (at stepn, 2n _{cellsare considered;therecursive procedureis explainednext),andthesuperscript 1atvariables u˜ and}˜f referstothecellnumber(20_{=1).MatrixK˜}0_{relatestheforcesanddisplacementsoftheedgesofonecell;twoofsuchcells}

can beplacednextto eachotherandcoupledwiththeportionofthedomainlimitedbyalignments iandii,obtaining in thiswayanewdomaincomprisingtwocells.Themotionofthistwo-cellsdomainisdescribedbythefollowingsystemof equations

⎡

⎢

⎢

⎢

⎣

˜ K0 ll K˜0lr ˜ K0 rl K˜ 0 rr+K˜couplell K˜ couple lr ˜ Kcouple_rl K˜0 ll+K˜ couple rr K˜0lr ˜ K0 rl K˜ 0 rr

⎤

⎥

⎥

⎥

⎦

⎡

⎢

⎢

⎣

˜ u1 l ˜ u1 r ˜ u2 l ˜ u2 r

⎤

⎥

⎥

⎦

=

⎡

⎢

⎢

⎢

⎣

˜_f1 ext,l ˜ f1 ext,r ˜_f2 ext,l ˜ f2 ext,r

⎤

⎥

⎥

⎥

⎦

(26)

whichrelatesforcesanddisplacementsattheextremitiesofthetwocells. MatrixK˜couple_{representsthedynamicstiffness}

(Eq.(18)), augmentedwiththedegreesoffreedom (andassociatedvalues)corresponding tothenodesofthe rail(oneon each side).Assumingagainthattherearenoforcesappliedinthemiddlealignments, i.e.,˜f1

ext,r=0and˜f2ext,l=0,thenthe

variablesu˜1

r andu˜2l canberemovedfromthesystemofequations,obtaininginthiswayanewsystemofthetype



_K˜1 ll K˜1lr ˜ K1 rl K˜ 1 rr



_u˜1 l ˜ u2 r



=



_˜f1 ext,l ˜f2 ext,r



(27) Ifthisprocedureisappliedrecursivelywiththenewlyobtainedmatrixntimes,oneobtainsasystemofequationsofthe type



_K˜n ll K˜nlr ˜ Kn rl K˜ n rr



_u˜1 l ˜ u2n r



=



_˜f1 ext,l ˜ f2n ext,r



(28) relatingthedisplacementsoftheleftalignmentofthefirstcellanddisplacementsoftherightalignmentofthecell2n_with the forcesapplied at thesame alignments. Thus, by applyingthisrecursive procedure ntimes, one isable todescribe a domaincontaining2n _{cellswithasystemofequationscontainingonlythedegreesoffreedomattheedges ofthedomain.} TheexpressionstorecursivelycalculateK˜n

αβ (

α

,

β

=l,r)are ˜ Kn ll=K˜ n−1 ll −



_˜ Kn−1 lr O



⎡

_⎣

_K˜n−1 rr +K˜couplell K˜ couple lr ˜ Kcouple_rl K˜n−1 ll +K˜ couple rr

⎤

⎦

−1

_

˜ Kn−1 rl O



˜ Kn lr=−



_˜ Kn−1 lr O



K˜n−1 rr +K˜couplell K˜ couple lr ˜ Kcouple_rl K˜n−1 ll +K˜ couple rr



−1



O ˜ Kn−1 lr



˜ Kn rl=−



O K˜n−1 rl



K˜n−1 rr +K˜ couple ll K˜ couple lr ˜ Kcouple_rl K˜n−1 ll +K˜ couple rr



−1



˜ Kn−1 rl O



˜ Kn rr=K˜nrr−1−



O K˜n−1 rl



K˜n−1 rr +K˜ couple ll K˜ couple lr ˜ Kcouple_rl K˜n−1 ll +K˜ couple rr



−1



O ˜ Kn−1 lr



(29)

Since each alignment contains few degreesof freedom (for mostcases not reaching the hundred figure), each evalu-ation of (29)isvery fast,and within afraction ofsecondsthousands ofcellscan be simulated.Also, aslong asthere is some damping(inpractice,thatisalways thecase),matricesK˜n

lr andK˜ n

rl willgoto zero,andmatricesK˜ n ll andK˜

n rr will no

longerchange,meaningthattheleftandrightedgesofthemodelleddomainbecomeuncoupled.Inotherwords,assensed fromeachedge,oneiseffectivelymodellingasemi-infinitedomain.Inthisway,afterconvergenceoftherecursivescheme describedinEq.(29),thedynamicflexibilitymatricesH˜_land/orH˜rarecalculatedwith

˜ Hl=



_˜ Kn rr



−1 ˜ Hr=



_˜ Kn ll



−1 (30)

Thisrecursiveprocedurerevealedtobequitefast,andfortheexampleshowninSection5,ittakeslessthanonesecond to calculate the dynamicflexibility matricesforeach frequency, withn ranging between11(2048 cells) and15 (32,768 cells)forconvergencetobereached;theotherapproachtakesmorethanoneminute.Therelativeerrorbetweenthetwo approacheswasintheorderof10−4,meaningthatthetwoprocedureswillgivethesamevaluesuptothethirdsignificant digit.

Now that the dynamicflexibility matriceshave been derived, the only unknown components left in Eq.(12) are the displacementsu_α,j_extattheedgesofthesemi-infiniteregionscausedbyexternalloads.Thesearereferredtoastheincident displacementsandarederivedinthefollowingsubsection.

3.3. Calculationofincidentdisplacements(movingload)

The incident displacements at the edges of the semi-infinite domains are needed to account for the external forces applied at thesedomains. Forthe problemtackled inthis article,theseforces correspondto loads travelling atconstant speed,firstalongtheleftregionandafteralongtherightregion.Iftheincidentdisplacementswerenotaccountedfor,then atthemomentthattheloadenteredthemid-region,atransientresponsewouldbegenerated,withwavespropagatingboth towards theleft andrightsidesofthe interface andpollutingtheresults (the samewouldhappen whentheloadexited to the rightdomain). Of course,it can be argued that it ispossible to get ridofthe transientresponses by making the mid-regionlongenough,butthatwoulddefeatthepurposeoftheapproachpresentedhere,whichistoreducethenumber ofdegreesoffreedomandthereforethesimulationtime.

Similarly totheprocedureemployed fortheimpulse responsematrices,the incidentdisplacementsu_α,j_extare obtained fromtheirfrequencycounterpart,u˜_α,ext

(ω)

,throughtheexpression

u_α,j_ext= 1 2

π

 _+∞

−∞ u˜α,ext

(

ω

)

e

Inturn,thefrequency-domainincidentdisplacementsu˜_α,ext

(ω)

canbecalculatedbycombiningtheresponsesinducedbya

loadmovingonaninfinitelyperiodicdomainwiththoseofpointloads,insuchawaythatthefreeboundaryisrespected. While theresponseofthe infinitedomainduetoamoving loadcanbe evaluated easily,the responseduetopoint loads still requiresthe evaluationof integralsover wavenumbers, rendering that method lessattractive. However, the stiffness matricesK˜n

ll andK˜ n

rr derivedintheprevioussubsectioncanbeusedtoobtaintheresponseoftheinfinitedomaintopoint

loadswhileavoidingtheintegralevaluation.Thewholeprocessisexplainednextforu˜l,ext

(ω)

.

ConsideraverticalforcemovingwithconstantspeedV alongtherailofthe(infinite)domaindepictedinFig.1.Consider thatthisloadcrossesthealignmenti(fortheleftdomainthatiswherethefreeboundaryhastobeimposed)attimet=t_l. This load induces displacements u˜ML

i at alignment i and displacements u˜MLii at alignment ii (these vectors contain both responsesoftherailandofthelatticemassesatthecorrespondingalignments,asexplainedinSection3.1;expressionsfor theseresponses canbe foundinreference[29]).Alongsidewiththedisplacements, themoving loadalsoinduces internal forcesatalignmenti(connectingforcesbetweentheleftandrightsidesofthealignment),whichcanbecalculatedvia

˜ fML i =



˜ Kcouple_ll K˜couple_lr



u˜ ML i ˜ uML ii



(32) MatricesK˜couple_αβ arethesameasexplainedinSection3.2,andvector˜fML

i containsbothshearforcesandmomentsattherail, andconnectingforcesatthelatticebetweenalignmentsiandii,asexplainedinSection3.1.BecausematricesK˜couple_αβ contain stiffnesscomponentsofadiscretized Eulerbeamrepresentingtherail,some approximationerrorswillbe introduced,but theseareverysmallduetothelongwavelengthoftherailbendingcomparedtothediscretizationlength.

Now,inordertoimposethefreeboundaryconditionatalignmenti,(fictitious)forcesofamplitudesafic_{mustbeapplied}

atalignmentiisuchthatthecombinedeffectresultsintherequiredboundarycondition.Thisstepisthesameasexplained inSection 3.1,withthedifference thatnowafic _{representsavector ofunknown amplitudesandnota matrix,dueto the}

factthatonlyoneloadingscenarioisconsidered(themovingload).Thesefictitiousforcesinducedisplacementsu˜fic i andu˜ficii atalignmentsiandii,andconnectingforces˜ffic

i givenby ˜ ffic i =



˜ Kcouple_ll K˜couple_lr



u˜ fic i ˜ ufic ii



(33) As mentioned already,duetoits computational cost,it isnot desiredto makeuseofexpressions in reference[29] to calculatethedisplacementsinducedbypoint forces.Instead,thematricesK˜n

llandK˜nrr describedintheprevious subsection

are addedtothestiffnessmatricesK˜couple_αβ inordertoobtainthedynamicstiffnessmatrixofthefulldomainassensedby alignmentsiandii,andinthisway,thedisplacementsu˜fic

i andu˜ficii arecalculatedvia(Iistheidentitymatrix,Oisamatrix ofzeros)



˜ ufic i ˜ ufic ii



=



˜ Kcouple_ll +K˜n rr K˜ couple lr ˜ Kcouple_rl K˜couplerr +K˜nll



−1



O I



afic ₍₃₄₎

After substitutingEq.(34) inEq.(33),theconnecting forces˜ffic

i induced bythefictitiousforces(ofyetunknown amplitudes

afic_)aregivenby ˜ ffic i =



_˜ Kcouple_ll K˜couple_lr



K˜ couple ll +K˜ n rr K˜ couple lr ˜ Kcouple_rl K˜couplerr +K˜nll



−1



O I



afic ₍₃₅₎

Thefreeboundaryconditionrequiresthatthereisnotransmissionofforcesbetweenleftandrightsidesofthedomain atalignmenti,whichisachievedbytheidentity˜fML

i +˜f fic

i =0.TakingintoconsiderationEqs.(32) and(35),thatidentityis satisfiedif afic₌₋

⎛

⎝



K˜couple_ll K˜couple_lr





˜ Kcouple_ll +K˜n rr K˜couplelr ˜ Kcouple_rl K˜couplerr +K˜nll



−1



O I

⎞

⎠

−1



˜ Kcouple_ll K˜couple_lr



u˜ ML i ˜ uML ii



(36)

The incident displacementsu˜_l,ext attheleft boundaryare nowcalculated bycombiningthe displacementsinduced by themovingloadandthoseinducedbythefictitiousforces,resultingin

˜ ul,ext=u˜MLi −



I O





˜ Kcouple_ll +_K˜n rr K˜ couple lr ˜ Kcouple_rl K˜couplerr +K˜nll



−1



O I



⎛

⎝



K˜couple_ll K˜couple_lr





˜ Kcouple_ll +K˜n rr K˜couplelr ˜ Kcouple_rl K˜couplerr +K˜nll



−1



O I

⎞

⎠

−1



˜ Kcouple_ll K˜couple_lr



u˜ ML i ˜ uML ii



(37)

Fortheincidentdisplacementsattherightboundary,u˜r,ext,theprocedureisverysimilar,withthedifferencethat

align-menti isswitchedwithalignmentii,andthatthedisplacementsu˜ML

i andu˜MLii arecalculatedforthemovingloadcrossing the alignmentii atthetime instant t=tr=tl+Lm/V,whereLm isthelength ofthe mid-region.The finalexpression for

˜ ur,extis ˜ ur,ext=u˜ML_ii −



O I





˜ Kcouple_ll +_K˜n rr K˜couplelr ˜ Kcouple_rl K˜couplerr +K˜nll



−1



I O



⎛

⎝



K˜couple_rl K˜couplerr





K˜couple_ll +K˜n rr K˜ couple lr ˜ Kcouple_rl K˜couplerr +K˜n_ll



−1



I O

⎞

⎠

−1



_˜ Kcouple_rl K˜couplerr



u˜ML i ˜ uML ii



(38) 4. Example1-Validation

The purposeofthe firstexample showninthiswork isto verifyall stepsof theprocedure.Forthat reason, a simple railwaystructure,consistingofrails,sleepersandballastrestingonarigidfoundationisconsidered.Also,inorderto com-pare theresultsobtainedusingthismethodwiththoseofreference[29],nostiffnessvariationisconsidered.The bending stiffnessandunit mass oftherail areEI_rail=12_.08_{× 10}6_N.m2 _andm

rail=120kg (correspondingtotwo UIC60 rails),the

width,massandrotationalinertiaofthesleepersareWs=0.3m,Ms=315kgandJs=4.73kg.m2,thesleeperspacing (cen-tre tocentre) isL=0.6m,andthethicknessoftheballastlayeris Hballast=0.3m.Railpadsare appliedbetweentherail

andthe sleepers,withverticalandrotationalstiffnesses Kv=108N/mandK_θ=2.14× 105N.m (thesevaluescorresponds

tothepadsunderbothrails,andarebasedonthevaluesreportedin[33];therotationalstiffnesswasobtainedassuming that thepad isequallydistributed alongthecontactwiththesleeper,of lengthWp=0.16m, resultinginK_θ=KvWp2/12),

andcorrespondingviscousdampingCv=105N.s/mandC_θ=4.25× 102N.m.s.Inaddition,under-sleeperpadsareused

be-tweensleepersandballast,totallingaverticalstiffnessofKusp=1.72× 108N/mpersleeper[34] andcorrespondingviscous

dampingCusp=1.72× 105N.s/m(the trackwidthisassumedtobe 2.6m,sameassleeperlength).Thedampingvaluesof

thepadsassume aviscous dampingratioofC_/K=0_.001s−1.Inthefrequencydomain,theviscousdampingisconsidered viacomplexstiffnessesofthetypeK˜

(ω)

=K+iC

ω

.

The propertiesofthe ballastare inferredfromthe boxtestresultsreportedby McDowelletal.[5].There,the authors subjected aballastsamplewhoseparticles sizesranged between25 and50mm (themedianwasapproximately30mm) tosuccessivecompressiveloads,andobserveda ballaststiffnessofabout300kPa/mm.Fromthereportedgranulometry,a particle sizeofd=0.03mwaschosen;thisresultsin11lattice massesper column(the bottomone beingfixed),N=11 lattice massesincontactwitheachsleeper (thestiffnessofthe under-sleeperpadsmustbe distributedby thesemasses) and M=9 free masses in between sleepers. On the other hand, the properties of the particle connections were deter-mined via an inverse problem, in which the springs of a lattice layer with the same thickness as the sample used in the box test (300 mm)were tuned to provide the reported vertical stiffness; together with assuming a Poisson ratio of

ν

=0_.2 [27] and aviscous damping ratioofC_/K=0_.001s−1,the followingvalues were obtained:normallattice stiffness anddampingKn

axi=9.68× 107N/mandCaxin =9.68× 104N.s/m; shearlatticestiffnessanddampingKaxis =1.21× 107N/m

andCs

axi=1.21× 10

4_{N.s/m;diagonallatticestiffnessanddampingK}n

diag=4.23× 10

7_N/mandCn

diag=4.23× 10

4_N.s/m.

Re-gardingthemasspropertiesoftheballast,anhomogenizeddensityof1800kg/m3 _{hasbeenassumed,whichtranslatesinto}

a particle massofm_b=4_.21kg (this mayseema largenumber, butthisvalue accountsforall theparticles atthe same levelalongthewidthofthetrack,whichis2.6m).Forthedefinitionofthelatticeparametersandrailstructureparameters, thereaderisreferredtothepreviouswork[29].

4.1. VerificationofimpulseresponsematricesF_αj

Before looking atthe response induced by a moving load, the intermediate steps are first verified, starting withthe calculation of the impulse response matrices F_αj. With that intention, the procedure explained in Section 3.2 is applied tocalculatethedynamicflexibilitiesH˜_α andthentheir inverseFouriertransformasinEq.(15).Thetimeresponseisthen comparedwiththatsameresponseobtainedwithatimedomainsolver.Forthetimesolver,100cellsoftheperiodicdomain representedinFig.1 aremodelledbasedonafiniteelementdiscretizationoftherail.The100cellsresultin60moftrack modelledoneach sideofthefreeboundaries, whichislongenoughtoabsorb thepossiblereflectionsattheedges ofthe modelleddomain.Atimestepof



t=0.001sisusedforbothapproaches.Thecomparisonsoftheresponses(components ofF_αj)aredisplayedinFig.5a-b(fortheverticaldeflectionoftherailduetoaverticalloadatthesamepoint)andFig.6 a-b (forthehorizontaldeflectionofthemiddlelattice massduetoa horizontalforceatthe sameposition).The subfigures c-d show thecorresponding componentsof thedynamicflexibilities H˜_α (which cannot be compared since they areonly calculatedinthefrequencydomain).

AscanbeseenfromFigs.5–6,thecorrespondencebetweenblueandredlinesinsubfiguresa-bistothepointthat no differences canbe detected between theresponses obtained withEq.(15) andwiththetime domainmethod.There are neverthelesssomesmalldifferencesbetweenthetwo,whichcanbefurtherminimizedbyfurtherdecreasingthetimestep.

Fig. 5. Vertical response of the rail due to vertical load on the rail. a) impulse response of left boundary; b) impulse response of right boundary; c) dynamic flexibility of left boundary; d) dynamic flexibility of right boundary.

Thegoodcorrespondencebetweentheresultsobtainedwiththedifferentmethodsverifiesthevalidityofthisintermediate step.

By comparing Figs. 5a,c with Figs. 5b,d,it is observed that the response is larger for the left boundary than for the right boundary. This is because thespan between the free edge of the rail andthe closest sleeper is larger forthe left boundarythanfortherightboundary(seeFig.4),resultinginamoreflexibleinterface.Therearealsodifferencesregarding thehorizontalresponseofthelatticemassinFig.6 (easiertoseewhencomparingredlinesinsubfigurescandd),which arejustifiedbytheasymmetriccutdecidedfortheboundaries.

4.2. Verificationofincidentdisplacementsu_α,j _ext

Toverifyiftheproceduretocalculatetheincidentdisplacementsattheboundariesiscorrect,theresultsobtainedwith theprocessexplainedinSection3.3 arecomparedwiththoseobtainedwiththesametimedomainmethodasinSection4.1. The externalloadisassumedvertical,ofunitmagnitudeandtravellingalongtherailatconstantspeedV=60m/s.Fig.7 shows the comparisonfor thevertical deflection of therail at the boundaries, andFig. 8 shows the comparisonfor the horizontalresponseofthemiddlelatticemass.

Once again, the match between blue lines and red lines in Figs. 7a-b and 8a-b reveals a very good correspondence betweenthemethoddescribedinSection3.3 andthetimedomainmethod,thusverifyingthecorrectnessoftheformer.In addition,bycomparingtheresponsesoftheleftboundaryandthoseoftherightboundary,verydistinctbehaviourscanbe discernedthatarenotjustifiedsolelybytheasymmetryoftheleftandrightdomains.Thereasonforthesebigdifferences isinthemovingnatureoftheload;whileintheleftdomaintheloadismovingtowardsthefreeedgeandthenleavesit,in therightdomaintheloadsuddenlyentersitandthenmovesawayfromthefreeedge.Thatiswhytheresponseoftheleft boundaryincreasesgradually andthenshowsafreedecaybehaviour(startingthemomenttheloadcrossesthefreeedge), andwhytheresponseoftherightboundaryiszerountil theloadreachestheedge,thenpresentsasuddenincrease,and decaysafterwards.

Fig. 6. Horizontal response of the middle mass of the lattice due to horizontal load at the same mass. a) impulse response of left boundary; b) impulse response of right boundary; c) dynamic flexibility of left boundary; d) dynamic flexibility of right boundary.

4.3. Verificationofthewholeprocedure

The wholeprocess asexplainedinSection2 isverifiednext.Withthat intention,theresponses ofaninfiniteperiodic trackduetoa loadmoving attheconstantspeed ofV=60m/s andunit magnitudeobtainedusingtheproposed method andusingtheequationspresentedinwork[29] arecompared.Fortheproposedmethod,themid-regionconsistsoften iden-ticalcells,andforthetime step,two casesareconsidered:



t=0.001s(asinthepreviousexamples)and



t=0.0001s. Theresultstobecomparedaretheverticaldeflectionsoftherailatthecontactpointswiththesleepers.Becausethe mag-nitudeofthemovingloadisconstantintime,theresponseofanygivenpointattwodistinctcellsmustbethesame(apart fromatimeshift),andthus,theexpressionsin[29] areusedtocalculatetheresponseoftherailonlyatoneposition.

Fig.9 showstheresultsobtainedwiththe twomethods;Figs. 9aand9cshowtheresponses atthe differentpositions over time, asobtained with the proposed method, for the time steps



t=0.001s and



t=0.0001s, respectively, and Figs. 9b and9d show thesame responses forthe sametime steps, butwith time shiftssuch that the periodicityof the response can be assessed.Figs. 9band9d alsodepict indashed blackline thetheoretical responseasobtained withthe expressions from[29].TheverticaldashedlinesinFigs.9aand9crepresenttheinstantsatwhichthemoving loadenters andleavesthemid-region.

FromFig.9aitisobservedthatatthemomentsthattheloadentersandleavesthemid-region,somedisturbancesoccur that affecttheresults.Thisisalsoobservable inFig.9b,whereitcanbe seenthat theresponseofthe railatthesleeper positions is not perfectlyperiodic,i.e., there aresome disturbances around theexpected response(in dashed black line). The reasonforthesedisturbancesistheincompatibilitybetweentheassumptionsmadefortheleftandrightregions and those madeforthemid-region;while fortheformeritis assumedconstantforce duringonetime interval, forthelatter itisassumedconstantacceleration.Byreducingthetimestep,itisexpectedthattheinconsistenciesbecomelessrelevant, and such is proven inFigs. 9c-d, whereit is seen that no disturbances can be seen whenthe load enters orleaves the mid-region,andalsothat thereisavirtuallyperfectoverlapoftheresponses oftherailatdifferentpositions.Despitethe differencesbetweenexpectedandobtainedresponsesfor



t=0.001s,themaximumdeflectionis wellcaptured(only at

Fig. 7. Vertical response of the rail due to vertical load moving along the rail. a) incident displacement at left boundary; b) incident displacement at right boundary; c) frequency response of left boundary; d) frequency response of right boundary.

thefirst sleeperthemaximumdeflectionisslightlysmaller)andtheoverallshapeaswell.Thus, thistime stepisusedin theexamplesolvedinthenextsection,sothatthecalculationsdonotbecometootimeconsuming.

5. Example2– culvertpassage

To demonstrate the capabilitiesof the model,the case ofa railway line passing over a culvert is considered next. It representsatransitionfromasoftfoundationtoaverystiff foundationandagaintoasoftfoundation,andisbasedonan existing stretch ofthe Dutchrailwayline nearthecity of Gouda, whichis thoroughlydescribedin previous publications by otherauthors[35–37].Here,onlytheparametersthatare essentialforthesimulationaredescribed,andthemodelling assumptionsareexplained.

Fig.10 schematizes theconditionsatthesite:a culvert,withlongitudinaldimensioncorrespondingto foursleepers, is placed atthe depth of60cm belowtheballast layer; theculvert rests onpiles, makingit practically immovable(rigid); twoconcreteslabswiththicknessof30cmandlengthcorrespondingto7sleeperslieoneachsideoftheculvertwiththe intentionofmitigatingtheabruptstiffnessvariation;thefoundationconsistsofasandembankmentontopofclayandsand layers.

Forthesuperstructure(rails/sleepers/ballast),theparametersandmodellingstrategies(latticewithcharacteristic diame-ter d₌3cm)arethesameasexplainedinSection4.Regardingtheunderlyingfoundation,itismodelledasalatticewith thesamecharacteristicdiameterd astheballastandwhosepropertiesarediscussedinthenextsubsection.Fortheculvert, since itispractically immovable,itismodelledbyconstrainingthemotionofalllatticemassesthatareatthelocationof theculvert(thussimulatingarigidandimmovablebody).Inturn,theapproachslabsaresimulatedbyassigningthe proper-tiesofconcretetothelatticemassesandlatticeconnectionsatcorrespondinglocations.Thepropertiesassumedforconcrete are Young’smodulus E=20GPa,Poisson’sratio

ν

=0.2andunitdensity

ρ

=2300kg/m3_{,andbasedon thesevaluesand}