Instability of vibrations of an oscillator moving at high speed through a tunnel embedded in soft soil

Delft University of Technology

Instability of vibrations of an oscillator moving at high speed through a tunnel embedded

in soft soil

Zhao, Mingjuan; de Oliveira Barbosa, João Manuel; Yuan, Jun; Metrikine, Andrei V.; van Dalen, Karel N.

DOI

10.1016/j.jsv.2020.115776

Publication date

2021

Document Version

Final published version

Published in

Journal of Sound and Vibration

Citation (APA)

Zhao, M., de Oliveira Barbosa, J. M., Yuan, J., Metrikine, A. V., & van Dalen, K. N. (2021). Instability of

vibrations of an oscillator moving at high speed through a tunnel embedded in soft soil. Journal of Sound

and Vibration, 494, 1-24. [115776]. https://doi.org/10.1016/j.jsv.2020.115776

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

Instability

of

vibrations

of

an

oscillator

moving

at

high

speed

through

a

tunnel

emb

e

dde

d

in

soft

soil

Mingjuan

Zhao

∗

,

João

Manuel

de

Oliveira

Barbosa,

Jun

Yuan,

Andrei

V.

Metrikine,

Karel

N.

van

Dalen

Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, 2628 CN, Delft, The Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 28 March 2020 Revised 3 September 2020 Accepted 7 October 2020 Available online 14 October 2020

Keywords:

Tunnel embedded in soft soil High-speed oscillator Vibration instability Indirect BEM

D-decomposition method Critical velocity for instability

a

b

s

t

r

a

c

t

Thispaperinvestigatestheinstabilityofverticalvibrationsofanobjectmovinguniformly throughatunnelembeddedinsoftsoil.UsingtheindirectBoundaryElementMethodin the frequencydomain, theequivalentdynamicstiffnessofthetunnel-soilsystematthe pointofcontactwiththemovingobject,modelledasamass-springsystemorasthe lim-iting caseofasinglemass,is computednumerically.Using theequivalentstiffness, the original2.5Dmodelisreducedtoanequivalentdiscretemodel,whoseparametersdepend onthe vibrationfrequency and the object’svelocity.The criticalvelocity beyondwhich theinstabilityoftheobjectvibrationmayoccurisfound,anditisthesameforboththe oscillatorandthesinglemass.Thiscriticalvelocityturnsouttobemuchlargerthanthe operationalvelocityofhigh-speedtrainsandultra-high-speedtransportationvehicles.This meansthatthemodeladoptedinthispaperdoesnotpredictthevibrationsofMaglevand Hyperloop vehicles tobecomeunstable. Furthermore, thecritical velocityfor resonance ofthesystemisfoundtobeslightlysmallerthanthevelocityofRayleighwaves,which is verysimilar tothatfor themodelofahalf-spacewitharegular trackplacedontop (withdamping).However,forthatmodel,thecriticalvelocityforinstabilityisonlyslightly largerthanthecriticalvelocityforresonance(oftheundampedsystem),whileforthe cur-rentmodelthecriticalvelocityforinstabilityismuchlargerthanthecriticalvelocityfor resonanceduetothelargestiffnessofthetunnelandtheradiationdampingofthewaves excited inthe tunnel.Aparametricstudy showsthat thethicknessand material damp-ing ratioofthe tunnel,the stiffness ofthe soiland the burial depthhave astabilising effect, whilethedampingofthe soilmay haveaslightlydestabilisingeffect(i.e., lower criticalvelocityforinstability).Inordertoinvestigatetheinstabilityofthemovingobject forvelocitieslargerthantheidentifiedcriticalvelocityforinstability,weemploythe D-decompositionmethodandfindinstabilitydomainsinthespaceofsystemparameters.In addition,the dependencyofthecriticalmassand stiffnessonthevelocityisfound.We concludethatthehigherthevelocity,thesmallerthemassoftheobjectshouldbeto en-surestability(singlemasscase);moreover,thehigherthevelocity,thelargerthestiffness ofthespringshouldbewhenaspringisadded(oscillatorcase).Finally,inviewofthe sta-bilityassessmentofMaglevandHyperloopvehicles,theapproachpresentedinthispaper

∗_{Corresponding author.}

E-mail address: m.zhao@tudelft.nl (M. Zhao).

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

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/ )

canbeappliedtomoreadvancedmodelswithmorepointsofcontactbetweenthemoving objectandthetunnel,whichresemblesrealityevenbetter.

© 2020TheAuthor(s).PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Dynamic effects are ofsignificance formodern high-speedtrains aspropagating wavesmay be generatedin the rail-waytrackandsubsoil.The studyofdynamictrain-track-soilinteractions hasbeen ofinterestforresearchers fordecades. Poppetal.[1]gaveacomprehensivereview oftheexisting modelsthatcanbeusedtostudythedynamictrain-track-soil interaction.

In general,studies onmoving trainsfallinto twocategories.The firstcategoryis environmentalvibrationsinduced by moving trainstoassessvibration hindranceandtoensure thesafetyofthenearbystructures. Thesecond categoryis the instability ofvibrations ofmoving trainsto ensurethesafetyandcomfortofthe passengersinthe trains.Forthe former category, the steady-stateregime is assumedwhen investigating the dynamicamplification dueto resonance [e.g.,2–4]. Other studiesinthefirstcategoryaredevotedtotransitionradiationwhichoccurswhenatrainpassesaninhomogeneity

[5,6].

Forstudiesfallingintothesecondcategory,thetrainisusuallymodelledasasingle-ormulti-degreeoffreedomsystem

[7].Wheninstability occurs,thefreevibration(i.e.,thevibrationintheabsenceofan externalforce)ofthetraingrows ex-ponentially,resultinginaninfinitedisplacementwhentimegoestoinfinity,whichimpliesthatasteady-statesolutiondoes not exist.Thisis verydifferentfromresonance,which happenswhenthe steady-stateresponseasinduced byan external movingloadisextreme(eitherboundedornotdependingonthepresenceofdamping).Bothphenomenacomewitha cer-taincriticalvelocity.Thecriticalvelocityforresonanceisdefinedasthevelocityatwhichthesteady-stateresponseinduced by a moving loadisextreme(i.e., resonancetakesplace atcertain specificvelocities) [8,9],whilethe criticalvelocity for instability isdefinedasthevelocitybeyondwhichinstabilitycanoccur(i.e.,instability occursinarangeofvelocities). An-other crucialdifferencebetweenresonanceandinstabilityisthatresonancecanbetotallyremovedbyincreasingdamping, whiledampingmostlyshiftstheinstabilitydomain,forexample,toaregionoflargervelocities[10,11].

The firststudyoninstability ofvibrationsisthat ofa massthatmovesuniformlyalong anelasticallysupported beam

[10].ThephysicalexplanationofinstabilitywasgivenbyMetrikine[12]whoarguedthattheinstabilityiscausedbythe ra-diatedanomalousDopplerwaves[13]which increase the energy of thevibrating object.In addition,the physical mechanism ofinstabilitywasdiscussedusingthelawsofconservationofenergyandmomentum[12,14].

Afterthepioneeringworksontheinstabilityphenomenon[10,15],severalaspectsthatinfluencethestabilityofanobject havebeendiscussed.Forexample,theeffectofthermalstressesinthestructurewasstudiedconsideringdifferentmodelsof movingoscillators[11,16,17].Otherpapersconsideredtheeffectofmorethanonecontactpointbetweentheobjectandthe structure[17],andofcontactnonlinearities[18].Moreover,amoreaccuratebeammodeltorepresenttherailwasconsidered andacomparisonbetweentheTimoshenkoandEuler-Bernoullibeammodelswasgiven[19].Fourdifferentbeamandplate modelswereconsideredin[14].Furthermore,Verichevetal.introducedothercomplexitiesintheirmodel,i.e.,abogiemodel whichconsistsofarigidbaroffinitelengthontwoidenticalsupports[20].Recently,Mazilustudiedtheinstabilityofatrain ofoscillatorsmoving alongan infiniteEuler-Bernoullibeamonaviscoelasticfoundation[21].Anotherworkfocusedonthe stabilityofa moving massincontactwitha systemoftwo parallelelasticallyconnected beams,withoneofthem being axiallycompressed[22].Later,thestabilityofvibrationsofarailwayvehiclemovingalonganinfinitethree-beam/foundation systemhasbeenconsidered,withan emphasisontheeffectofthedampingandstiffnessofthesecondarysuspension of therailwayvehicle[23].Foramoresimplemodel,Dimitrovová presentedasemi-analyticalsolutionfortheevolutionofthe beamdeflectionshapesandoscillatorvibrations[24].Inthatpaper,notonlytheonsetofinstability,butalsotheseverityis addressed.

All theabove worksuseda one-dimensionalortwo-dimensional modeloftherailwaytrack, whichmaybe less accu-ratethan three-dimensionalmodels [25,26]to predict theinstability ofmoving trains. However,they all conveythevery important messagethat, in thepresence ofdamping, theinstability of moving trainsmayhappen atspeeds that exceed the critical velocity forresonance of theundamped system(which is equal to the minimumphase velocity of waves in thestructure);that is,the criticalvelocity forinstability islargerthanthecriticalvelocityforresonance.Thefewexisting worksrelatedtoinstabilityanalysisemployingthree-dimensionalmodelsoftherailwaytrackconsidertrainsmovingona trackfoundedonthegroundsurface[25,26].Ithasbeenshownthatthecriticalvelocityforinstabilityofthemovingobject is closetotheRayleigh wave speedinthe soil.Instability oftrainsmoving throughan underground tunnelhasnot been analysed yet.Inthispaper,wethereforeaimtoconducttheinstabilityanalyses foranoscillatorandthelimitingcaseofa singlemassmovingthroughatunnelembeddedinsoftsoil.Wewillinvestigatewhetherthecriticalvelocityforinstability of themoving objectis alsocloseto theRayleigh wave speed inthesoil, forboth a shallowandadeep tunnel.The

re-Fig. 1. A 2.5D model of an oscillator (i.e., mass-spring system) moving through a tunnel embedded in an elastic half-space and the associated coordinate systems.

sultsareofpracticalrelevanceespeciallyforcontemporaryhigh-speedrailwaytracksaswellasupcomingultra-high-speed transportationsystemssuchasMaglevandHyperloop,respectively[27–29].

Thepaperisorganisedasfollows.WepresentthemodelandaframeworktoconducttheinstabilityanalysisinSection2.

Section3discussesthe2.5DGreen’sfunctionsofafull-spaceandahalf-space[30,31],presentstheGreen’sfunctionsofthe shell andtheformulationof theindirectBoundary ElementMethod(BEM). Validationsofthe proposedindirect BEMare giveninSection4.InSection5,weconducttheinstabilityanalysisofthesinglemassandthemass-springoscillator.Tothis end, westudytheequivalentdynamicstiffnesstofindthe criticalmassandstiffness,andanalysetheeffectofthetunnel thickness, thematerial dampingratiosinthetunnel-soil system,the Lamé parametersofthesoil andtheburialdepth of thetunnelonthecriticalvelocityforinstability.Moreover,thedependencyofthecriticalmassandstiffnessonthevelocity isinvestigated.ConclusionsaregiveninSection6.

2. Modelandsolutionoftheproblem

2.1. Modeldescription

In this paper,we study thevibrations of an object moving through a tunnel embeddedin softsoil using a so-called 2.5D model.Thesoilismodelledasan elasticcontinuum, whereasthetunnelismodelledbytheFlüggeshelltheory[32]. Both thesoil and tunnelare assumed to be linear,elastic, homogeneous andisotropic. The soil is characterisedby den-sity

ρ

₁ ,Poisson’sratio

ν

₁ ,andcomplexLamé parameters

λ

∗₁ =

λ

1

(

1+2isgn

(

ω

)

ξ

1

)

and

μ

∗1 =

μ

1

(

1+2isgn

(

ω

)

ξ

1

)

,where

i is the imaginary unit,

ω

is the frequency and

ξ

1 the material damping ratio of the soil related to the adopted

hys-teretic damping model. The parameters of the tunnel are density

ρ

₂ , Poisson’s ratio

ν

₂ , andcomplex Lamé parameters

λ

∗

2 =

λ

2

(

1+2isgn

(

ω

)

ξ

2

)

and

μ

∗2 =

μ

2

(

1+2isgn

(

ω

)

ξ

2

)

,with

ξ

2 beingthehystereticmaterialdampingratioofthetunnel.

The burial depthofthetunnelis H,andits innerandouter radiiare R_i andRo .The objectismodelledby amass-spring oscillator (seeFig. 1), whichischaracterisedby its mass MandspringstiffnessK, andmovesthroughthe tunnelwitha constantvelocityV.Notethatthereisnoverticalexternalforceactingonthemassbecausethepresenceofsuchaforceis irrelevantforthedynamic-instabilityanalysis.

Shallowanddeepembeddedtunnelsareconsidered inthispaper.Fortheshallowtunnel,thesoil mediumismodelled asa half-space,whileforthe deeptunnel,thesoilis modelledasafull-space.Fig.1only showsthe configurationofthe shallowtunnel.IfH→∞,itessentiallybecomesadeeptunnel.

Thegoverningequationsoftheshellarepresentedlater,inSection3.Thecurrentsection onlypresentstheframework toconductinstabilityanalysis.

2.2. Methodofsolution

Toanalyseinstability ofvibrations ofthemoving object,theconceptoftheequivalentstiffness(alsoreferredto as dy-namicstiffness)isemployed[33,34].TheprocedureisillustratedinFig.2andgoesasfollows.First,wecomputethe steady-state response of the system shown inFig. 2 (a) which is subjectto a uniformlymoving oscillatory point loadapplied atthe tunnelinvert(r₁ =R_i ,

θ

1 =−π2 ,x=Vt). The oscillatoryloadhas theformofP

(

t

)

=P0 exp

(

i

t

)

, inwhichP0 isthe

amplitude,

₌2

π

f₀ istheangularfrequency;f0istheloadfrequencyinHz.P(t)essentiallyrepresentsaharmonic

interac-tionforcebetweenthemovingobjectandthetunnel-soilsystem.Thesteady-stateradialdisplacementattheloadingpoint can be expressed asUr₁

(

r1 =Ri ,

θ

1 =−π2 ,x=Vt

)

=U0

(

,V

)

exp

(

i

t

)

,where U0 (

, V) is the complexamplitude of this

Fig. 2. (a) A tunnel embedded in an elastic half-space subject to a uniformly moving oscillatory point load P(t) = P 0 exp ( i t) , (b) an equivalent discrete

model consisting of a mass-spring oscillator resting on an equivalent spring K eq ( , V ); is the angular frequency of vibrations and V the velocity of the

oscillator.

presentedindetailinSection3.Fromtheresult,weobtaintheequivalentstiffnessofthetunnel-soilsystemattheloading pointusingthefollowingrelation:

Keq

(

,V

)

=_U P0

0

(

,V

)

. (1)

By doing so, theoriginal 2.5D modelcan be reduced to an equivalentdiscrete model, shownin Fig.2 (b),consisting of a mass-springsystem resting on an equivalentspring with a complex-valuedstiffness K_eq (

, V), which depends on the frequencyandvelocityoftheoscillator.

Tostudytheinstabilityofamovingoscillator,weapply,inaccordancewithpreviousdynamic-instabilitystudies[11,24], theLaplaceintegraltransformwithrespecttotimet(sdenotestheLaplaceparameter)

W

(

s

)

=

∞

0 w

(

t

)

exp

(

−st

)

dt (2)

tothewell-known governingequationoftheverticalmotionoftheoscillator.Assumingzeroinitialconditions(whichcan be done asthey donot influence the stability[11,24]), thefollowing characteristicequation forthefree vibration of the oscillatorisobtained:

Ms2 +_KKKeq

(

s,V

)

+Keq

(

s,V

)

=0. (3)

TherootsofEq.(3)determinethe(complex)eigenfrequencies,

= − is,oftheverticalmotionoftheoscillatorasitinteracts withthe tunnel-soilsystem. Ifone oftherootssofthe characteristicequation hasapositive realpart,theresponse will grow exponentially, which implies that the vertical vibration of the oscillator is unstable [11,19,25,26,35]. Obviously,the equivalentstiffnessmustbesingle-valuedforall

inorderforEq.(3)tobemeaningful;seealsoSection2.3.

It has beenshown in [11]that theinstability of a moving object mayoccur if andonly ifthe imaginary partof the equivalentstiffnessK_eq isnegativeinafrequencyband.In[11],asinglemovingmassisconsidered,themotionofwhichis evennecessarilyunstableassoonasIm(Keq )<0atanyfrequencyband.Theimaginarypartoftheequivalentstiffnesscan beconsideredtobethedampingcoefficientofthedashpotintheequivalentmass-springsystem.Anegativeimaginarypart oftheequivalentstiffnessindicatesanegativedamping,whichmakesthevibrationofthemovingmassunstable.

It would be very laborious to determine all the roots ofthe characteristic equation andcheck whether one of these roots hasapositive realpart.Alternatively, wefollowa convenientmethodofrootanalysis, namelytheD-decomposition method,todeterminethenumberof‘unstableroots’.Thismethodhasbeenusedinseveralpapers[11,19,25,26,35].Theidea ofthismethodistomap theimaginaryaxisofthecomplexsplane(i.e.,theborderbetweenstabilityandinstability)onto theplane ofa systemparameter,MorK,whichisallowed tobecomplex.ThemappedlinedividestheMorKplane into domainswithdifferentnumbersofunstableroots.Itisnotedthattheimaginarypartofthecomplexsystemparameterhas no physicalmeaning.Onlythepositive realpartofthesystemparameter isphysical,andthe crucialquestioniswhether oneoftheso-calledinstabilitydomainsoverlaythepositiverealaxis.

The procedureisasfollows.Considers=i

,where

serves astheparameterofthemapping,isrealvaluedandhas the meaning offrequency(same as introduced above),andhasto be varied fromminus to plus infinity. Wediscuss the followingtwocasesinthispaper.Thefirstoneisthelimitcaseofasinglemassmovingthroughthetunnel,thusassuming K→∞.ThecharacteristicequationforasinglemassisreducedfromEq.(3)to

Ms2 +Keq

(

s,V

)

=0. (4)

Substitutings=i

intoEq.(4)givesthefollowingruleforthemapping:

M= Keq

(

₂ ,V

)

. (5)

Thesecondcaseweconsideristhemoregeneraloneofthemovingoscillator,takingintoaccountboththemassandthe springoftheoscillator.Inthiscase,thestiffnessKwillbeusedastheparameterfortheD-decompositionassumingMtobe

constantbecauseitisofpracticalrelevance.Takingthelimitcaseofthemovingmassasthestartingpoint,itisinteresting toknowwhattheaddedstiffnessofthespringshouldbetorendertheoscillatorvibrationunstable(seealsoSection5.2). Substitutings=i

intoEq.(3),wegetthefollowingmappingruleforthecomplexKplane:

K=M

2 Keq

(

,V

)

Keq

(

,V

)

− M

2 . (6)

Byreplacingsbyi

inKeq (s,V),whichessentiallyentailsconsideringthelimitcaseofs→i

,onecanusethe equiva-lentstiffness(Eq.(1))whichisdeterminedbasedonthesteady-stateresponsetotheharmonicloading.EmployingEq.(5)or

(6)asthemappingrule,onecanplottheD-decompositioncurve,forexample,Im(M)versusRe(M)(asshowninFig.11,for example), where

istherunning parameteralongthiscurve.Oneside oftheD-decompositioncurve isshaded,andthis sideisrelatedtotheright-handsideoftheimaginaryaxisinthesplane.Crossingthecurveinthedirectionoftheshading onceindicatesthat thereisanadditionalunstableroot.Thus,onecanfindinformationontherelativenumberofunstable rootsindomainsofthecomplexMorKplanes.Thenumberofunstablerootsinallthedomainscanbedeterminedifthe absolutenumberofthoseisknownforanyarbitraryvalueoftheconsideredsystemparameter.Bydoingso,theinstability domainscanbefoundintheMorKplane,whichgenerallyallowstoidentifythecriticalvelocityforinstability(definedas thevelocityatwhichaninstabilitydomainfirstoverlaysthepositiverealaxiswhenincreasingthevelocity).Theinstability analysisisconductedinSection5.

2.3. Responsetoamovingoscillatorypointloadandderivationofequivalentstiffness

AsshowninSection2.2,itiscustomarytoemploytheequivalentstiffnessKeq toconducttheinstabilityanalysis.Inthe currentsection, we aimto derive theexpression fortheequivalentstiffness. Tothisend, we firstderive the steady-state response toamoving oscillatorypointloadattheloadingpoint.Additionally,theresponseatafixed observationpoint is derived,whichisneededforvalidatingtheindirectBEMinSection4.Alltheresponsescanbecomputedusingtheindirect BEMpresentedinSection3.Herewesummarisetheimportantstepsandoutcomeinviewofthespecifiedaim.

The shearstresses

σ

r₁θ₁ and

σ

r1x at theinner surfaceof thetunnelwall induced bythe movingoscillatory point load (seeFig.2(a))arezero.Thenon-zeronormalstress

σ

r1r1

(

Ri ,

θ

1 ,x,t

)

canbeexpressedas

σ

r1r1

(

Ri ,

θ

1 ,x,t

)

= P0 Ri

δ



θ

1 +

π

₂



δ

(

x− Vt

)

exp

(

i

t

)

, (7)

where

δ

(.)istheDiracdeltafunction.

Astheconsideredproblemislinear,weapplytheFourierTransformtoderivetheresponseofthesystemsubjecttothe uniformlymovingoscillatory pointloadinthewavenumber-frequency(kx,

ω

)domain.TheFourierTransformappliedwith respecttotimetandspatialcoordinatexisdefinedinthefollowingform(foranarbitraryfunctiong(r₁ ,

θ

₁ ,x,t)):

˜ ˜ g

(

r1 ,

θ

1 ,kx,

ω

)

= ∞ −∞  ∞ −∞ g

(

r1 ,

θ

1 ,x,t

)

exp

(

− i

(

ω

t− kxx

)

)

dxdt (8)

withtheinverseFourierTransformgivenby g

(

r1 ,

θ

1 ,x,t

)

=₄

_π

1₂ _∞ −∞  _∞ −∞ ˜ ˜ g

(

r1 ,

θ

1 ,kx,

ω

)

exp

(

+i

(

ω

t− kxx

)

)

dkxd

ω

. (9)

TheFourierseries,whichisusedtoderivetheresponseinthe(kx,

ω

)domain,ofageneralresponsequantityf(

θ

1 )reads

f

(

θ

₁

)

= n_= ∞ n= −∞ fnexp

(

in

θ

1

)

, fn= 1 2

π

_{2 π} 0 f

(

θ

1

)

exp

(

− in

θ

1

)

d

θ

1 . (10)

Expandingtheterm

δ

(

θ

₁ +π₂

)

inEq.(7)intoaFourierseries,thenormalstresscanberewrittenas

σ

r1r1

(

Ri ,

θ

1 ,x,t

)

= n_= ∞ n= −∞ P0 2

π

Ri exp



in



θ

1 +

π

₂



δ

(

x− Vt

)

exp

(

i

t

)

. (11)

Applying the FourierTransformdefinedbyEq.(8) toEq.(11), thenormalstress inthe wavenumber-frequencydomain is obtainedas: ˜ ˜

σ

r1r1

(

Ri ,

θ

1 ,kx,

ω

)

= n_= ∞ n= −∞ P0 2

π

R_i exp



in



θ

1 +

π

₂



2

πδ

(

ω

−

− kxV

)

=

σ

˜˜aux

(

Ri ,

θ

1 ,kx,

ω

)

2

π

P0

δ

(

ω

−

− kxV

)

. (12) The responseinduced by theauxiliary stress

σ

˜˜aux

(

R_i ,

θ

1 ,kx,

ω

)

,whichrelatesto aradial stress intheformof

δ



θ

1 +π2



·

δ

(

x

)

δ

(

t

)

,canbecomputedusingtheindirectBEM(Section3)andisdenotedasU˜˜_{1 ,aux}

(

r₁ ,

θ

1 ,kx,

ω

)

.Thereafter,wegetthe expressionoftheactualdisplacementvectorexcitedbythestress

σ

˜˜r₁r₁

(

Ri ,

θ

1 ,kx,

ω

)

showninEq.(12):

˜ ˜

We obtain the space-time domain response by applying the inverseFourier Transformover wavenumber kx and fre-quency

ω

toEq.(13): U₁

(

r₁ ,

θ

1 ,x,t

)

= ₄

_π

1₂  _∞ −∞ _∞ −∞ ˜ ˜ U₁

(

r₁ ,

θ

1 ,kx,

ω

)

exp

(

+i

(

ω

t− kxx

)

)

dkxd

ω

= P0 2

π

 ∞ −∞ ∞ −∞ ˜ ˜ U1 ,aux

(

r1 ,

θ

1 ,kx,

ω

)

δ

(

ω

−

− kxV

)

exp

(

+i

(

ω

t− kxx

)

)

dkxd

ω

= P0 2

π

 ∞ −∞ 1 V ˜ ˜ U1 ,aux



r1 ,

θ

1 ,

ω

−_V

,

ω

exp



− i

ω

−

V x

exp

(

i

ω

t

)

d

ω

. (14)

InEq.(14),theinverseFourierTransformoverwavenumberkx hasbeenevaluatedanalytically,whereastheinverseFourier Transformoverfrequency

ω

needstobeevaluatednumerically.

The radial displacement componentofthe steady-stateresponse atthe loading pointcan be obtainedby substituting x=Vt intoEq.(14): Ur1



R_i ,−

π

2,Vt,t



=P0exp

(

i

t

)

2

π

∞ −∞ 1 V ˜ ˜ Ur1,aux



R_i ,−

π

2,

ω

−

V ,

ω

d

ω

. (15)

InaccordancewithEq.(15),thecomplexamplitudeofthisharmonicresponse,whichisrelevantforthecomputationofthe equivalentstiffness,isgivenas

U0

(

,V

)

=₂P

_π

0 ∞ −∞ 1 V ˜ ˜ Ur1,aux



Ri ,−

π

₂,

ω

−_V

,

ω

d

ω

. (16)

UsingEq.(16),theequivalentstiffnessKeq definedinEq.(1)isobtained:

Keq = 1 1 2 π ∞ −∞ V1 U˜˜r1,aux



R_i ,−π₂ ,ω− V ,

ω



d

ω

. (17) Wenotethatthisresultissingle-valuedforall

astheGreen’sfunctions(seeSection3)usedintheindirectBEM compu-tationsareuniquelydefined.

Wealsoconsiderthesteady-stateresponseatafixedobservationpointx=0,whichisneededforthevalidationofthe indirectBEM.Substitutingx₌0intoEq.(14)givesthecorrespondingdisplacementvector:

U1

(

r1 ,

θ

1 ,0,t

)

=₂P

_π

0 ∞ −∞ 1 V ˜ ˜ U1 ,aux



r1 ,

θ

1 ,

ω

−_V

,

ω

exp

(

i

ω

t

)

d

ω

. (18)

Eq.(18)containstheresponsesobservedatx=0;foranobservationpointatthetunnelinvert,t<0indicatesthatx=0> Vt,whichmeansthatthemovingloadhasnotreachedtheobservationpointyet;t=0indicatesthatthemovingloadisat theobservationpoint;t>0indicatesthatx=0<Vt,whichmeansthatthemovingloadhaspassedtheobservationpoint. For thecaseof astationary (i.e.,non-moving)harmonic point load, whichis alsoused inSection 4 forvalidation,an expressionfortheinduceddisplacementsisgiveninAppendixA.

3. Indirectboundaryelementmethod

In this paper, the indirectBEM is employed to compute the response ofthe tunnel-soil systemin the wavenumber-frequency(r₁ ,

θ

₁ ,kx,

ω

)domain.Tothisend,theGreen’sfunctionsofthesoilandtunnelareneeded;notethattheindirect BEM uses the Green’s functionsof the soil withoutcavity (full-space or half-space).In Sections 3.1and3.2, theGreen’s functionsofthesoilandtunnelarepresented.TheindirectBEMisformulatedinSection3.3.

3.1. Green’sfunctionsofthesoil

Theso-calledtwo-and-a-halfdimensionalGreen’sfunctionsofanelastodynamicfull-space[30]andahalf-space[31]are usedinourwork.Thesourceconsideredinthementionedpapersisaspatiallyvaryinglineloadinthelongitudinal direc-tion, havingthe formof fj

(

y,z,x,t

)

=Fj

δ

(

y− ys

)

δ

(

z− zs

)

exp

(

i

(

ω

t− kxx

)

)

, where j=y,z,xindicates the directionofthe load,subscript“s” indicatesthecoordinatesofthesourcepoint,andFjistheamplitudeofthesource(Green’sfunctionscan beobtainedbysettingF_j=1).

The Green’s functionsofthe half-spaceconsist ofsourceterms whichare the sameasthose ofthefull-space, andof surfacetermswhicharenecessarytosatisfythestress-freeboundaryconditionsatthesurfaceofthehalf-space[31]. How-ever, we found that the stress-freeconditionsare not satisfied usingthe Green’s functionspresented in[31], whilethey are satisfiedwhenthe sourcetermsarereplaced bythe onespresentedin[30]that containsthefull-spaceGreen’s func-tions.Therefore,theGreen’sfunctionsofthehalf-spaceusedinthecurrentpaperconsistofthesourcetermspresentedin

Fig. 3. A cylindrical shell, the associated coordinate system and displacement components of the shell.

Becausethereferenceframesin[31]aredifferentfromthatinthecurrentpaper,wehavetotransformtheGreen’s func-tions fordisplacements andstressesgivenin[31]through thefollowingrelationsG˜˜_u,1 =T₁ G˜˜u,ref TT 1 andG˜˜σ,1 =[T1 ]G˜˜σ,ref TT 1 ,

wheresubscripts“1” and“ref” denotetheresponses definedinthecoordinatesystemsofthecurrentandreferencepapers, respectively.Thetransformationmatrixreads:

T₁ =



_{1 0 0} 0 −1 0 0 0 1



. (19) ˜ ˜

Gu,1 andG˜˜σ,1 are theGreen’s functionsfor displacements andstressesof the soil withouttunnel/cavity (i.e.,of the

full-space orhalf-space),andare3 × 3matrices. InmatrixG˜˜u,1 ,thefirst, secondandthird rowsrepresentthe displacement

components u˜˜r₁,u˜˜_θ1 andu˜˜x,whilethefirst, secondandthird columnscorrespondtothe spatiallyvarying unit lineloads actinginy,zandxdirections,respectively.InmatrixG˜˜_σ,1 ,thefirst,secondandthirdrowsrepresentthestresscomponents

˜ ˜

σ

r1r1,

σ

˜˜r1θ1 and

σ

˜˜r1x,whilethecolumnsalsocorrespondtotheloadsindifferentdirections.

3.2. Green’sfunctionsandresponsesofacylindricalshell

ThetunnelismodelledbyaninfinitelylongcylindricalFlüggeshell.TheassociatedcoordinatesystemisshowninFig.3. Theequationsofmotionoftheshellread[32]:

−ρ2 R



1−

ν

2 2



E2

∂

2 _u¯

∂

t2 −

ν

2

∂

w¯

∂

x2 − 1 R

∂

v

¯

∂θ

2 − 1 Ru¯− h2 12

R

∂

4 _u¯

∂

x4 ₂ + 2 R

∂

4 _u¯

∂

x2 ₂

∂θ

₂ 2 + 1 R3

∂

4 _u¯

∂θ

4 2



+ R



1−

ν

2 2



E2 h qr2 +h2 12

∂

3 _w¯

∂

x3 ₂ −

(

1−

ν

2

)

2R2

∂

3 _w¯

∂

x2

∂

θ

2 2 +

(

3−

ν

2

)

2R

∂

3

_v

_¯

∂

x2 ₂

∂

θ

2 − 1 R3 u¯− 2 R3

∂

2 _u¯

∂θ

2 2



=0, (20) −

ρ

2 R



1−

ν

2 2



E

∂

2

_v

_¯

∂

t2 +

(

1+

ν

2

)

2

∂

2 _w¯

∂

x2

∂

θ

2 +R

(

1−

ν

2

)

2

∂

2

_v

_¯

∂

x2 ₂ + 1 R

∂

2

_v

_¯

∂θ

2 2 +1 R

∂

u¯

∂θ

2 +R



1−

ν

2 2



E2 h qθ2 +h2 12

3

(

1−

ν

2

)

2R

∂

2

_v

_¯

∂

x2 ₂ −

(

3−

ν

2

)

2R

∂

3 _u¯

∂

x2 ₂

∂

θ

2



=0, (21) −

ρ

2 R



1−

ν

2 2



E2

∂

2 _w¯

∂

t2 +R

∂

2 _w¯

∂

x2 ₂ +

(

1−

ν

2

)

2R

∂

2 _w¯

∂θ

2 2 +

(

1+

ν

2

)

2

∂

2

_v

_¯

∂

x2

∂

θ

2 +

ν

2

∂

u¯

∂

x2 +R



1−

ν

2 2



E2 h qx2 +h2 12

(

1−

ν

2

)

2R3

∂

2 _w¯

∂θ

2 2 −

∂

3 u¯

∂

x3 ₂ +

(

1−

ν

2

)

2R2

∂

3 _u¯

∂

x2

∂

θ

₂ 2



=0, (22)

where u¯,

v

¯ andw¯ are the mid-surfacedisplacements in directionsr₂ ,

θ

₂ and x₂ , respectively.E₂ is theYoung’s modulus of theshell andhits thickness.Theradii ofthe inner andouter surfaceoftheshell canbe expressed asR_i =R−h

2 and Ro =R+h

2 ,respectively.qr₂,q_θ2 andqx₂ arethenetexternalstressesactingontheshell,namelythedifferencebetweenthe stressesactingattheinnerandoutersurfaces.

Thegoverningequationsoftheshellcanberewrittenintomatrixformas

A¯u₂ = ¯q₂ , (23)

where ¯u₂ =

(

u¯,

v

¯,w¯

)

isthedisplacementvector, ¯q₂ =



¯qr₂,¯q_θ2,¯qx₂



isthenetstressvectorcorrespondingtothemid-surface oftheshell,andAisanoperatormatrixgiveninAppendixB.

Thestressvector ¯q2isrelatedtothestressvectorsqo ₂ andqi 2correspondingtotheouterandinnersurfacesoftheshell,

respectively,throughthefollowingrelations[36]:

¯q2 =

⎛

⎝

1 h 2 R∂θ∂2 h 2 ∂∂x2 0 1+ h 2 R 0 0 0 1

⎞

⎠

qo ₂ =Bo qo ₂ , ¯q2 =

⎛

⎝

1 −h 2 R∂θ∂2 −h 2 ∂∂x2 0 1+−h 2 R 0 0 0 1

⎞

⎠

qi ₂ =Bi qi ₂ . (24)

AccordingtoLove’ssimplificationintheshelltheory[37],thelongitudinalandtangentialdisplacementsvarylinearlyacross the shell’sthickness, whereas theradial displacementisindependent ofradial coordinate.Therefore, themid-surface dis-placement vector ¯u₂ is relatedtothedisplacement vectoruo ₂ corresponding tothe outersurfaceofthe shellthroughthe followingrelation uo ₂ =

⎛

⎝

−h1 0 0 2 R∂θ∂2 1+ h 2 R 0 −h 2 ∂∂x2 0 1

⎞

⎠

¯u₂ =D¯u_{2 .} (25)

AfterapplyingtheFourierTransformovertimetandspatialcoordinatex₂ toEq.(23),computingtheFouriercoefficients of the circumferential harmonics(i.e., thesecond relationin Eq.(10)), andconsidering Eqs.(24) and(25),the governing equationsoftheshellcanbewrittenas

˜ ˜

AnD˜˜−1 n u˜˜2 o ,n=B˜˜o nq˜˜o 2 ,n+B˜˜ni q˜˜i 2 ,n, (26) whichisessentiallyasetofalgebraicequations;ndenotesthenumberofthecircumferentialharmonic,andmatricesA˜˜n,B˜˜o n,

˜ ˜

Bi _n andD˜˜naregiveninAppendixB;u˜˜o _{2 ,n}andq˜˜{₂ o _,n,i }containtheFouriercoefficientsofu˜˜o ₂ andq˜˜₂ {o ,i },respectively.Ifq{₂ o ,i } are takenasq₂ {o ,i }₌[

δ

(

θ

2

)

,

δ

(

θ

2

)

,

δ

(

θ

2

)

]T

δ

(

x₂

)

δ

(

t

)

_,then ˜q˜{₂ o _,n,i }₌



₂ 1 _π,2 1 _π,2 1 _π



T .Theassociated Green’sfunctions ofthe shellcan be derivedbysolvingEq.(26)foreachoftheloadcomponents,andsubsequentlyaddingthesolutionsforall components intheFourierseries(seeEq.(10)):

˜ ˜ go

(

θ

2 ,kx2 ,

ω

)

= n_= ∞ n= −∞ ˜ ˜ go _nexp

(

in

θ

2

)

= n_= ∞ n= −∞ 1 2

π

˜ ˜ DnA˜˜−1 n B˜˜o nexp

(

in

θ

2

)

, ˜ ˜ gi

(

θ

2 ,kx2 ,

ω

)

= n_= ∞ n= −∞ ˜ ˜ gi nexp

(

in

θ

2

)

= n_= ∞ n= −∞ 1 2

π

˜ ˜ DnA˜˜−1 n B˜˜i nexp

(

in

θ

2

)

, (27)

where g˜˜o _n interrelates u˜˜o _{2 ,n} andq˜˜o _{2 ,n}, _g˜_˜i

n interrelates u˜˜o 2 ,n andq˜˜2 i ,n, andg˜˜o ,g˜˜i , g˜˜o n andg˜˜i n are 3× 3 matrices. Thepositive directionsofthelongitudinalaxesintheglobalcoordinatesystem(seeFig.1)andthelocalcoordinatesystemfortheshell (seeFig.3)areoppositetoeachother(i.e.,x_{2 =−x}).Therefore,therelationbetweenthelongitudinalwavenumberskx2 and kx (which,forthemovingoscillatorypointloadconsideredinSection2.3,isdefinedas ω−_V )isasfollows:kx2 =−kx;this relationisusedbelow.

Usingtheconvolutionrule,thedisplacementvectoroftheshellunderanarbitraryloadcanbeobtainedas ˜ ˜ uo ₂

(

θ

2 ,kx2 ,

ω

)

= Ro  2 π 0 n_= ∞ n= −∞ 1 2

π

Ro ˜ ˜ DnA˜˜−1 n B˜˜o nexp



in



θ

2 −

θ

2 



_˜ ˜ qo ₂



θ

₂ ,kx2 ,

ω



d

θ

2  +R_i _{2 π} 0 n_= ∞ n= −∞ 1 2

π

R_i ˜ ˜ DnA˜˜−1 n B˜˜i nexp



in



θ

₂ −

θ

 2



_˜ ˜ qi ₂



θ

₂ ,kx2 ,

ω



d

θ

2 . (28)

The displacements in Eq. (28) are definedin the local coordinate system of the shell (r2,

θ

2, x2, see Fig. 3). To satisfy

the continuityofdisplacementsandstressesattheshell-soilinterface, thedisplacementandstressvectorsu˜˜o ₂ ,_q˜_˜o 2 andq˜˜i 2

definedinthe localcoordinatesystemoftheshell havetobe transformedto U˜˜₂ o ,Q˜˜o ₂ andQ˜˜i ₂ ,respectively, definedinthe globalcylindricalcoordinatesystemofthesoil(r₁ ,

θ

₁ ,x),whichhasoriginatthecenterofthetunnel(seeFig.1),through relations ˜ ˜ U₂ o

(

θ

1 ,kx,

ω

)

=T2 u˜˜o ₂

(

θ

2 ,−kx,

ω

)

, ˜ ˜ qo ₂

(

θ

2 ,−kx,

ω

)

=TT 2Q˜˜o 2

(

θ

1 ,kx,

ω

)

, q˜˜i 2

(

θ

2 ,−kx,

ω

)

=T2T Q˜˜i 2

(

θ

1 ,kx,

ω

)

, (29) inwhich

θ

₁ =

θ

2 and T2 =



_{1 0 0} 0 1 0 0 0 −1



. (30)

Fig. 4. The elastic half-space with the fictitious cavity. L r (solid line) and L s (dashed line) are surfaces at which the receiver (filled circles) and source (open

circles) points are located, respectively.

SubstitutingEq.(29)intoEq.(28),weobtainthedisplacementsoftheshellintheglobalcylindricalcoordinatesystem: ˜ ˜ Uo ₂

(

θ

1 ,kx,

ω

)

= Ro  _{2 π} 0 n_= ∞ n= −∞ 1 2

π

RoT2 ˜ ˜ DnA˜˜−1 n B˜˜no TT 2 exp



in



θ

1 −

θ

2 



Q˜˜o 2



θ

2 ,kx,

ω



d

θ

₂  +R_i  _{2 π} 0 n_= ∞ n= −∞ 1 2

π

R_i T2 ˜ ˜ DnA˜˜−1 n B˜˜ni TT 2 exp



in



θ

1 −

θ

2 



Q˜˜i 2



θ

2 ,kx,

ω



d

θ

₂ . (31)

3.3. Formulationoftheindirectboundaryelementmethod

Inthissection,theformulationoftheemployed indirectBEMispresented.Theexcitationstress vectorsare [

σ

˜˜aux ,0,0] and[

σ

˜aux ,0,0] forthemoving andstationarypoint loadcases,respectively,which comeintoplay through Q˜˜i ₂ (asshown below).Theexpressions for

σ

˜˜aux and

σ

˜aux are giveninEq.(12)andEq.(A.3),respectively.Fortheconsideredproblem,we havecontinuousdisplacementsandstressesatthetunnel-soilinterfaceL_r :

˜ ˜ U1

(

Ro ,

θ

1 ,kx,

ω

)

=U˜˜o ₂

(

θ

1 ,kx,

ω

)

, (32) ˜ ˜



1

(

Ro ,

θ

1 ,kx,

ω

)

=



˜˜o ₂

(

θ

1 ,kx,

ω

)

, (33)

where U˜˜₁ andU˜˜₂ are the displacementvectors of thesoil andtunnel, respectively, and



˜˜₁ and



˜˜₂ their stress vectors. The displacementandstress vectorsofthesoil andshellatthetunnel-soil interface areexpressed asU˜˜1=[U˜˜r1,U˜˜θ1,

˜ ˜ Ux]T , ˜ ˜



1 =[

σ

˜˜r₁r₁,

σ

˜˜r₁θ₁,

σ

˜˜r₁x]T , U˜˜o ₂ =[U˜˜ro ₂,U˜˜_θo 2, ˜ ˜

Uxo ₂]T and



˜˜₂ o =[

σ

˜˜ro ₂r₂,

σ

˜˜ro ₂θ₂,

σ

˜˜ro ₂x₂]T . Note that all thesedisplacements and stresses aredefinedintheglobalcylindricalcoordinatesystem(r₁ ,

θ

₁ ,x).

AccordingtotheindirectBEM,thedisplacementandstressvectorsinthesoilaregivenas[36]

˜ ˜ U1

(

xr ,kx,

ω

)

=  Ls ˜ ˜ Gu,1

(

xr ,xs ,kx,

ω

)

F

(

xs

)

dl

(

xs

)

, (34) ˜ ˜



1

(

xr ,kx,

ω

)

=  Ls ˜ ˜ G_σ,1

(

xr ,xs ,kx,

ω

)

F

(

xs

)

dl

(

xs

)

, (35) where F(xs ) is the yet unknown vector of the source amplitudes placed inside the fictitious cavity which is commonly usedfortheindirectBEM(seeFig.4).Vectorsxr =[xr ,yr ,zr ]andxs =[xs ,ys ,zs ]arecoordinatesofthereceiverandsource points, respectively. Ls isthe surfaceatwhich thesource pointsare located,andthe radius ofthesurfaceLs is takenas Rs =Ro − 3

(

2

π

Ro /Nr

)

,whereNr denotesthenumberofreceiverpoints,andNr ≥ 20assuggestedin[36].ThesurfaceLr at whichthereceiverpointsarelocatedliesatr₁ =Ro ,whichistheoutersurfaceoftheactualtunnel.

AnexpressionforthedisplacementvectoroftheshellhasbeenobtainedinEq.(31),andcanberewrittenas ˜ ˜ Uo ₂

(

xr,kx,

ω

)

=  Lr ˜ ˜ Go u,2



xr,x_r ,kx,

ω



Q˜˜o ₂



x_r ,kx,

ω



dl



x_r



+  L ˜ ˜ Gi u,2



xr,x,kx,

ω



Q˜˜i ₂



x,kx,

ω



dl



x



, (36) whereG˜˜{_uo _,2 ,i }containtheGreen’sfunctionsforthedisplacementsoftheshelldefinedintheglobalcylindricalreferenceframe. TheintegrationsurfacesinEq.(36)areL_r andL(i.e.,sourcepointslocatedatr_{1 =}R_o andr_{1 =}R_{i ,}respectively)corresponding totheGreen’sfunctionsoftheshellrelatedtotheforcesactingatitsouterandinnersurfaces.Thestressvectorsactingat theouterandinnersurfacesoftheshell,respectively,forthecurrentproblemread,employingEq.(33)

˜ ˜ Qo ₂



x_r ,kx,

ω



=

_

˜˜₂



_x r ,kx,

ω



=

_

˜˜₁



_x r ,kx,

ω



, _Q˜˜i 2



x,kx,

ω



=_P˜˜



_x_,kx,

_ω



_{, (37)}

Table 1

Displacement components (20 ·log 10 | U i |) at different locations ( r 1 , θ1 , x ) for a tunnel embedded in an elastic

full-space subjected to a stationary harmonic point load with excitation frequency f 0 = 10 Hz obtained using

different numbers of source and receiver points ( N s , N r ). In each row, the displacements are normalised by the

corresponding response obtained using (Ns , N r) = (20 , 40) .

Displacements (Ns , N r) = (20 , 40) (Ns , N r) = (30 , 60) (Ns , N r) = (40 , 80) (Ns , N r) = (60 , 60) Ur1 ( R o , −π2 , 0) 1.0000 1.0002 1.0002 1.0002 Ur1 ( R o , π2 , 0) 1.0000 1.0008 1.0008 1.0008 Ur1 ( R o , π, 0) 1.0000 1.0000 1.0000 1.0000 Uy1 (20 m, π, 20 m) 1.0000 1.0000 1.0000 1.0000 Uz1 (20 m, π, 20 m) 1.0000 1.0000 1.0000 1.0000 Ux (20 m, π, 20 m) 1.0000 1.0000 1.0000 1.0000

whereP˜˜ istheexcitationstressvector(i.e.,[

σ

˜˜aux ,0,0]T or[

σ

˜aux ,0,0]T )inducedbytheexternalpointload(seeFig.2(a)). SubstitutingEqs.(34)and(36)intoEq.(32),andconsideringEqs.(35)and(37),wederivetheboundaryintegralequation intermsoftheunknownsourceamplitudevector:

 Ls ˜ ˜ K

(

xr ,xs ,kx,

ω

)

F

(

xs

)

dl

(

xs

)

=R˜˜

(

xr ,kx,

ω

)

, (38) where ˜ ˜ K

(

xr ,xs ,kx,

ω

)

=G˜˜u,1

(

xr ,xs ,kx,

ω

)

−  Lr ˜ ˜ Go _u,2



xr ,xr ,kx,

ω



_˜˜ G_σ,1



x_r ,xs ,kx,

ω



dl



x_r



, ˜ ˜ R

(

xr ,kx,

ω

)

=  L ˜ ˜ Gi _u,2



xr ,x,kx,

ω



_˜˜ P



x,kx,

ω



dl



x



. (39)

Inordertocompute thesourcevectorforeverykxand

ω

combination,Eq.(38)isdiscretised(i.e.,surfacesLr ,LandLs ),as indicatedabove.NotethatthesizeofthematrixrelatedtoK˜˜is(3Nr × 3Ns ),whereNs denotesthenumberofsourcepoints; thisimpliesthat itisasmallmatrix,andtheinversionofthematrixdoesnottake muchtime.Themosttimeconsuming part(even intheentireprocedureofinstabilityanalysis)istheevaluationofintegralsoverthehorizontalwavenumberky intheGreen’sfunctionsofthehalf-space,forwhichweusethe“quadv” routineinMatlab.

4. Validations

TovalidatetheaccuracyofthepresentedindirectBEM,wecomparetheresultsobtainedbytheproposedmethodwith thosecalculatedby Yuanetal.[38].Thefirstvalidationisperformedforthecaseofatunnelembeddedinanelastic full-space subjectto a stationary harmonicpoint load atthe tunnel invert. The excitation for theindirect BEM computation in thiscaseis

σ

˜aux (Eq.(A.3)), withP0=1N, andthe steady-stateresponseis giveninEq.(A.5).Theelastic full-spaceis

characterisedby its longitudinalwave speedC_{P ,1} =944 m/s,shearwave speedC_{S ,1} =309m/s,density

ρ

₁ =2000kg/m3 andmaterialdampingratio

ξ

₁ =0.03.TheelasticparametersforthetunnelaretheYoung’smodulusE₂ =50GPa,Poisson’s ratio

ν

₂ =0.3, density

ρ

₂ =2500 kg/m3 and material dampingratio

ξ

₂ =0.The inner andouter radii of the tunnel are R_{i =}2_.75mandRo=3m.

Before showing the results, we first present a convergence test forthe proposed method for the considered loading case. AsshowninEq.(A.5),theinverseFourierTransformoverlongitudinalwavenumberkx hastobeevaluatedtogetthe harmonicresponseinthespace-timedomain.Theintegral wascomputednumericallyusinganinversefastFourier Trans-formalgorithminMatlab.Theconvergencewastestedregardingthediscretisationofkx (i.e.,

kx andkmax x ),themaximum numberofcircumferentialmodesoftheshellN_shell max inEq.(31)andthemaximumnumberofFouriercomponentsN_load max in

Eq.(A.3),andthenumberofsourceandreceiverpoints(Ns ,Nr ).Wefoundthatitissufficienttousekmax x =2

π

,

kx=₁₀₂₃ 2 π , Nmax _{shell =}20andN_load max ₌20.Theconvergencetestforthenumberofsourceandreceiverpoints(N_s ,N_o )atdifferentlocations (r₁ ,

θ

₁ ,x)isgiveninTable1.Responsesatthetunnelinvert,tunnelapex(Ro , π₂ ,0),tunnelside(Ro ,

π

,0) andatapoint farfromtheload(20m,

π

,20m)arepresented;theloadischaracterisedbyP₀ =1Nand f₀ = _{2 π} =10Hz.Itisclearthat convergedresultscanindeedbeobtainedusing(N_s ,N_r )₌(20,40).

Fig. 5 showsthe converged vertical displacements atthe tunnel invert, tunnel apex andtunnel side as a function of frequencyforthefirstvalidationcase.Agoodagreementcanbeobservedbetweentheresultsobtainedbydifferentmethods, whichvalidatestheproposedmethodanditsimplementation.

Thesecondvalidationcaseisthatofashallowtunnelembeddedinanelastichalf-spacesubjecttoastationaryharmonic load. The soil is characterisedby its longitudinal wave speed C_{P ,1} =400 m/s, shear wave speed C_{S ,1} =200 m/s, density

Fig. 5. Vertical displacements ( 20 · log 10| U z1| ) at locations of (a) tunnel invert ( r 1 = R o , θ1 = −π2 , x = 0 ), (b) tunnel apex ( r 1 = R o , θ1 = π2 , x = 0 ) and (c)

tunnel side ( r 1 = R o , θ1 = π, x = 0 ) for the case of a tunnel embedded in an elastic full-space subject to a stationary harmonic point load.

Fig. 6. Displacement components (20 · log 10 | U i |) at a point on the ground surface ( y = −20 m , z = 0 , x = 20 m ) for the case of a tunnel embedded in an elastic half-space ( H = 5 m ) subject to a stationary harmonic point load acting at the tunnel invert: (a) horizontal displacement, (b) vertical displacement and (c) longitudinal displacement.

Table 2

Velocities ( V i ( y, z, x )) and displacement (Ur1(r1 , θ1 , x )) at different locations for a tunnel embedded in an elastic half- space subject to a uniformly moving point load ( V = 75 m / s , f 0 = 0 ) using different numbers of source and receiver

points ( N s , N r ). In each row, the responses are normalised by the corresponding response obtained using (Ns , N r) =

(20 , 40) . Responses Time (s) (Ns , N r) = (20 , 40) (Ns , N r) = (30 , 60) (Ns , N r) = (40 , 80) (Ns , N r) = (60 , 60) Vz t = 0 1.0000 1.0000 1.0000 1.0000 (0, 0, 0) t = 1 1.0000 1.0000 1.0000 1.0000 Vx t = 0 1.0000 1.0000 1.0000 1.0000 (0, 0, 0) t = 1 1.0000 1.0000 1.0000 1.0000 Vz t = 0 1.0000 1.0000 1.0000 1.0000 ( −20 m , 0, 0) t = 1 1.0000 1.0000 1.0000 1.0000 Vx t = 0 1.0000 1.0000 1.0000 1.0000 ( −20 m , 0, 0) t = 1 1.0000 1.0000 1.0000 1.0000 Ur1 t = 0 1.0000 0.9998 0.9998 0.9998 ( R o , −π2 , 0) t = 1 1.0000 1.0000 1.0000 1.0000

except that

ξ

₂ =0.015, and the burial depth of the tunnel is H=5 m. It is notedthat in the referencepaper [38], the mentionedmaterialdampingratioshouldbethelossfactor(thereisadifferenceofafactor2),whichisindicatedinpaper

[39].Thisalsoholdsforthenextvalidationcase.Thedisplacementcomponents(againobtainedusingEq.(A.5))atapointon thegroundsurface(y=−20m,z=0,x=20m)arepresentedinFig.6,whereagainagoodmatchbetweentheresultsis observed.Theminordifferencescanbeattributedtotheuseofacontinuumtomodelthetunnelin[38],insteadofashell. The third validation comprisesthe caseof atunnel embeddedin an elastic half-spacesubjectto a uniformlymoving constant point load. The excitation forthe indirect BEM computationin thiscaseis

σ

˜˜_aux (Eq.(12)) andthe steady-state response isgiveninEq.(14).The parametersforthesoilandtunnelare asfollows:

μ

₁ =1.154× 107 _N/m2 _,

_λ

1 =1.731×

107 N/m2 ,

ρ

1 =1900 kg/m3 ,

ξ

1 =0.025,

μ

2 =1.042× 1010 N/m2 ,

λ

2 =6.944× 109 N/m2 ,

ρ

2 =2400 kg/m3 ,

ξ

2 =0.01, R_{i =}2_.75 m_, R_{o =}3 m_,H₌15 m.The moving loadis characterised by velocityV₌75 m_/s_, excitation frequency f_{0 =}0 andP₀ =1N.TheinverseFourierTransformoverfrequenciesneedstobeevaluatedtogetthespace-timedomainresponse (see Eq.(14)).The convergenceforthemoving point loadcasewastestedregardingthediscretisation of

ω

(i.e.,

ω

and

ω

max _{), N}max

shell in Eq.(31) and Nload max in Eq.(12),and the numberof source andreceiver points (Ns , Nr ). Numerical results

relatedtotwopointsonthegroundsurfaceandoneatthetunnelinvertarepresentedinTable2fortheconsideredcase ofthemovingloadusingdifferentnumbersofsourceandreceiverpoints(Ns ,Nr ).Wefoundthatconvergedresultscanbe obtainedusing fmax = ω2maxπ =15Hz,

f= ω2π =0.05Hz,Nmax shell =20,Nmax load =20and(Ns ,Nr )=(20,40).Thisisclearfrom

Fig. 7. Velocities ( V i ) at the origin ( y = 0 , z = 0 , x = 0 ) on the ground surface for the case of a tunnel embedded in an elastic half-space ( H = 15 m ) subject to a uniformly moving point load ( V = 75 m / s , f 0 = 0 ): (a) vertical velocity and (b) longitudinal velocity.

Table2whichpresentstheresponsesobservedatx₌0forvaryingtimemoments:t₌0meansthattheloadisrightbelow theobservationpoint,whereast=1sindicatesthattheloadhaspassedthatpoint.Fig.7presentsthecomparisonbetween theresultsobtainedbytheproposedmethodandthoseshownintheliterature.Thegoodagreementgivesconfidenceabout theaccuracyoftheproposedmethod.

The convergencerequirementsforthe computationofthe equivalentdynamicstiffnessfordifferentvelocities are pre-sentedinAppendixC.

5. Instabilityofvibrations

ThemainframeworktoconductinstabilityanalysishasbeengiveninSection2.2.Inthecurrentsection,wepresentthe resultsforboththefull-spaceandhalf-space.Thebase-caseparametersofthetunnel-soilsystemarelistedinTable3;itis notedthatthebasecaseassumesthattheburial depthH→∞.Theparameters presentedinTable3representasoftsoil andaconcretetunnel.

We chosethefull-spacecaseasthebasecasesimplybecausetheresultsforthehalf-spaceare prettysimilar, andthe computationofthe two-and-a-halfdimensionalGreen’s functionsofthefull-space[30] islessexpensivethan thatofthe half-spaceGreen’sfunctions[31].TheGreen’sfunctionsofthefull-spaceareavailableanalyticallyandcanbeevaluatedvery fast;however,thesurface-termspartoftheGreen’sfunctionsofthehalf-space(seeSection3.1)arenotavailableanalytically, andintegralsoverthehorizontalwavenumberkyneedtobeevaluatednumerically.

Table 3

Base-case parameters of the tunnel-soil system.

Soil Tunnel μ1 = 1 . 154 × 10 7 N / m 2 μ2 = 1 . 042 × 10 10 N / m 2 λ1 = 1 . 731 × 10 7 N / m 2 λ2 = 6 . 944 × 10 9 N / m 2 ρ1 = 1900 kg / m 3 ρ2 = 2400 kg / m 3 ξ1 = 0 . 05 ξ2 = 0 . 02 CS,1 = 77 . 94 m / s CS,2 = 2083 . 7 m / s CP,1 = 145 . 80 m / s CP,2 = 3402 . 5 m / s CR,1 = 72 . 29 m / s Ri = 2 . 75 m , R o = 3 m , h = 0 . 25 m , H → ∞

5.1. Criticalvelocityforinstabilityofthemovingobject

As has been discussedin Section 2.2, theimaginary part ofthe equivalentstiffness beingnegativeindicates that the vibrationoftheobject(i.e.,movingmassoroscillator)canbecomeunstable.Therefore,wefirststudytheequivalentstiffness tofindthecriticalvelocityforinstabilityofthemovingobject(heredefinedasthevelocityatwhichIm(K_eq )_<0firsttakes place).Notethatthecriticalvelocityforinstabilitygenerallydiffersfromtheclassical criticalvelocityatwhichthe steady-stateresponseinducedbyamovingloadisextreme(i.e.,resonance).Inthegeneralcase,wheretheoscillatorhasdissipative components,thecriticalvelocityforinstabilityshouldbeidentifiedfromtheD-decompositioncurveinthecomplexMorK plane(seeSections2.2and5.2),notfromtheanalysisofIm(Keq )alone[19,26].Im(Keq )<0isonlyanecessaryconditionfor instability.Asthemovingoscillatorandmovingmassconsideredinthispaperdonothaveintrinsicdissipativecomponents, theircriticalvelocitiesforinstabilityarethesame.

Fig. 8. The real part of the equivalent stiffness for different velocities: (a) V = 0 . 90 V inst

cr , (b) V = 1 . 00 V crinst , (c) V = 1 . 01 V crinst , (d) V = 1 . 02 V crinst , (e) V =

1 . 03 V inst

cr , (f) V = 1 . 04 V crinst , (g) V = 1 . 05 V crinst , (h) V = 1 . 06 V crinst and (i) V = 1 . 20 V crinst . These results are related to the base case presented in Table 3 , and

Vinst

cr = 942 m / s .

In previous studies [11,19] where beamon elastic foundation models (without damping) are considered, it is shown that the critical velocity for instability V_cr inst of the moving object is equal to the critical velocity for resonance (of the undamped system),which inturnisequal tothe minimumphasevelocityV_ph min ofwavesinthe system.Fora half-space modelwitha regulartrackontop[25,26],thecriticalvelocityforinstabilityinthepresence ofdampingisslightlylarger thanVmin

ph .Additionally,Vph min,whichis closetoandsmallerthan thevelocity ofRayleighwaves,iseasily found fromthe

dispersionrelationofthesystem[12,40].However,forthetunnel-soilsystemconsideredinthispaper,itisverydifficultto getthedispersioncurves,asthedispersioncharacteristicsofthesystemareconsiderablymorecomplicated.Therefore,the minimum phasevelocity cannot be easily computed.We can, however,compute thesteady-stateresponse ofthe tunnel-soil systemsubjecttoa uniformlymovingnon-oscillatoryloadandcheck thefeaturesofresponses fordifferentvelocities to determine thecritical velocity forresonance(for the systemwithdamping,strictly speaking, buttheinfluence ofthe dampingonV_cr res issmall).ThisanalysisispresentedinAppendixC,anditshowsthatV_cr res ≈ 70m/sforthecurrent tunnel-soilsystem,whichisalsoclosetoandsmallerthanthevelocityofRayleighwaves,likefortheabove-mentionedhalf-space model.

Forthesystemwiththebase-caseparameters,we findthat theimaginarypartoftheequivalentstiffnessstartshaving a negative signfor at least a small frequency range ata velocity of Vinst

cr =942 m/s. Based on this critical velocity for

instability,westudythebehaviorofKeq (

,V)fordifferentvelocitiesintherangeof

(

0.9− 1.2

)

V_cr inst .Therealandimaginary parts are shownin Figs.8 and9,respectively. Nine differentvelocitieswere chosen toshow the features ofRe(Keq ) and Im(K_eq )asafunctionoftheloadfrequency

.Notethat,ifits massisrelativelysmall,thevibrationoftheobjectcanstill bestablewhenitmovesfasterthanV_cr inst ,aswillbedemonstratedinthenextsection(seealsoFig.11).

In Fig.8, we observethat the realpart ofthe equivalentstiffnessis positive, andthe decayingtrend ofRe(Keq ) with frequency

issimilarforeachvelocity.Thedecayingtrendmayberelatedtotheeffectofinertia.Thetroughcanprobably