圓形洞室在徑向非均勻荷載下的瞬態響應

耿大新， 陶 彪， 胡文韜

(華東交通大學 土木與建筑學院，南昌 330013)

隨著地下空間深入開發，地下管線、地鐵隧道、水下隧道等結構日益增多，彈性介質中含空腔或殼狀結構的動力響應問題一直是研究的熱點。然而這些結構多受一個隨機的內部動態荷載作用，并且結構中大多數內源荷載并不是環向均勻的。有的結構由于受外部約束，荷載傳遞不均，向某方向集中匯聚。因此確定非均勻瞬態荷載引起的結構動力響應是地下工程領域中十分重要的問題。

針對圓形洞室與薄壁殼體等在軸對稱瞬態荷載下的響應問題，目前主要的分析方法有解析法與數值法。解析法方面已有許多學者利用波函數展開法、積分變換法等方法進行研究。早在20世紀90年代初Senjuntichai等[1]采用波函數展開法，推導了全空間圓柱形腔體在三種不同類型軸對稱荷載下的徑向位移、應力、孔隙壓力的精確通解，并通過數值反演拉普拉斯解，而得到時域解。在其基礎上，Gao等[2-3]進一步研究了在襯砌洞室內部作用三種荷載情況下的動力響應。隨后，Engin等[4]和陸建飛等[5]將洞室穩態下的解答，分別推廣到半空間及任意洞室中，研究了彈性空間內土骨架的位移應力表達。此外，也有類似研究考慮了瞬態彈性波入射的情況，Karinski[6]、王瀅[7]和李偉華[8]分別研究了襯砌洞室在不同種瞬態波散射下應力、位移、動應力集中的時域解，并考慮了剛度、襯砌厚度等因素對動應力集中的影響。翟朝嬌等[9]針對反平面沖擊荷載作用下洞室的瞬態響應問題進行探討，研究分析了沿z軸方向瞬態荷載對土體動力響應的變化規律。然而當動荷載強度非常大，應變水平較高，對于這類問題一般常采用數值方法進行研究。Feldgun等[10]開創性的利用戈杜諾夫變分差分方法，分析了彈性、塑性及多孔介質中襯砌洞室的動力響應問題，并就解答的正確性與Glen等[11]的結果相印證。目前來說，已有研究大都針對洞室在環向均勻沖擊荷載下的動力響應解析解求解[12-13]或者利用數值法結合有限元、邊界元等對洞室內部均勻爆炸荷載作用下的動力響應分析[14]。但對徑向非均勻均布荷載下的動力響應與波動特性，卻一直未見類似研究，為此建立一種針對徑向非均勻荷載作用下圓形洞室瞬態響應的計算方法，對隧道工程與地下空間領域具有深遠意義。

本文將基于彈性介質波動理論，運用波函數展開法與Laplace變換法，根據圓形洞室內表面非均勻應力邊界條件，求解出單位脈沖荷載下圓形洞室的動力響應解答。給出了全空間洞室中應力和位移場在時域內的數值解，并通過算例，分析了徑向非均勻瞬態荷載下的波動特性以及剪切模量、不同角度對應力位移場的影響。

1 圓形洞室模型及波場求解

假定巖體為單相彈性介質，無限長圓柱形洞室埋置其中。因而，洞室內表面作用非對稱瞬態荷載的動力響應問題可簡化為平面應變問題。徑向非均勻荷載隨時間t,環向角度θ變化，如圖1所示，r為極軸，a為圓形洞室內半徑。

圖1 徑向非均勻荷載洞室模型Fig.1 Local concentrated load cavern model

圖2 瞬態三角形脈沖荷載Fig.2 Transient triangular pulse load

2 巖體的控制方程

非均勻瞬態荷載F(t)從洞室內部傳遞至洞室邊界后，在巖體中產生向外傳播的膨脹波，洞室外的巖體視為單相彈性介質，其幾何方程為[15]：

(1)

假定符合理想線彈性模型關系為：

σij=λ1δijεkk+2μ1εij

(2)

σr=(λ1+2μ1)εr+λ1εθ

(3a)

σθ=(λ1+2μ1)εθ+λ1εr

(3b)

σrθ=2μ1εrθ

(3c)

在極坐標系下，其振動方程可表示為：

(4)

(5)

式中：εr、εθ、εrθ分別為巖體徑向應變、環向應變、切向應變；σr、σθ、σrθ為巖體徑向應力、環向應力、切向應力；λ1為巖體的Lame常數；δij為Kronecker參數，當i≠j時δij=0，i=j時δji=1；ur、uθ為巖體介質的徑向位移、環向位移。

將式(1)～(3)代入方程式(4)～(5)，得到以位移表示巖體的控制方程為：

(6)

(7)

由于ur(r,θ,t)和uθ(r,θ,t)是相互耦合的，為了解耦引入巖體部分的位移標量勢函數φ(r,θ,t)和矢量勢函數ψ(r,θ,t)，根據Helmholtz矢量分解定理有：

(8)

位移表示為：

(9)

對時間t進行Laplace變換和逆變換為：

(10)

(11a)

(11b)

由式(11)整理可得下式：

(12)

由式(12)可得到勢函數Laplace變換后，滿足如下的Helmholtz方程：

(13a)

(13b)

(14a)

(14b)

在極坐標下，Laplace算子與Laplace變換后的勢函數可表示為：

(15a)

(15b)

采用分離變量法，對于線性系統中變換后勢函數的解可表達為如下形式[16]：

(15c)

將式(15)代入式(13)，整理后可分解為兩個不同的波數方程：

(16a)

(16b)

式(16)是n階虛宗量Bessel函數，ki為波數，下標i=1,2。

式(16)中，勢函數在極坐標下的通解可用Bessel函數線性組合的形式表達：

In(k1r)(A2cosnθ+B2sinnθ)]

(17a)

In(k2r)(C2cosnθ+D2sinnθ)]

(17b)

式中：Ai、Bi、Ci、Di，i=1,2為待定系數，In(·)為第一類虛宗量Bessel函數，Kn(·)為第二類虛宗量Bessel函數。

根據本文假設，在無限空間中，在r→∞時，須滿足u,v→0，In不滿足假設，因此待定系數A2、B2、C2、D2=0，則巖體勢函數可表示為：

(18)

將式(18)代入表達式(9)，并考慮本構關系式(2)，可得極坐標下巖體中位移、應力以勢函數的表達如下：

(19)

3 邊界條件以及數值求解

本文研究在無限空間中半徑為a的圓柱形洞室在內部受非均勻性沖擊荷載作用如圖1所示。確定待定系數與波場關系后，利用邊界條件求解上述勢函數中的待定系數，考慮圓形洞室與巖體交界面的邊界條件可得：

(20)

本文為求徑向非均勻瞬態荷載下的動力響應表達，荷載形式如下所示，將脈沖荷載表達式進行Laplace變換，得到Laplace變換域下表達，其表達式為：

(21)

式中：f(θ)可為任意形式的徑向非均勻動荷載；b,c為徑向非均勻荷載形狀參數；T為三角形脈沖荷載的周期。

(22)

式中：P11、P12、P13、P14、E11、E12、E13、E14為系數項，具體表達詳見附錄。通過矩陣求解出待定系數后，代入式(19)中即可求出頻域下應力與位移的解。

求得頻域解答后，由于半解析解的形式直接進行Laplace逆變換較為困難，利用Laplace數值逆變換，轉換為時域中的解，本文采取的是Durbin[17]數值逆變換方法，其變換表示為：

(23)

根據收斂準則確定NSUM的范圍：

(24)

4 計算結果與算例分析

4.1 結果驗算

(1)為了驗證本文計算的合理性與正確性，將本文計算結果與文獻[1]結果進行對比。為此將本文求解的徑向非均勻瞬態荷載退化為徑向均勻的瞬態荷載，即選取b=c=1，F0=0.1 MPa，λ=1.73×108Pa，a=3。計算所得無量綱環向應力、徑向位移隨時間分布圖，如圖3所示。本文在t*=8與t*=18峰值位置結果稍小，原因是參考文獻[1]采用的數值逆變換方法與本文不同，在數值逆變換時無法完全擬合出對應的實部項α值。其次激勵函數局部時間換算至全局時間過程取值偏大，將造成波數到達峰值時間略微推后。由圖3可知，本文在環向應力的變化趨勢上與文獻[1]基本一致，說明了本文公式推導結果的合理性。

圖3 本文退化為均勻的瞬態荷載與文獻[1]結果比較Fig.3 Comparison of variation of hoop stress with time between present work and Ref.[1]

(2)另將本文計算結果與文獻[18]數值模擬結果相對比，進一步驗證推導結果的正確性。本文模型取值令F0=0.068 7 MPa，a=0.02 m，b=5,c=4，土體物理力學參數等與文獻[18]一致，計算所得沿洞室徑向不同位置處的正應力，如圖4所示。圖中r表示θ=0°位置上距洞室內表面的距離，由圖可知正應力分布曲線與文獻[18]數值模擬結果吻合較好，由此說明了本文推導結果的正確性。

圖4 本文計算結果與文獻[18]結果比較Fig.4 Result comparisons between calculation results presented in this paper and those in Ref.[18]

4.2 算例分析

考慮洞室r=a處(內表面)各個角度θ=0°,30°,60°,90°對土體應力位移響應的影響，土體的基本參數，見表1。

表1 計算參數Tab.1 Calculation parameters

圖5表示洞室內壁，不同角度上的應力與位移響應值。從圖5(a)曲線可知，t*>10時各角度位移響應曲線基本趨近于0。隨著θ角0°→90°變化，環向位移先增大后減小，0°、90°的響應值為零。由于環向相互擠壓變形，導致uθ隨時間增大，各角度到達峰值的時間不同并且數值上差距較大。特別地，當接近t*=2時，30°處位移開始減小，60°位置處的位移剛到達最大值。表明在30°位置位移開始減小時，使得同一內徑上60°位置的環向位移值有一小幅上升的階段，每個θ角的振動是獨立且異步的。圖5(b)中可以看出不同角度對徑向位移的影響十分顯著。當θ從0°→30°時，徑向位移峰值減小近20%。而θ持續增加到90°時，位移峰值衰減速率逐漸放緩，愈接近90°位置，衰減越慢，θ對徑向位移的影響越小。與環向位移不同是各角度到達峰值的時間相同，振動是同步的。其中θ=0°時位移的幅值最大，與假設激勵函數性質相同。在考慮徑向非均勻的脈沖荷載最大位移響應時，應充分考慮0°即荷載集中位置的情況。

圖5(c)、(d)應力在t*=1附近到達最大值，且各角度到達峰值的時間都相同，環向與徑向均處于受壓狀態。隨時間推移，環向應力由相互擠壓狀態變為環向拉伸狀態，直至趨于穩定。θ角0°→90°時應力幅值逐漸減小，變化規律基本相同。并且θ改變時不僅影響應力峰值，也導致不同角度上應力衰減的速率有所差異。由圖(c)、(d)可知θ從0～90°方向上其衰減速率逐漸減小，在θ=0°方向上能量擴散速度最快。圖5(e)所示剪切應力隨時間波動遞減，呈往復態勢。應力減小至t*=2時有一顯著增大過程，越靠近90°位置往復性越明顯，剪切響應越小，直到趨近0。此現象類似于環向應力，不同之處是環向應力隨時間推移沿洞室環向拉壓狀態改變，而剪應力雖呈往復態勢，但剪切方向不會改變。

圖5 不同角度對洞室內表面(r=a)位移應力響應的影響Fig.5 Displacement stress response of the chamber surface at different angles

考慮不同徑向距離對土體應力位移響應的影響，圖6給出r*=r/a,r分別取1倍、1.2倍、1.4倍、和1.6倍洞徑，角度取θ=0°或30°的應力位移隨時間變化分布曲線。由于θ=0°時環向位移為零，無法考慮環向位移的波動關系，故任取一角度研究其波動關系，此算例取30°時的環向位移進行分析，下同。值得注意的是，從圖(a)、(b)位移曲線來看，盡管r*增加位移響應進入峰值的時間各不相同，但是不同r*對應的曲線幾乎是同一時間衰減至0。距離荷載中心位置越遠，振動到達峰值的時間越長，振動衰減的越快。

從徑向與環向應力圖來看，圖6(c)、(d)曲線都是隨r*增大響應逐漸減小。徑向應力波動周期比環向周期明顯要短，在t*=2左右就衰減至零。環向應力波動時間之所以更長，是因為在同一土環內，環向會受到相互擠壓作用，造成波動持續的時間更長。然而無限空間下，徑向不存在相互作用，應力隨著入射波不斷向遠處傳播而越來越小，波動周期的也更短。

圖6 不同徑向距離對同一角度下位移應力響應的影響Fig.6 Displacement stress response at different radial distances at the same angle

對不同拉梅常數下的應力位移時域響應進行分析，剪切模量μ0=1.15×108Pa，μ取0.1、0.3、0.5倍的剪切模量，θ分別取0°或30°，其他參數如表1所示。圖7(a)、(b)計算結果表明不同剪切模量對響應有顯著影響。對于三角形脈沖荷載，剪切模量增加不僅使位移響應值提前到達最大值，并且也導致峰值大幅度降低。這表明介質剛度越大，對變形吸收越多，位移響應也逐漸減弱。而剪切模量變化對位移響應的衰減速率沒有影響。當剪切模量較小時，位移響應的周期也越長，剪切模量越大，波動的時間越短，響應越早趨于平穩。

圖7(c)為剪切模量對應力響應的影響，由圖可知剪切模量的增加使環向應力值提前到達最大值，當剪切模量增大至0.3μ0后，剪切模量對響應到達峰值的時間的影響降低，0.3μ0以后的曲線幾乎同一時刻到達峰值。而應力峰值的變化規律與位移響應規律相反，隨剪切模量增大應力峰值越來越大，響應也愈加明顯。從圖中曲線可知，剪切模量對其衰減速率有顯著影響，剪切模量增大時，應力衰減速率逐漸增大，相對于位移曲線沒有表現出這種規律。并且位移與應力曲線都是在t*=1左右達到峰值，符合時域特征激勵函數的假設。

圖7 不同剪切模量對同一角度下位移應力響應的影響Fig.7 Displacement stress response of different shear modulus at the same angle

考慮不同時刻下的沿洞室環向應力位移時域響應，取洞室內表面r=a處。土體參數見表1，t*分別取0.5 s、1 s、3 s、5 s、10 s五個瞬時點。從圖8(a)可知，在波動初始階段t*=0.5時，環向在接近15°位置位移最大，在t*=1時環向位移最大值在30°。當t*=3時，環向位移最大值在60°位置出現，這說明隨時間推移環向位移極值會隨角度發生改變，各時刻下每個角度上的振動都是獨立的，并不是同步到達峰值后衰減。在施加荷載初始階段t*<3時，uθ出現峰值的角度隨時間增大而增大。徑向位移曲線分布如圖(b)所示，不同t*下，最大位移出現在洞室θ=0°和θ=180°位置。時間變化并不影響徑向位移峰值的位置，并且越接近90°位置位移響應衰減越多。從位移響應可知，在t*>3后位移波動趨于平穩，徑向位移的響應值遠大于環向位移。

圖8(b)、(c)為環向應力與徑向正應力沿洞室內表面的分布。在三角形脈沖荷載作用下，環向應力狀態會發生改變，初始階段環向受壓，隨環向響應增大，環向應力逐漸變為受拉。并且環向應力的極值隨時間轉動。當t*=0.5，環向應力在0°和180°取極大值，當t*=3時，極大值發生在90°和270°位置。與環向應力規律不同的是徑向正應力絕大數都是受壓狀態，時間不會改變徑向拉壓狀態。由此可以看出在洞室內表面位置，不同瞬時的差別，直接影響曲線出現極值的位置。

圖8 不同瞬時對洞室內表面(r=a)位移應力響應的影響Fig.8 Displacement stress response of the chamber surface at different time

5 結 論

本文基于彈性動力學理論，建立無限介質中圓形隧道的徑向非均勻瞬態荷載模型。采用波函數展開法，并利用三角函數正交性，求得洞室在非均勻瞬態荷載下頻域內的半解析解，通過Laplace數值逆變換得到時域下應力位移的解答。此外，通過算例還分析了不同角度、介質模量等對響應的影響。并得到了以下結論：

(1)徑向非均勻荷載作用下，環向位移響應隨時間推移，各個角度的振動都是異步的，在各角度的振動到達峰值的時間各不相同，其中角度為0、π/2、π等位置處的環向位移為零。環向應力在拉壓狀態改變后峰值衰減也具有異步性，在徑向應力和位移中并未出現類似現象。且在徑向非均勻荷載作用時，由于環向響應的異步性，也使得結構發生破壞的幾率增大。

(2)在徑向非均勻荷載作用時，洞室環向應力與位移的極值位置隨時間推移發生旋轉。在t*<1荷載施加初始階段，環向受到相互擠壓，應力極值所在位置在0-π內隨時間增大而增大。t*>1荷載釋放過程中，環向由受壓狀態變為受拉狀態，位移響應逐漸增強，位移的極值位置也開始旋轉，位移旋轉速率明顯慢于應力。

(3)剪切模量對圓形洞室響應的影響較大，不僅使應力位移響應的峰值提前到達，并且極大的影響響應的幅值。剪切模量越大時，位移幅值越低，應力的幅值越高。

(4)環向位移與徑向相比，其響應恢復到穩定的時間更長，環向振動的周期大于徑向。由于應力拉壓狀態發生改變，當t*>4時，洞室內表面出現微小的負位移。θ在0與π荷載徑向集中處，不論是徑向或環向應力其波動幅值明顯大于其他方向，且徑向應力、位移的最值大于環向。

附錄：