圓形洞室在徑向非均勻荷載下的瞬態(tài)響應(yīng)

耿大新， 陶 彪， 胡文韜

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

隨著地下空間深入開發(fā)，地下管線、地鐵隧道、水下隧道等結(jié)構(gòu)日益增多，彈性介質(zhì)中含空腔或殼狀結(jié)構(gòu)的動力響應(yīng)問題一直是研究的熱點(diǎn)。然而這些結(jié)構(gòu)多受一個(gè)隨機(jī)的內(nèi)部動態(tài)荷載作用，并且結(jié)構(gòu)中大多數(shù)內(nèi)源荷載并不是環(huán)向均勻的。有的結(jié)構(gòu)由于受外部約束，荷載傳遞不均，向某方向集中匯聚。因此確定非均勻瞬態(tài)荷載引起的結(jié)構(gòu)動力響應(yīng)是地下工程領(lǐng)域中十分重要的問題。

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

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

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

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

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

圖2 瞬態(tài)三角形脈沖荷載Fig.2 Transient triangular pulse load

2 巖體的控制方程

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

(1)

假定符合理想線彈性模型關(guān)系為：

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

在極坐標(biāo)系下，其振動方程可表示為：

(4)

(5)

式中：εr、εθ、εrθ分別為巖體徑向應(yīng)變、環(huán)向應(yīng)變、切向應(yīng)變；σr、σθ、σrθ為巖體徑向應(yīng)力、環(huán)向應(yīng)力、切向應(yīng)力；λ1為巖體的Lame常數(shù)；δij為Kronecker參數(shù)，當(dāng)i≠j時(shí)δij=0，i=j時(shí)δji=1；ur、uθ為巖體介質(zhì)的徑向位移、環(huán)向位移。

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

(6)

(7)

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

(8)

位移表示為：

(9)

對時(shí)間t進(jìn)行Laplace變換和逆變換為：

(10)

(11a)

(11b)

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

(12)

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

(13a)

(13b)

(14a)

(14b)

在極坐標(biāo)下，Laplace算子與Laplace變換后的勢函數(shù)可表示為：

(15a)

(15b)

采用分離變量法，對于線性系統(tǒng)中變換后勢函數(shù)的解可表達(dá)為如下形式[16]：

(15c)

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

(16a)

(16b)

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

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

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

(17a)

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

(17b)

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

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

(18)

將式(18)代入表達(dá)式(9)，并考慮本構(gòu)關(guān)系式(2)，可得極坐標(biāo)下巖體中位移、應(yīng)力以勢函數(shù)的表達(dá)如下：

(19)

3 邊界條件以及數(shù)值求解

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

(20)

本文為求徑向非均勻瞬態(tài)荷載下的動力響應(yīng)表達(dá)，荷載形式如下所示，將脈沖荷載表達(dá)式進(jìn)行Laplace變換，得到Laplace變換域下表達(dá)，其表達(dá)式為：

(21)

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

(22)

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

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

(23)

根據(jù)收斂準(zhǔn)則確定NSUM的范圍：

(24)

4 計(jì)算結(jié)果與算例分析

4.1 結(jié)果驗(yàn)算

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

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

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

圖4 本文計(jì)算結(jié)果與文獻(xiàn)[18]結(jié)果比較Fig.4 Result comparisons between calculation results presented in this paper and those in Ref.[18]

4.2 算例分析

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

表1 計(jì)算參數(shù)Tab.1 Calculation parameters

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

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

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

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

從徑向與環(huán)向應(yīng)力圖來看，圖6(c)、(d)曲線都是隨r*增大響應(yīng)逐漸減小。徑向應(yīng)力波動周期比環(huán)向周期明顯要短，在t*=2左右就衰減至零。環(huán)向應(yīng)力波動時(shí)間之所以更長，是因?yàn)樵谕煌镰h(huán)內(nèi)，環(huán)向會受到相互擠壓作用，造成波動持續(xù)的時(shí)間更長。然而無限空間下，徑向不存在相互作用，應(yīng)力隨著入射波不斷向遠(yuǎn)處傳播而越來越小，波動周期的也更短。

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

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

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

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

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

圖8(b)、(c)為環(huán)向應(yīng)力與徑向正應(yīng)力沿洞室內(nèi)表面的分布。在三角形脈沖荷載作用下，環(huán)向應(yīng)力狀態(tài)會發(fā)生改變，初始階段環(huán)向受壓，隨環(huán)向響應(yīng)增大，環(huán)向應(yīng)力逐漸變?yōu)槭芾２⑶噎h(huán)向應(yīng)力的極值隨時(shí)間轉(zhuǎn)動。當(dāng)t*=0.5，環(huán)向應(yīng)力在0°和180°取極大值，當(dāng)t*=3時(shí)，極大值發(fā)生在90°和270°位置。與環(huán)向應(yīng)力規(guī)律不同的是徑向正應(yīng)力絕大數(shù)都是受壓狀態(tài)，時(shí)間不會改變徑向拉壓狀態(tài)。由此可以看出在洞室內(nèi)表面位置，不同瞬時(shí)的差別，直接影響曲線出現(xiàn)極值的位置。

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

5 結(jié) 論

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

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

(2)在徑向非均勻荷載作用時(shí)，洞室環(huán)向應(yīng)力與位移的極值位置隨時(shí)間推移發(fā)生旋轉(zhuǎn)。在t*<1荷載施加初始階段，環(huán)向受到相互擠壓，應(yīng)力極值所在位置在0-π內(nèi)隨時(shí)間增大而增大。t*>1荷載釋放過程中，環(huán)向由受壓狀態(tài)變?yōu)槭芾瓲顟B(tài)，位移響應(yīng)逐漸增強(qiáng)，位移的極值位置也開始旋轉(zhuǎn)，位移旋轉(zhuǎn)速率明顯慢于應(yīng)力。

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

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

附錄：