翟朝嬌，夏唐代，陳煒昀，杜國(guó)慶

（1.浙江大學(xué)軟弱土與環(huán)境土工教育部重點(diǎn)實(shí)驗(yàn)室，浙江 杭州310058；2.浙江大學(xué)巖土工程研究所，浙江 杭州310058；3.南京工業(yè)大學(xué)巖土工程研究所，江蘇 南京210000；4.安徽省電力設(shè)計(jì)院，安徽 合肥230601）

深埋洞室常用于地下隧道、流體運(yùn)輸。隨著地下交通的發(fā)展，地下洞室在外力作用下的動(dòng)力響應(yīng)成為諸多領(lǐng)域的重要課題之一。相關(guān)研究者［1－5］利用波函數(shù)展開(kāi)法、積分變換和數(shù)值方法等研究了圓柱形洞室對(duì)波的散射和衍射。然而，洞室表面受到突加荷載作用而引起的波的傳播問(wèn)題越來(lái)越引起注意。G.Eason［6－7］、張慶元等［8］和劉國(guó)利等［9］運(yùn)用 Laplace變換及其逆變換，得到了球形空腔及圓柱形洞室表面受均勻且僅為時(shí)間的函數(shù)的外力作用所引起的彈性波的傳播規(guī)律及動(dòng)力響應(yīng)；H.Herman等［10］、T.L.Geers［11］和楊俊等［12］基于積分變換研究了無(wú)限長(zhǎng)圓柱形殼體在突加壓力作用下的瞬態(tài)動(dòng)力響應(yīng)；M.J.Forrestal等［13］求得了內(nèi)源Heaviside荷載作用下，以初等函數(shù)形式表示的流體介質(zhì)中無(wú)限長(zhǎng)彈性柱殼動(dòng)力響應(yīng)的精確解；T.B.Moodie等［14－15］運(yùn)用積分變換及數(shù)值方法，研究了圓柱形洞室表面受突加外荷載作用時(shí)膨脹波的傳播，得到了漸進(jìn)的波震面展開(kāi)式；N.Akkas等［16］，U.Zakout等［17］運(yùn)用殘余變量法將二階擴(kuò)散方程降為一階，得到了Heaviside荷載作用下，不同模態(tài)下圓柱孔洞的表面位移；V.R.Feldgun等［18］采用變分差分的方法研究了內(nèi)爆炸荷載作用下球形和圓柱形孔洞動(dòng)力響應(yīng)問(wèn)題，分析了彈性波在殼體和土體中的傳播；高盟等［19－20］和蔡袁強(qiáng)等［21］采用拉普拉斯變換，得到圓柱形孔洞在內(nèi)源荷載作用下的動(dòng)力響應(yīng)問(wèn)題，并運(yùn)用逆拉普拉斯變換的數(shù)值方法，給出了問(wèn)題的數(shù)值解。

上述研究均是圓柱形或球形洞室表面受到突加的徑向外力作用而引起的膨脹波的傳播，對(duì)洞室表面受軸向力作用下引起的剪切波的傳播的影響比較少。J.B.Haddow等［22－23］用數(shù)值的方法研究了非線性軸向剪切波的傳播，分析了洞室表面剪切應(yīng)力引起的有限振幅方位剪切波的傳播；K.Watanabe等［24］研究了圓柱形各向異性固體在點(diǎn)源SH波入射下的瞬態(tài)響應(yīng)，得到了任意圓柱形各向異性體的一般解決辦法；D.W.Barclay［25－27］采用波震面展開(kāi)的方法研究了圓柱形洞室表面在軸向剪切應(yīng)力的作用下軸向剪切波傳播，得到了近似的解析解答。由于數(shù)學(xué)分析上的困難，對(duì)剪切波傳播的研究都是基于數(shù)值方法或是半解析的，但是解析解答仍然非常必要。

本文中，將列車(chē)急剎車(chē)時(shí)對(duì)隧道產(chǎn)生的沖擊簡(jiǎn)化為無(wú)限彈性體中施加在圓柱形洞室表面沿軸線方向的均布荷載，利用Laplace變換，討論無(wú)限彈性土體內(nèi)圓柱形孔洞在軸線方向沖擊荷載的作用下的瞬態(tài)響應(yīng)，采用圍道積分［1］進(jìn)行Laplace逆變換，求解土體位移和應(yīng)力的一般解析表達(dá)式，并將計(jì)算結(jié)果與靜力情況下的響應(yīng)和采用F.Durbin［28］提出的拉普拉斯數(shù)值逆變換得出的結(jié)果進(jìn)行比較。

1 計(jì)算模型與求解

將土體視為無(wú)限彈性勻質(zhì)各向同性材料。地下洞室受到剎車(chē)荷載力，荷載先傳遞給襯砌，再由襯砌傳遞給外部土體，襯砌的模量遠(yuǎn)大于土體的模量，因此可以認(rèn)為沖擊通過(guò)襯砌均勻地傳遞給外部土體，列車(chē)急剎車(chē)時(shí)對(duì)洞室產(chǎn)生的沖擊近似沿軸線方向的均布荷載。該問(wèn)題簡(jiǎn)化為無(wú)限空間彈性土體中一半徑為r0的無(wú)限長(zhǎng)圓柱形孔洞內(nèi)部有一沿洞室軸線方向的均布突加荷載τi（t），如圖1所示。

圖1 計(jì)算模型Fig.1Computation model

1.1 控制方程

作為反平面問(wèn)題，以位移表示的控制方程為：

式中：r為土體中一點(diǎn)到洞室圓心的距離，w（r，t）為土體的軸線方向的位移，vs為剪切波在彈性土體中的傳播速度，vs＝（G/ρs）1/2，G 為材料的Lame常數(shù)，ρs為土體的密度。

1.2 邊界條件和初始條件

該問(wèn)題中，應(yīng)滿足在無(wú)窮遠(yuǎn)處應(yīng)力為0，因此應(yīng)力τ滿足邊界條件：

初始時(shí)刻，土體的位移和速度均應(yīng)為0，因此有初始條件：

1.3 求解頻域解

式（1）中同時(shí)含有空間變量r和時(shí)間變量t，邊界條件（3）中包含時(shí)間函數(shù)τi（t），在物理平面上不易求解。利用Laplace變換法，可將式（1）寫(xiě)為：

式中：W（r，s）為w（r，t）關(guān)于t的Laplace變換，記為W（r，s）＝L［w（r，t）］，s為變換參數(shù)，s＝a＋i b，a為任意正值，且a＞Re W。式（4）為零階虛宗量的Bessel方程，令k＝s/vs，則式（4）的通解［29］為：

式中：I0（kr）、K0（kr）分別為第一類(lèi)和第二類(lèi)零階虛宗量 Bessel函數(shù)，C1（s）、C2（s）為待定系數(shù)，由虛宗量Bessel函數(shù)的漸近性質(zhì)與遞推關(guān)系［29］，可知C1（s）＝0。

以 Heaviside荷載為例，τi（t）＝τ0H（t）（其中 H（x）為 Heaviside函數(shù)），L［τi（t）］＝τ0/s。經(jīng) La －place變換后，結(jié)合應(yīng)力位移關(guān)系，式（2）變?yōu)椋?/p>

結(jié)合（5）式可得，未知系數(shù)C2（s）＝τ0/［ksGK1（kr0）］，K1（kr0）為第二類(lèi)一階虛宗量Bessel函數(shù)。

由以上得到無(wú)限彈性土體中圓柱形孔洞反平面內(nèi)源問(wèn)題的頻域表達(dá)式：

式中：Τ（r，s）＝L［τ（r，t）］。

1.4 解析法求解時(shí)域解

為了在物理空間內(nèi)獲得時(shí)域解，須對(duì)頻域解施以Laplace逆變換，記：

式中：ζ＝kr0，γ＝vst/r0。為了計(jì)算式（7）的 Laplace逆變換，要選擇合適的積分路線?？疾霥（ζ）和E（ζ）的奇異性可知，ζ＝0是D（ζ）和E（ζ）的支點(diǎn)，另外I0（kr）和K0（kr）沒(méi)有零點(diǎn)，因此采用圍道積分的方法計(jì)算式（7）的積分。圖2給出了所選取的積分圍道，其中L1和L5為半徑無(wú)窮大的圓弧，L3為半徑趨于0的圓，L2和L4為分支割線。

圖2 計(jì)算模型Fig.2Computation model

在積分圍道內(nèi)部，D（ζ）和E（ζ）是解析的，由留數(shù)定理知：

可以證明［1］：

沿L3計(jì)算積分。ζ＝ρeiθ。應(yīng)用虛宗量Bessel函數(shù)的近似公式［30］，當(dāng)ρ→0時(shí)，K0（ζ）～－ln（ζ/2），K1（ζ）～1/ζ。根據(jù)留數(shù)定理［29］，D（ζ）和E（ζ）沿L3的積分可以寫(xiě)成：

計(jì)算D（ζ）和E（ζ）沿分支割線L2和L4的積分。沿L2和L4，ζ＝ρe±iπ＝－ξ，對(duì)于ρ＞0，Kν（－ξ）＝（－1）－νKν（ξ）?iπIν（ξ）［1］，其中Kν（ξ）和Iν（ξ）為ν階虛宗量貝賽爾函數(shù)。由此得到：

根據(jù)式（9）～（12）可以得到時(shí)域內(nèi)土體位移和應(yīng)力的解析表達(dá)式：

2 數(shù)值計(jì)算

用MATLAB計(jì)算工具可得到數(shù)值解，計(jì)算參數(shù)如下：孔洞半徑r0＝6m，τ0＝100kPa，G＝76 MPa，ρs＝1 900kg/m3。為方便討論計(jì)算結(jié)果，將位移、應(yīng)力和時(shí)間均進(jìn)行量綱一化處理：

圖3為r＝10，20，30m處土體的應(yīng)力時(shí)程曲線。在反平面階躍荷載作用下，波的起跳時(shí)間為波由波源傳播至該點(diǎn)所需要的時(shí)間，當(dāng)波傳播到該點(diǎn)，此處土體的應(yīng)力和位移瞬間增大并達(dá)到最大值，應(yīng)力波離開(kāi)后，該處應(yīng)力逐漸減小并逐漸趨于穩(wěn)定值，該穩(wěn)定值與靜力計(jì)算結(jié)果相吻合。3種不同半徑處應(yīng)力波起跳時(shí)間間隔為彈性波在此段距離傳播需要的時(shí)間。由于應(yīng)力波的發(fā)散傳播，距離波源越遠(yuǎn)，應(yīng)力值越小，這符合彈性波的性質(zhì)。

圖4為r＝10，20，30m處土體位移時(shí)程曲線。在波到達(dá)之前，該點(diǎn)的位移為零，當(dāng)波到達(dá)時(shí)，此處位移瞬間增大到最大值，隨后逐漸減小并趨向于靜力響應(yīng)狀態(tài)下的位移而保持不變。在r＝10，20，30 m 處最大量綱一位移分別為0.125 04、0.081 80和0.064 90，相對(duì)應(yīng)的量綱一時(shí)間為1.63、3.23和4.90，由此同樣可以得出與圖3相同的結(jié)論：3種不同半徑處位移達(dá)到最大值的時(shí)間間隔為彈性波在此段距離傳播需要的時(shí)間。結(jié)合圖3可知，應(yīng)力和位移的變化是同步的，這也符合彈性波動(dòng)理論中彈性波的傳播性質(zhì)。

圖3 應(yīng)力時(shí)程曲線Fig.3Stress －time curves

圖4 位移時(shí)程曲線Fig.4Displacement －time curves

由圖3～4還可以看出，土體應(yīng)力和位移最后的穩(wěn)定值分別為被積函數(shù)D（ζ）和E（ζ）沿L3的積分值，早期解的貢獻(xiàn)主要是沿分支割線的積分，這與文獻(xiàn)［1］中關(guān)于壓縮波入射的結(jié)論一致。

圖5和圖6分別為土體位移和應(yīng)力的最大值以及靜力值隨半徑的變化規(guī)律。由圖5～6可以看出，土體所受動(dòng)應(yīng)力和動(dòng)位移的最大值隨距離波源距離的增大呈現(xiàn)衰減趨勢(shì)，這是由于彈性波的發(fā)散傳播的特性，其趨勢(shì)與靜力下的衰減趨勢(shì)相近，且相同情況下動(dòng)力響應(yīng)的衰減更慢，隨著傳播距離的增大，動(dòng)力值與靜力值之比越來(lái)越大，因此在實(shí)際工程中動(dòng)荷載的破壞性要大得多。

圖5 應(yīng)力隨傳播距離的衰減曲線Fig.5Stress attenuation against propagation distance

圖6 位移隨傳播距離的衰減曲線Fig.6Displacement attenuation against propagation distance

3 計(jì)算結(jié)果驗(yàn)證

為了驗(yàn)證本文中計(jì)算結(jié)果的正確性，用本文中的計(jì)算方法和F.Durbin［28］提出的拉普拉斯數(shù)值逆變換分別計(jì)算r＝10m處的土體的應(yīng)力和位移，并將2種方法得出的結(jié)果做了比較，，如圖7～8所示，可以看出兩者吻合較好。

圖7 應(yīng)力時(shí)程曲線比較Fig.7Comparison of stress －time curves

圖8 位移時(shí)程曲線比較Fig.8Comparison of displacement －time curves

4 結(jié) 論

利用Laplace變換和圍道積分，討論了土體內(nèi)圓柱形孔洞的瞬態(tài)響應(yīng)問(wèn)題，求得了問(wèn)題的解析表達(dá)式。可以得出以下結(jié)論：（1）通過(guò)將本文的計(jì)算結(jié)果與用F.Durbin［28］提出的拉普拉斯數(shù)值逆變換的計(jì)算結(jié)果和靜力結(jié)果相比較，驗(yàn)證了本文方法在研究洞室在反平面沖擊荷載作用下的瞬態(tài)響應(yīng)時(shí)的適用性和正確性。（2）土體應(yīng)力和位移最后的穩(wěn)定值由被積函數(shù)沿原點(diǎn)附近小圓的積分值貢獻(xiàn)，早期解的貢獻(xiàn)主要是沿分支割線的積分，即初始時(shí)刻，分支割線上各點(diǎn)的貢獻(xiàn)最大。（3）在波到達(dá)之前，該處土體應(yīng)力和位移都保持為零，波到達(dá)之后，應(yīng)力和位移瞬間增大到最大值，接著慢慢減小并逐漸趨于穩(wěn)定值，不同半徑處應(yīng)力和位移到達(dá)峰值的間隔即為彈性波在此段距離傳播需要的時(shí)間。由于波沿徑向發(fā)散傳播，距離波源越遠(yuǎn)，應(yīng)力和位移值越小，相同的半徑處應(yīng)力和位移值相同?？锥吹脑摲€(wěn)定值與靜力下的受力相接近，且應(yīng)力和位移的變化是同步的。（4）在反平面動(dòng)荷載作用下，土體動(dòng)應(yīng)力和動(dòng)位移的最大值隨距離波源距離的增大呈現(xiàn)衰減趨勢(shì)，其趨勢(shì)與靜力值的衰減相近，且相同情況下動(dòng)力響應(yīng)的衰減更慢，因此在實(shí)際工程中危害性更大。

［1］ Pao Yih －h(huán)sing，Mow C C.Diffraction of elastic waves and dynamic stress concentrations［M］.New York，USA：Adam Hilger Limited，1973.

［2］ Lee J，Mal A K.A volume integral equation technique for multiple scattering problems in elastodynamics［J］.Applied Mathematics and Computation，1995，67（1/2/3）：135 －159.

［3］ Davis C A，Lee V W，Bardet J P.Transverse response of underground cavities and pipes to incident SV waves［J］.Earthquake Engineering and Structural Dynamics，2001，30（3）：395 －410.

［4］ Iakovlev S.Interaction of a spherical shock wave and a submerged fluid －filled circular cylindrical shell［J］.Journal of Sound and Vibration，2002，255（4）：615 －633.

［5］ Manolis G D.Elastic wave scattering around cavities in inhomogeneous continua by the BEM［J］.Journal of Sound and Vibration，2003，266（2）：281 －305.

［6］ Eason G.Propagation of waves from spherical and cylindrical cavities［J］.The Journal of Applied Mathematics and Physics，1963，14（1）：12 －23.

［7］ Eason G.The propagation of waves from a cylindrical cavity［J］.Journal of Composite Materials，1973，7（1）：90 －99.

［8］ 張慶元，戰(zhàn)人瑞.爆轟載荷作用下球形空腔的動(dòng)力響應(yīng)［J］.爆炸與沖擊，1994，14（2）：182 －185.Zhang Qing－yuan，Zhan Ren －rui.Dynamic response of a spherical cavity subjected to blast loads［J］.Explosion and Shock Waves，1994，14（2）：182 －185.

［9］ 劉國(guó)利，趙會(huì)濱，許貽燕.階躍SH波作用下半圓形凹陷地形的瞬態(tài)反應(yīng)：長(zhǎng)期解［J］.地震工程與工程振動(dòng)，1995，15（1）：92 －99.Liu Guo －li，Zhao Hui －bin，Xu Yi －yan.Transient response of semi －circular canyon under step SH wave：Long term solution［J］.Earthquake Engineering and Engineering Vibration，1995，15（1）：92 －99.

［10］ Herman H，Klosner J M.Transient response of a periodically supported cylindrical shell Immersed in a fluid medium［J］.Journal of Applied Mechanics，1965，32（3）：562 －568.

［11］ Geers T L.Excitation of an elastic cylindrical shell by a transient acoustic wave［J］.Journal of Applied Mechanics，1969，36（3）：459 －469.

［12］ 楊峻，宮全美，吳世明，等.飽和土體中圓柱形孔洞的動(dòng)力分析［J］.上海力學(xué)，1996，17（1）：37 －45.Yang Jun，Gong Qun －mei，Wu Shi －ming，et al.Dynamic analysis of cylindrical holes in saturated soil［J］.Shanghai Mechanics，1996，17（1）：37 －45.

［13］ Forrestal M J，Sagartz M J.Radiated pressure in an acoustic medium produced by pulsed cylindrical and spherical shells［J］.Journal of Applied Mechanics，1971，38（4）：1057 －1060.

［14］ Moodie T B，Barclay D W.Wave propagation from a cylindrical cavity［J］.Acta Mechanica，1977，27（1/2/3/4）：103 －120.

［15］ Moodie T B，Haddow J B，Mioduchowski A，et al.Plane elastic waves henerated by dynamical loading applied to edge of circular hole［J］.Journal of Applied Mechanics，1981，48（3）：577 －581.

［16］ Akkas N，Erdogogan F.The residual variable method applied to the diffusion equation in cylindrical coordinates［J］.Acta Mechanic，1989，79（3/4）：207 －219.

［17］ Zakout U，Akkas N.Transient response of a cylindrical cavity with and without a bonded shell in an infinite elastic medium［J］.International Journal of Engineering Science，1997，35（12/13）：1203 －1220.

［18］ Feldgun V R，Kochetkov A V，Karinski Y S，et al.Internal blast loading in a buried lined tunnel［J］.International Journal of Impact Engineering，2008，35（3）：172 －183.

［19］ 高盟，高廣運(yùn)，王瀅，等.內(nèi)部荷載作用下圓柱形孔洞的動(dòng)力響應(yīng)解答［J］.力學(xué)季刊，2009，30（2）：266 －272.Gao Meng，Gao Guang－yun，Wang Ying，et al.Solution on dynamic response of a cylindrical cavity under internal load［J］.Chinese Quarterly of Mechanics，2009，30（2）：266 －272.

［20］ 高盟，高廣運(yùn)，王瀅，等.飽和土與襯砌動(dòng)力相互作用的圓柱形孔洞內(nèi)源問(wèn)題解答［J］.固體力學(xué)學(xué)報(bào)，2009，30（5）：481 －488.Gao Meng，Gao Guang－yun，Wang Ying，et al.Dynamic solutions of cylindrical cavities with the lining under internal load in the saturated soil［J］.Chinese Journal of Solid Mechanics，2009，30（5）：481 －487.

［21］ 蔡袁強(qiáng)，陳成振，孫宏磊.黏彈性飽和土中隧道在爆炸荷載作用下的動(dòng)力響應(yīng)［J］.浙江大學(xué)學(xué)報(bào)：工學(xué)版，2011，45（9）：1657 －1663.Cai Yuan －qiang，Chen Cheng －zhen，Sun Hong －lei.Dynamic response of tunnel in viscoelastic saturated soil subjec －ted to blast loads［J］.Journal of Zhejiang University：Engineering Science，2011，45（9）：1657 －1663.

［22］ Haddow J B，Lorimer S A，Tait R J.Nonlinear axial shear wave propagation in a hyperelastic incompressible solid［J］.Acta Mechanica，1987，66（1/2/3/4）：205 －216.

［23］ Haddow J B，Jiang L.Finite amplitude azimuthal shear waves in a compressible hyperelastic solid［J］.Journal of Applied Mechanics，2001，68（2）：145 －152.

［24］ Watanabe K，Payton R G.SH wave in a cylindrically anisotropic elastic solid a general solution for a point source［J］.Wave Motion，1996，25（2）：197 －212.

［25］ Barclay D W .Wavefront expansion for non －linear axial shear wave propagation［J］.International Journal of Nonlinear Mechanics，1998，33（2）：259 －274.

［26］ Barclay D W.Shock front analysis for axial shear wave propagation in a hyperelastic incompressible solid［J］.Acta Mechanica，1999，133（1/2/3/4）：105 －129.

［27］ Barclay D W.Shock calculations for axially symmetric shear wave propagation in a hyperelastic incompressible solid［J］.International Journal of Non －linear Mechanics，2004，39（1）：101 －121.

［28］ Durbin F.Numerical inversion of Laplace transformation：An efficient improvement to Durbin and Abatep’s method［J］.The Computer Journal，1974，17（4）：371 －376.

［29］ 梁昆淼.數(shù)學(xué)物理方法［M］.北京：高等教育出版社，2010.

［30］ Watson G N.Theory of Bessel function［M］.Cambridge，UK：Cambridge University Press，1995.