隧道地震預(yù)報(bào)波場(chǎng)的有限元數(shù)值模擬

王朝令，劉爭(zhēng)平，黃云艷，黃顯彬

1.四川農(nóng)業(yè)大學(xué)土木工程學(xué)院，成都 611830 2.西南交通大學(xué)土木工程學(xué)院，成都 610031

隧道地震預(yù)報(bào)波場(chǎng)的有限元數(shù)值模擬

王朝令1，2，劉爭(zhēng)平2，黃云艷1，黃顯彬1

1.四川農(nóng)業(yè)大學(xué)土木工程學(xué)院，成都 611830 2.西南交通大學(xué)土木工程學(xué)院，成都 610031

在勘探地球物理領(lǐng)域中，能快速、有效地?cái)?shù)值模擬2D和3D的地震全波場(chǎng)，包括同時(shí)模擬P波、S波、面波多分量波場(chǎng)特征的軟件并不是很多；但在結(jié)構(gòu)力學(xué)領(lǐng)域中，已有多種成熟的多波多分量有限元數(shù)值模擬大型商用軟件用于求解彈性本構(gòu)方程，如ANSYS，PLAXIS等。從方程結(jié)構(gòu)來(lái)看，彈性本構(gòu)方程與完全彈性波動(dòng)方程相比，只多一項(xiàng)對(duì)時(shí)間的一階微分項(xiàng)——阻尼項(xiàng)。因此，調(diào)整合理的介質(zhì)參數(shù)——Rayleigh系數(shù)，令阻尼項(xiàng)近似為0，彈性本構(gòu)方程就蛻化為完全彈性波動(dòng)方程。隧道地震波場(chǎng)是在巖體中傳播，具有彈性參數(shù)較好、無(wú)常規(guī)地震勘探中覆蓋層等低速帶干擾等的優(yōu)點(diǎn)，采用ANSYS軟件對(duì)隧道地震波場(chǎng)進(jìn)行數(shù)值模擬，研究了頻散、阻尼與吸收邊界等數(shù)值模擬的相關(guān)技術(shù)。針對(duì)介質(zhì)吸收，比較了數(shù)值模擬中有無(wú)阻尼的時(shí)間記錄和波場(chǎng)快照，對(duì)隧道地震預(yù)報(bào)來(lái)說(shuō)，將阻尼系數(shù)設(shè)定為0是一種合理的假設(shè)；比較了實(shí)例中不同網(wǎng)格長(zhǎng)度與頻散的關(guān)系，當(dāng)網(wǎng)格長(zhǎng)度小于波長(zhǎng)的1/π倍時(shí)，才能消除頻散；通過(guò)由波方程的推導(dǎo)，引入了黏彈性邊界條件，通過(guò)實(shí)例計(jì)算證明它可以有效地吸收邊界反射；最后對(duì)隧道地震預(yù)報(bào)進(jìn)行了實(shí)例計(jì)算，算例表明采用ANSYS軟件可以有效地模擬隧道復(fù)雜地質(zhì)條件下全波場(chǎng)的激發(fā)傳播過(guò)程。

數(shù)值模擬；ANSYS；有限元；隧道地震預(yù)報(bào)；地震波場(chǎng)；阻尼；頻散；邊界條件

0 引言

隧道屬于地下建筑工程，其圍巖地質(zhì)情況復(fù)雜多變［1］，常遇到滑坡、崩坍、巖堆、巖洞、高應(yīng)力高強(qiáng)度地層、松散地層、軟土等不良地質(zhì)地段，以及膨脹地層、含水未固結(jié)圍巖、溶洞、斷層、巖爆、流沙和瓦斯溢出地層等特殊地質(zhì)地段，它們具有不同的性質(zhì)、形狀和規(guī)模，且成因和變異條件非常復(fù)雜，地質(zhì)條件具有突發(fā)性。設(shè)計(jì)和勘察所提供的地質(zhì)資料不可能自始至終符合實(shí)際情況，其詳細(xì)程度也不能完全達(dá)到隧道施工的具體要求，如果沒(méi)有可靠的超前地質(zhì)預(yù)報(bào)，潛在的地質(zhì)隱患就可能變成災(zāi)難。因此，隧道施工中安全和工程事故時(shí)有發(fā)生，輕則延長(zhǎng)工期、增加投資，重則毀壞機(jī)械設(shè)備，甚至造成人員傷亡，而且事故處理難度大。在隧道地震超前預(yù)報(bào)方面：20世紀(jì)90年代初，瑞士Amberg測(cè)量技術(shù)有限公司基于地震波反射原理開(kāi)發(fā)研制了一套超前預(yù)報(bào)系統(tǒng)設(shè)備——隧道地震預(yù)報(bào)（TSP）系統(tǒng)，可探測(cè)和預(yù)報(bào)圍巖軟硬變化，斷層、破碎帶、節(jié)理裂隙發(fā)育情況，含水、溶洞及圍巖類別等［2］，TSP在瑞士、德國(guó)、法國(guó)、意大利、奧地利、日本等發(fā)達(dá)國(guó)家的隧道施工中，得到了廣泛的應(yīng)用；Hanson等［3］采用真反射成像法（true reneetion tomo，TRT），該方法在觀測(cè)方式上采用的是空間多點(diǎn)接收和激發(fā)系統(tǒng)，檢波器和激發(fā)的炮點(diǎn)呈空間分布，以充分獲得空間波場(chǎng)信息，提高對(duì)前方不良地質(zhì)體的定位精度。

地震模擬技術(shù)可以分為物理模擬和數(shù)值模擬2種［4］，數(shù)值模擬技術(shù)主要包括波動(dòng)方程法和射線法兩大類。射線法模擬精度較低；波動(dòng)方程法實(shí)質(zhì)是求解地震波波動(dòng)方程，它能比較準(zhǔn)確地模擬任意復(fù)雜介質(zhì)中的地震波場(chǎng)，為研究地震波傳播機(jī)理和復(fù)雜底層的解釋提供更多的數(shù)學(xué)及物理依據(jù)。鑒于隧道地震波場(chǎng)的復(fù)雜性，開(kāi)展波場(chǎng)數(shù)值模擬技術(shù)是非常有意義的。

常用的地震波場(chǎng)數(shù)值模擬方法主要包括有限差分法（FDM）、有限元法（FEM）等，目前大多采用基于有限差分的彈性或黏彈性雙程波動(dòng)方程進(jìn)行計(jì)算。Boore［5］將有限差分法用于SH波非均勻介質(zhì)地震波的CT成像模擬；Carcione等［6］提出了線性黏彈性介質(zhì)中地震波傳播的模擬方法；Talezcr等［7］進(jìn)行了線性黏彈性介質(zhì)中地震波傳播的方法研究；Robertsson等［8］給出了黏彈性波有限差分模擬方法；Graves［9］給出了三維速度－應(yīng)力方程交錯(cuò)網(wǎng)格有限差分法彈性波傳播的模擬方法；Carcione等［10］研究了孔隙彈性介質(zhì)中Biot縱波的數(shù)值模擬問(wèn)題；Ozdenvar等［11］給出了孔隙彈性介質(zhì)中地震波的交錯(cuò)網(wǎng)格有限差分模擬方法；Carcione等［12］提出了孔隙黏彈性介質(zhì)中地震波傳播的交錯(cuò)網(wǎng)格有限差分模擬方法；何樵登等［13］研究了橫向各向同性介質(zhì)中地震波數(shù)值模擬方法；王延光等［14］給出了一種適于強(qiáng)橫向變速的高階差分正演模擬方法；裴正林等［15］給出了三維任意各向異性介質(zhì)中彈性波傳播的交錯(cuò)網(wǎng)格高階有限差分模擬方法；董良國(guó)［16］利用交錯(cuò)網(wǎng)格高階差分方法對(duì)起伏地表情況下的彈性波傳播進(jìn)行了數(shù)值模擬；戴志陽(yáng)等［17］利用褶積微分算子實(shí)現(xiàn)了對(duì)地震波場(chǎng)的數(shù)值模擬，并利用模型驗(yàn)證了其有效性和正確性；周曉華等［18］利用交錯(cuò)網(wǎng)格有限差分方法模擬了微動(dòng)信號(hào)，結(jié)果表明采用這種方法得到的實(shí)測(cè)結(jié)果更符合實(shí)測(cè)結(jié)果。

這些方法具有較強(qiáng)的地震波傳播模擬能力，也能很好地描述小尺度不均勻各向異性介質(zhì)（巖層的孔隙、裂隙等）中地震波的傳播特征。

在隧道波場(chǎng)數(shù)值模擬方面：Bohlen等［19－21］通過(guò)數(shù)值模擬在隧道中開(kāi)展集成地震成像系統(tǒng)（ISIS）進(jìn)行隧道超前預(yù)報(bào)，并通過(guò)3D有限差分的方法模擬了隧道超前預(yù)報(bào)，詳細(xì)地模擬了隧道邊墻上橫波與面波的轉(zhuǎn)換特征；劉晶波等［21］采用非線性有限元軟件LS-DYNA模擬了爆破對(duì)隧道結(jié)構(gòu)的影響，并與現(xiàn)場(chǎng)爆破進(jìn)行了對(duì)比，證明模擬方法的有效性；魯光銀等［22］采用高階交錯(cuò)網(wǎng)格差分方法進(jìn)行了隧道超前預(yù)報(bào)的全波場(chǎng)數(shù)值模擬計(jì)算。

目前，在勘探地球物理領(lǐng)域中，能快速、有效地?cái)?shù)值模擬2D和3D地震全波場(chǎng)，包括同時(shí)模擬P波、S波、面波多波多分量波場(chǎng)特征的軟件，并不是很多。但在結(jié)構(gòu)力學(xué)領(lǐng)域中，已有多種成熟的多波多分量有限元數(shù)值模擬軟件用于解彈性本構(gòu)方程，如ANSYS，PLAXIS等。從方程結(jié)構(gòu)來(lái)看，彈性本構(gòu)方程與完全彈性波動(dòng)方程相比，只多一項(xiàng)對(duì)時(shí)間的一階微分項(xiàng)——阻尼項(xiàng)。因此，調(diào)整合理的介質(zhì)參數(shù)——Rayleigh系數(shù)，令阻尼項(xiàng)近似為0，彈性本構(gòu)方程就蛻化為完全彈性波動(dòng)方程。且這些軟件中模擬震源可以用任意自定義數(shù)據(jù)加載，因此可方便地用現(xiàn)場(chǎng)實(shí)驗(yàn)和物理模擬提取的噪聲震源函數(shù)加載模擬，顯然比一般地球物理數(shù)值模擬中，主要用Ricker子波進(jìn)行模擬更接近于實(shí)際情況。數(shù)值模擬中的關(guān)鍵技術(shù)是介質(zhì)Rayleigh系數(shù)調(diào)整和吸收邊界條件的合理加載。

在本文中，筆者首先引入彈性力學(xué)的本構(gòu)、Cauchy、Navier3個(gè)基本方程，這3個(gè)方程涉及應(yīng)力、應(yīng)變和位移3個(gè)物理量，是彈性波傳播理論中的基礎(chǔ)方程；然后以這3個(gè)方程構(gòu)建波動(dòng)方程；最后基于由以上步驟的理論準(zhǔn)備，利用ANSYS有限元軟件實(shí)現(xiàn)對(duì)地震波場(chǎng)的數(shù)值模擬。

1 理論基礎(chǔ)

人工地震中，當(dāng)介質(zhì)（即物體）受到非零的外力（如人工激發(fā)震源）時(shí)，該外力要轉(zhuǎn)化為介質(zhì)內(nèi)部的應(yīng)力，并使彈性介質(zhì)內(nèi)部發(fā)生應(yīng)變和位移，形成地震波（即彈性波）場(chǎng)，介質(zhì)的應(yīng)力、應(yīng)變以及能量的改變都是動(dòng)態(tài)的運(yùn)動(dòng)過(guò)程。

1）本構(gòu)方程的矩陣表達(dá)形式［23］為

式中：σ6×1和ε6×1分別為單角標(biāo)形式的應(yīng)力和應(yīng)變矩陣；C6×6是物性矩陣，也稱為介質(zhì)的彈性性質(zhì)矩陣，它描述了介質(zhì)的彈性和相關(guān)物理性質(zhì)。

2）Cauchy方程的一維矩陣形式可以寫(xiě)成

式中：u3×1是位移矩陣；L6×3是偏導(dǎo)數(shù)算符矩陣。

3）Naiver方程：

其中：ρ是介質(zhì)密度；算符L2t表示對(duì)時(shí)間變量求二階偏導(dǎo) 數(shù)；ux，uy，uz表 示x，y，z方 向 的 位 移；［LxLyLz］T＝［?/?x?/?y?/?z］T；［FxFyFz］T表示力矢量。式（3）就是振動(dòng)介質(zhì)的“平動(dòng)運(yùn)動(dòng)方程”，即Navier方程。

將Cauchy方程（2）代入本構(gòu)方程式（1），再將得到的方程式代入Navier方程式（3）得

式中：位移矢量u的矩陣表達(dá)式為u3×1＝［uxuyuz］T；¨u3×1為節(jié)點(diǎn)的加速度矢量。式（4）即在外力作用下用位函數(shù)表示的波動(dòng)方程。在地震波場(chǎng)中的外力就是震源的激發(fā)力，求解波動(dòng)方程的問(wèn)題就是波的傳播問(wèn)題。

在有限元求解中，存在下列矩陣微分方程：

式中：m是質(zhì)量矩陣；R是阻尼矩陣；K是剛度矩陣；˙u為節(jié)點(diǎn)速度矢量。式（5）是將單元?jiǎng)偠认禂?shù)和單元力按照前面給出的方式相加得到的，而單元力包括外部力、初始應(yīng)力等，此式相比于波動(dòng)方程式（4）多了一個(gè)阻尼項(xiàng)R˙u，新的矩陣R和m也是按照一般方法由單元子矩陣組裝得到的，即式中：Re與me分別是單元阻尼矩陣和單元質(zhì)量矩陣；Ni為形狀函數(shù)矩陣；μ為材料的黏滯矩陣；Ω為積分區(qū)域。

一般將式（5）定義的質(zhì)量矩陣叫作一致性質(zhì)量矩陣，是由Zienkiewicz和Cheung提出的［24］，其中的矩陣Re和R也稱為一致性阻尼矩陣。在計(jì)算過(guò)程中，集中質(zhì)量矩陣更方便、更高效。實(shí)際上，阻尼矩陣R很難確定，通常將阻尼矩陣定義為剛度矩陣和質(zhì)量矩陣的線性組合，即

式中，質(zhì)量阻尼系數(shù)α和阻尼剛度系數(shù)β可通過(guò)實(shí)驗(yàn)確定。這樣的阻尼稱為Rayleigh阻尼。

與波動(dòng)方程式（4）形式類比，有限元問(wèn)題方程（5）除去一阻尼項(xiàng)外，二者具有完全相同的形式，表明二者求解的彈性波場(chǎng)問(wèn)題是等價(jià)的。

2 數(shù)值模擬實(shí)現(xiàn)策略

基于上述分析可知，瞬態(tài)動(dòng)力學(xué)分析的基本方程式如式（5）所示，加入時(shí)間變量后變換為如下形式：

相比于式（5），將其中的力矩陣F變?yōu)楹奢d矢量F（t）。

在ANSYS軟件中，采用Newmark時(shí)間積分方法（包括改進(jìn)的算法HHT）和中心差分時(shí)間積分法求解瞬態(tài)動(dòng)力學(xué)基本方程。

在采用ANSYS軟件實(shí)現(xiàn)彈性波場(chǎng)模擬中，主要牽涉以下幾個(gè)方面的問(wèn)題：1）單元選擇；2）介質(zhì)吸收；3）頻散效應(yīng)；4）邊界條件。下面分別進(jìn)行討論。

2.1 單元選擇

由于計(jì)算機(jī)能力的限制，筆者首先研究2D的數(shù)值模擬，ANSYS中針對(duì)結(jié)構(gòu)的二維單元包括Plane25、Plane42、Plane83、Plane182、Plane183。其中：Plane183是高階二維八節(jié)點(diǎn)單元，具有二次位移項(xiàng)，適于生成不規(guī)則網(wǎng)格模型，因此不選用此單元；Plane25、Plane83均是諧響應(yīng)對(duì)稱單元，在本文研究的模型中，由于震源加載的限制，不能滿足對(duì)此類單元加載載荷的要求，因此不選擇此2種單元。比較之下，可供選擇的平面單元只有2個(gè)：Plane42和Plane182。

Plane42和Plane182單元非常相近。Plane182單元具有選擇簡(jiǎn)化積分這一選項(xiàng)，這對(duì)于不可壓縮情況有助于防止網(wǎng)格被鎖死；此外Plane182單元具有4個(gè)用戶不可調(diào)節(jié)的自由度，相比于Plane42單元，在計(jì)算中具有較高的精確度。在本文模擬中，選擇Plane182單元進(jìn)行彈性波場(chǎng)的數(shù)值模擬。圖1是在分別基于Plane42和Plane182單元采用相同震源激發(fā)等條件時(shí)接收到的地震記錄，Plane182單元計(jì)算結(jié)果更好，但這種差異是比較微小的。

圖1 不同單元數(shù)值模擬Fig.1 Modeling in different elements

2.2 介質(zhì)吸收

ANSYS可以進(jìn)行加載阻尼的運(yùn)算。雖然筆者在本文中只研究彈性介質(zhì)，不需要考慮阻尼因素，但當(dāng)進(jìn)行地震波場(chǎng)模擬時(shí)，需要考慮Rayleigh阻尼所造成的介質(zhì)吸收問(wèn)題，否則會(huì)出現(xiàn)意想不到的結(jié)果。據(jù)式（8）知，要計(jì)算Rayleigh阻尼，必須定義α和β；但一般情況下，α和β不能直接獲得，通常由模態(tài)阻尼比ξ計(jì)算得出。模態(tài)阻尼比是特定振動(dòng)模式下實(shí)際阻尼與臨界阻尼的比值。若ωi是第i階的自然圓頻率，那么α和β滿足如下關(guān)系式：

關(guān)于模態(tài)阻尼比ξ，不少文獻(xiàn)中的研究表明取0.02～0.05是比較合適的；在一次模擬中，往往取一個(gè)常數(shù)；對(duì)于巖石等，一般取0.05。在許多實(shí)際問(wèn)題中，往往忽略質(zhì)量阻尼系數(shù)α。因此，將ωi＝2πfi代入式（11）可以得到

式中：f為頻率。

若加載主頻為100Hz，則由式（13）可以得到剛度阻尼系數(shù)β為1.59×10－4。大部分體系的剛度阻尼系數(shù)都是落在0～0.05這一范圍的，在此范圍內(nèi)進(jìn)行動(dòng)力學(xué)計(jì)算時(shí)，位移的變化不超過(guò)5%；β為零時(shí)，系統(tǒng)的位移最大。因此，出于保守的目的，當(dāng)系統(tǒng)承受沖擊荷載時(shí)，取β為0是很合理的逼近。

下面通過(guò)數(shù)值計(jì)算來(lái)研究在加載阻尼系數(shù)β為1.59×10－4和0（即不加載阻尼系數(shù)）時(shí)的模擬效果。模型如圖2所示，網(wǎng)格長(zhǎng)度為0.5m，震源主頻為300Hz，提取偏移距為10m、道間距1m的時(shí)間記錄。

圖2 阻尼加載模型圖Fig.2 Model figure of damping loading

圖3是加載阻尼和無(wú)阻尼情況下的波場(chǎng)快照?？梢钥闯觯?種情況下傳播時(shí)的差異比較小，且相同時(shí)刻波場(chǎng)也非常接近；圖3a中由于施加了阻尼，在波場(chǎng)擴(kuò)散中能量不斷減小，對(duì)比50ms時(shí)刻的波場(chǎng)快照尤其明顯。

圖4為提取的時(shí)間記錄，用來(lái)比較加載前后的影響。由圖4可以看出，加載阻尼和無(wú)阻尼的時(shí)間記錄非常接近，無(wú)阻尼時(shí)得到的剖面信噪比更高，反射波的同相軸更清晰。

由此可見(jiàn)，ANSYS可以很好地模擬阻尼介質(zhì)的波場(chǎng)。在地震超前預(yù)報(bào)的數(shù)值模擬中，將質(zhì)量阻尼系數(shù)α和阻尼剛度系數(shù)β都設(shè)定為0。

2.3 網(wǎng)格剖分的頻散效應(yīng)

對(duì)模型進(jìn)行網(wǎng)格劃分時(shí)，由于要進(jìn)行瞬態(tài)分析來(lái)模擬波場(chǎng)傳播，必須考慮由于波長(zhǎng)與網(wǎng)格長(zhǎng)度所產(chǎn)生的頻散效應(yīng)，這主要取決于網(wǎng)格長(zhǎng)度和時(shí)間步長(zhǎng)的大小。從波動(dòng)力學(xué)帶寬定理［25］得到

圖3 阻尼波場(chǎng)快照Fig.3 Snapshot of damping

式中：k為波數(shù)；Δt為時(shí)間步長(zhǎng)；Δx為網(wǎng)格長(zhǎng)度。式（15）表明，能量所在空間或時(shí)間越集中，其頻譜就越寬。

國(guó)內(nèi)外有不少學(xué)者從不同角度提出離散化準(zhǔn)則，如 Kausel［26］，Costantino等［27］和 Miller等［28］分別給出：

式中：λ為波長(zhǎng)；λmin為最小波長(zhǎng)。

圖4 阻尼加載時(shí)間記錄Fig.4 Time recording of damping loading

據(jù)宗福開(kāi)［25］的研究結(jié)果，由λmin＝v/fmax（v為速度，fmax為最高頻率），當(dāng)波在傳播方向劃分的網(wǎng)格數(shù)n大于10時(shí)，則得出有限元離散化準(zhǔn)則為

式（18）具有明確的物理意義，即要求網(wǎng)格長(zhǎng)度小于波長(zhǎng)的1/π倍。應(yīng)當(dāng)指出，根據(jù)波傳播邊界條件的性質(zhì)，通常Δx≤（1/π～1/8）λmin。

據(jù)以上結(jié)論，筆者進(jìn)行了網(wǎng)格長(zhǎng)度不同的有限元數(shù)值模擬，模型示意圖如圖5所示。為減小邊界對(duì)模擬結(jié)果的影響，模型4個(gè)邊界均加載黏彈性吸收邊界。模型背景設(shè)定橫波波速為2 000m/s，縱波波速為4 000m/s，密度為1 800kg/m3。由實(shí)際資料分析，震源主頻取300Hz。由于橫波的波長(zhǎng)較縱波的小，故模型的λmin取為橫波的主波長(zhǎng)，約為20 m。分別設(shè)定網(wǎng)格長(zhǎng)度為10.0、5.0、2.5、1.0和0.5 m，分析波場(chǎng)數(shù)值模擬的頻散效應(yīng)。

圖6是不同網(wǎng)格長(zhǎng)度、同一節(jié)點(diǎn)的位移數(shù)值。根據(jù)公式（15），當(dāng)網(wǎng)格長(zhǎng)度為10.0m＞20/π＝6.366m時(shí)，波場(chǎng)的位移數(shù)值解震蕩很嚴(yán)重。隨著網(wǎng)格長(zhǎng)度越來(lái)越短，波的震蕩現(xiàn)象逐漸減弱。當(dāng)網(wǎng)格長(zhǎng)度為5.0m時(shí)，相比于10.0m的計(jì)算值，波的震蕩變小了；當(dāng)網(wǎng)格長(zhǎng)度為2.5m時(shí)，仍然存在輕微的震蕩；當(dāng)網(wǎng)格長(zhǎng)度由1.0m進(jìn)一步變小到0.5 m時(shí)，模擬計(jì)算結(jié)果并沒(méi)有太大的改進(jìn)?？紤]到計(jì)算量的問(wèn)題，對(duì)于實(shí)際模型模擬時(shí)，采用最大網(wǎng)格長(zhǎng)度為1.0m的剖分網(wǎng)格就可以最大程度地減小頻散效應(yīng)，很好地模擬波場(chǎng)的傳播了。因此，在進(jìn)行模擬時(shí)，按照式（18）所要求的網(wǎng)格長(zhǎng)度進(jìn)行劃分，即可消除頻散的影響。

圖5 網(wǎng)格剖分模型圖Fig.5 Model figure of grid meshing

圖6 采樣點(diǎn)位移隨網(wǎng)格長(zhǎng)度變化圖Fig.6 Displacement of sampling point varied with meshing

2.4 邊界條件

采用有限元模擬時(shí)，通常采用時(shí)間域的吸收邊界，這種邊界由一維或二維條件導(dǎo)出，并使得沒(méi)有能量從邊界反射。其中：應(yīng)用最廣泛的是Lysmer等［29］引入的黏彈性邊界，它采用黏彈性阻尼來(lái)代替遠(yuǎn)場(chǎng)邊界；Novak等［30］提出了平面應(yīng)變邊界，這種邊界主要應(yīng)用到嵌入基礎(chǔ)和樁基中。但是，平面應(yīng)變邊界包含依賴于頻率的項(xiàng)，這在加載瞬態(tài)載荷的瞬態(tài)響應(yīng)分析中會(huì)變得比較復(fù)雜；相比之下黏彈性邊界應(yīng)用更廣泛，加載簡(jiǎn)單，且對(duì)整個(gè)計(jì)算所增加的計(jì)算量可以忽略不計(jì)，在采用黏彈性邊界時(shí)，極少會(huì)產(chǎn)生這種反射。

當(dāng)波穿過(guò)介質(zhì)時(shí)，波前的應(yīng)力產(chǎn)生在法向和切向2個(gè)分量上。假定在平面應(yīng)力條件下，徑向位移的波動(dòng)方程為［31］

其中：ur表示徑向位移；G為壓縮模量；r為半徑；Λ為L(zhǎng)ame常數(shù)。

式中：E為Young模量；v是Possion比。

將徑向位移ur用位移勢(shì)函數(shù)φ表示：

則式（19）可表示為

兩邊都對(duì)徑向r進(jìn)行積分，則波方程用φ表示為

可以得到估計(jì)表達(dá)式：

用f′表示頻率f的導(dǎo)數(shù)，則徑向位移ur、徑向應(yīng)力σr、徑向應(yīng)變?chǔ)舝、橫向應(yīng)變?chǔ)纽瓤梢杂妙l率f和半徑r表示：

對(duì)于任意的半徑rb，對(duì)頻率f的任意時(shí)刻導(dǎo)數(shù)都是其對(duì)幅角導(dǎo)數(shù)的相反值，則半徑rb處的位移值、速度值和加速度值可以表示為

在半徑rb處二階時(shí)間導(dǎo)數(shù)為

將式（33）代入式（29）可以得到f（rb，t）在邊界的應(yīng)力：

對(duì)式（34）做時(shí)間微分：

聯(lián)立式（34）和式（35），可以消去未知函數(shù)f：

此種在有限元中集合時(shí)間的邊界，在實(shí)現(xiàn)中，若采用彈簧形式將會(huì)非常方便。Wolf［32］提供了如圖7所示的彈簧。

圖7 彈簧模型Fig.7 Spring model

圖7中彈簧模型的波動(dòng)方程如下：

由式（37）可以得到

式（39）對(duì)時(shí)間求導(dǎo)：

將式（39）和式（40）代入到式（38）可以得到第一個(gè)自由度的波動(dòng)方程：

聯(lián)立式（36）—式（41）可以得到等價(jià)的剛度、阻尼和質(zhì)量：

可以將式（42）（43）（44）所得的值賦予彈簧，來(lái)實(shí)現(xiàn)黏彈性邊界。

歸納二維人工邊界等效物理系統(tǒng)的彈簧系數(shù)K和阻尼系數(shù)C在針對(duì)不同邊界時(shí)分別如下。

切向邊界：

法向邊界：

式中：KBN、KBT分別為彈簧的法向與切向剛度；CBT、CBN分別為彈簧的法向與切向阻尼；rb為波源至人工邊界點(diǎn)的距離；γT和γN分別為切向與法向黏彈性人工邊界參數(shù)。根據(jù)劉晶波等［33］的研究結(jié)果，人工邊界比較合適的參數(shù)γT取值范圍為［0.35，0.65］，γN的取值范圍為［0.8，1.2］。

在ANSYS軟件中，適合用來(lái)做黏彈性邊界的單元選用Combin14，它在1維、2維或3維應(yīng)用中具有軸向或扭轉(zhuǎn)的性能?？v向阻尼彈簧是單軸壓縮張力的單元，在每個(gè)節(jié)點(diǎn)有3個(gè)自由度：x、y和z軸。不考慮扭轉(zhuǎn)和彎曲，另外彈簧阻尼單元沒(méi)有質(zhì)量矩陣，圖8為Combin14單元的示意圖和此單元在ANSYS數(shù)值模擬計(jì)算中的加載方式。

圖8 Combin14單元示意圖（a）和邊界加載方式（b）Fig.8 Combin14unit schematic（a）and boundary loading mode（b）

圖9 邊界比較模型圖Fig.9 Model figure of boundary comparing

按照式（19）和式（21）將Combin14在模型周邊布設(shè)，利用如圖9所示的模型計(jì)算四周是Dirichlet邊界和黏彈性邊界的數(shù)值模擬。震源為100Hz的Ricker子波，加載方向豎直向下，網(wǎng)格長(zhǎng)度0.5m，觀測(cè)系統(tǒng)以震源為中心，兩邊等間距，道間距為2 m，道數(shù)為48道，分別取得震源兩邊的時(shí)間記錄，如圖10所示。

在模型四周施加了黏彈性邊界時(shí)，直達(dá)P波和直達(dá)S波都受到了抑制，且相較之下，S波受影響更大。這是由于黏彈性邊界彈簧的存在使除直達(dá)P、S波之外的邊界反射被吸收的緣故，證明了黏彈性邊界可以有效地吸收反射。在默認(rèn)位移為0的Dirichlet邊界所得到的記錄中（圖11），各邊界的反射波返回到接收排列，使得對(duì)波的辨識(shí)和同相軸的判別都變得比較困難。由圖10時(shí)間記錄和圖11波場(chǎng)快照的相同時(shí)刻的對(duì)比可以看出，Dirichlet邊界的波場(chǎng)快照中各種反射波夾雜在一起，使得波場(chǎng)變得非常復(fù)雜，不易辨識(shí)，增加了波場(chǎng)分析的難度。

圖10 2種邊界的時(shí)間記錄Fig.10 Time record of two boundaries

圖11 2種邊界的波場(chǎng)快照Fig.11 Wavefield snapshot of two boundaries

3 計(jì)算實(shí)例

圖12為邊墻激發(fā)二維隧道數(shù)值模型。模型隧道開(kāi)挖工作面前方具有一模擬斷層垂直傾角的低波阻抗反射帶。隧道模型區(qū)域I中的縱、橫波速度分別設(shè)為4 000m/s和2 000m/s，低波阻抗區(qū)域II的縱、橫波速度為2 000m/s和1 000m/s，另一側(cè)區(qū)域III的縱、橫波速度設(shè)為3 000m/s和1 500m/s；其觀測(cè)系統(tǒng)中激發(fā)震源位于邊墻上，接收排列的48道檢波器也設(shè)置在邊墻上，并與震源在同一直線上，最小偏移距10m，道間距1m。由互換原理可知，若將接收與激發(fā)換位，即為T(mén)SP觀測(cè)系統(tǒng)。模擬震源為主頻300Hz的Ricker。采樣間隔為0.1ms，采樣視窗長(zhǎng)0.1s，網(wǎng)格長(zhǎng)度0.5m。模型邊界均采用黏彈性邊界條件。

圖12 邊墻激發(fā)數(shù)值模型Fig.12 Numerical model of sidewall excited

圖13 實(shí)例波場(chǎng)快照Fig.13 Snapshot of the example

圖13是邊墻激發(fā)所產(chǎn)生的波場(chǎng)快照。由圖13可見(jiàn)，震源激發(fā)的橫波能量遠(yuǎn)大于縱波能量，所以波場(chǎng)快照中縱波傳播和反射特征表現(xiàn)不如橫波那么明顯。由圖13分析，P波和S波的隧道波場(chǎng)特征有很大的差異，因此對(duì)二者的特征分述如下。首先，當(dāng)t＝5ms時(shí)，傳播中的P、S波還沒(méi)有完全分離；t＝10ms、t＝15ms時(shí)，由于P波比S波速度快，二者已經(jīng)完全分離。對(duì)于P波：當(dāng)t＝23ms時(shí)，P波傳播到掌子面；t＝30ms時(shí)，P波經(jīng)過(guò)掌子面，并繼續(xù)向前方介質(zhì)傳播；t＝37ms時(shí)，P波到達(dá)第一波阻抗變化界面，部分能量反射回來(lái)，另一部分則透射過(guò)界面繼續(xù)向前傳播，反射回來(lái)的P波將會(huì)被沿邊墻布置的檢波器觀測(cè)排列接收。對(duì)于S波，激發(fā)的S波一部分能量以體波的形式向邊墻外的介質(zhì)體內(nèi)傳播，另一部分能量則轉(zhuǎn)換為瑞利（Rayleigh）面波，沿邊墻向前傳播。t＝40ms時(shí)，S波到達(dá)掌子面。t＝46ms、t＝50ms時(shí)，瑞利面波的部分能量在掌子面轉(zhuǎn)換為S波，新產(chǎn)生的S波以體波的波陣面繼續(xù)向前傳播［20］，部分能量則繼續(xù)以面波的波陣面沿邊墻向后傳播，可被沿邊墻布置的檢波器觀測(cè)排列接收。這類來(lái)自于掌子面的反射面波被稱為SRR波，顯然對(duì)于提取掌子面前方介質(zhì)的反射S波信息來(lái)看，是干擾波。t＝55ms、t＝58ms時(shí)，轉(zhuǎn)換S波繼續(xù)向前傳播，直到t＝75ms時(shí)刻，轉(zhuǎn)換S波到達(dá)第一反射界面，部分能量發(fā)生反射回傳，部分發(fā)生透射繼續(xù)向前傳播。反射回傳的S波到達(dá)隧道邊墻時(shí)，部分能量將再次轉(zhuǎn)換為面波，并會(huì)被沿邊墻布置的檢波器觀測(cè)排列接收。分析表明，隧道邊墻激發(fā)條件下，觀測(cè)排列記錄的反射S波實(shí)際上是由S波轉(zhuǎn)換為面波，面波再轉(zhuǎn)換為S波，最后再轉(zhuǎn)換為面波的復(fù)雜轉(zhuǎn)換波，被稱為SRSR波。波場(chǎng)快照還顯示，SRSR波的能量遠(yuǎn)大于反射P波。

圖14是邊墻激發(fā)數(shù)值模型位移x分量、y分量的時(shí)間記錄。圖14中的同相軸主要包括反射縱波、反射橫波和面波轉(zhuǎn)換橫波，與圖13的波場(chǎng)快照相對(duì)應(yīng)，圖14中的負(fù)視速度同相軸就是面波在掌子面轉(zhuǎn)換為橫波傳播到接收排列所形成的同相軸SRSR波。

圖14 實(shí)例時(shí)間記錄Fig.14 Profiles of the example

以上分析表明，有限元數(shù)值方法可以有效地模擬隧道復(fù)雜地質(zhì)條件下的全波場(chǎng)的激發(fā)傳播過(guò)程。

4 結(jié)論

1）由于隧道地震預(yù)報(bào)的波場(chǎng)傳播主要是在巖體中進(jìn)行，相比于常規(guī)地面地震勘探，具有更好的彈性參數(shù)。對(duì)隧道地震波場(chǎng)的模擬而言，將質(zhì)量阻尼系數(shù)α和阻尼剛度系數(shù)β都設(shè)置為0是比較合理的近似，能滿足數(shù)值模擬的需要；若模擬地表地震勘探，則需要考慮介質(zhì)阻尼。

2）針對(duì)網(wǎng)格長(zhǎng)度，為滿足計(jì)算要求，根據(jù)網(wǎng)格長(zhǎng)度不大于（1/π～1/8）λmin的原則，針對(duì)主頻為300 Hz的震源而言，取網(wǎng)格長(zhǎng)度為1m就可以消除頻散效應(yīng)的影響。

3）加載黏彈性邊界條件可以很好地吸收邊界反射，并且其有良好的魯棒性。

4）算例表明，采用大型有限元軟件可以有效地模擬隧道復(fù)雜地質(zhì)條件下的全波場(chǎng)的激發(fā)傳播過(guò)程。

5）本文只是對(duì)二維隧道地震波場(chǎng)進(jìn)行了數(shù)值模擬研究，相比于實(shí)際需求，開(kāi)展三維隧道地震波場(chǎng)的數(shù)值模擬研究更具實(shí)際意義。

本文采用數(shù)值模擬方法雖然只研究了隧道地震反射波場(chǎng)的傳播規(guī)律，根據(jù)不同模型，它也可以模擬地震勘探領(lǐng)域中其他地震波場(chǎng)的傳播規(guī)律，這為地震波場(chǎng)的數(shù)值模擬提供了一種解決手段。

（References）：

［1］魯光銀.隧道地質(zhì)災(zāi)害反射波法探測(cè)數(shù)值模擬及圍巖f－Ahp分級(jí)研究［D］.長(zhǎng)沙：中南大學(xué)，2009.

Lu Guangyin.Numerical Simulation for Tunnel Geological Hazard Detection by Reflected Wave Method and Rock Dynmanic Classification by f-Ahp Method Research［D］.Changsha：South China University，2009.

［2］劉志剛，劉秀峰.TSP隧道地震勘探在隧道隧洞超前預(yù)報(bào)中的應(yīng)用與發(fā)展［J］.巖石力學(xué)與工程學(xué)報(bào)，2003，22（8）：1399-1402.

Liu Zhigang，Liu Xiufeng.TSP Application and Development in Tunnel Lead Forecast［J］.Chinese Journal of Rock Mechanics and Engineering，2003，22（8）：1399-1402.

［3］Hanson D R，Haramy K Y，Neil D M.Seismic Tomography Applied to Site Characterization［M］.Denver：American Society of Civil Engineers，2000.

［4］佘德平.波場(chǎng)數(shù)值模擬技術(shù)［J］.勘探地球物理進(jìn)展，2004，27（1）：16-21.She Deping.On Wave Field Numerical Modeling Techniques［J］.Progress in Exploration Geophysics，2004，27（1）：16-21.

［5］Boore D M A.Note on the Effect of Simple Topography on Seismic SH Waves［J］.Bulletin of the Seismological Society of America，1972，62（1）：275-284.

［6］Carcione M，Kosloff，Kosloff R.Wave Propagation Simulation in a Linear Viscoelastic Medium［J］.Geophysical Journal，1988，95（3）：597-611.

［7］Talezer H，Carcione J M，Kosloff D.An Accurate and Efficient Scheme for Wave-Propagation in Linear Viscoelastic Media［J］.Geophysics，1990，55（10）：1366-1379.

［8］Robertsson.Numerical Free-Surface Condition for Elastic Viscoelastic Finite-Difference Modeling in the Presence of Topography［J］.Geophysics，1996，61（6）：1921-1934.

［9］Graves R W.Simulating Seismic Wave Propagation in 3DElastic Media Using Staggered-Grid Finite Differences［J］.Bulletin of the Seismological Society of America，1996，86（4）：1091-1106.

［10］Carcione J，Quiroga-Goode G，Cavallini F.Wavefronts in Dissipative Anisotropic Media：Comparison of the Plane-Wave Theory with Numerical Modeling［J］.Geophysics，1996，61（3）：857-861.

［11］Ozdenvar T，McMechan G.Algorithms for Staggered-Grid Computations for Poroelastic，Elastic，Acoustic，and Scalar Wave Equations［J］.Geophysical Prospecting，1997，45（3）：403-420.

［12］Carcione J，Helle H.Numerical Solution of the Poroviscoelastic Wave Equation on a Staggered Mesh［J］.Journal of Computational Physics，1999，154（2）：520-527.

［13］何樵登，張中杰.3D橫向各向同性介質(zhì)中的體波［J］.石油地球物理勘探，1990，25（2）：137-147.

He Qiaodeng，Zhang Zhongjie.3DTransversely Isotropic Medium Body Wave［J］.Oil Geophysics Exploration，1990，25（2）：137-147.

［14］王延光，李振春，穆志平，等.一種適于強(qiáng)橫向變速的高階差分正演模擬方法［J］.石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2002，26（5）：37-39.

Wang Yanguang，Li Zhenchun，Mu Zhiping，et al.A High-Order Finite-Difference Forward-Modeling Fit for Strong Lateral Velocity Variation［J］.Journal of the University of Petroleum，China，2002，26（5）：37-39.

［15］裴正林.地震波標(biāo)量方程的小波自適應(yīng)網(wǎng)格有限差分法數(shù)值模擬［J］.石油地球物理勘探，2005，40（3）：267-272.

Pei Zhenglin.Wavelet Scalar Adaptive Mesh Seismic Wave Equation Finite Difference Numerical Simulation［J］.Oil Geophysics Exploration，2005，40（3）：267-272.

［16］董良國(guó).復(fù)雜地表?xiàng)l件下地震波傳播數(shù)值模擬［J］.勘探地球物理進(jìn)展，2005，28（3）：187-194.

Dong Liangguo.Complex Surface Conditions Numerical Simulation of Seismic Wave Propagation［J］.Progress in Exploration Geophysics，2005，28（3）：187-194.

［17］戴志陽(yáng)，孫建國(guó)，查顯杰.地震波場(chǎng)模擬中的褶積微分算子法［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2006，35（4）：520-524.

Dai Zhiyang，Sun Jianguo，Zha Xianjie.Seismic Wave Field Modeling with Convolutional Differentiator Algorithm［J］.Journal of Jilin University：Earth Science Edition，2006，35（4）：520-524.

［18］周曉華，陳祖斌，曾曉獻(xiàn)，等.交錯(cuò)網(wǎng)格有限差分法模擬微動(dòng)信號(hào)［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2012，42（3）：852-857.

Zhou Xiaohua，Chen Zubin，Zeng Xiaoxian，et al.Simulation of Microtremor Using Staggered-Grid Finite Difference Method［J］.Journal of Jilin University：Earth Science Edition，2012，42（3）：852-857.

［19］Bohlen T.Parallel 3-D Viscoelastic Finite Difference Seismic Modelling［J］.Computers ＆ Geosciences，2002，28（8）：887-899.

［20］Bohlen T，Lorang U，Rabbel W，et al.Rayleigh-to-Shear Wave Conversion at the Tunnel Face：From 3D-FD Modeling to Ahead-of-Drill Exploration［J］.Geophysics，2007，72（6）：67-79.

［21］Liu Jingbo，Yan Qiushi，Wu Jun.Analysis of Blast Wave Propagation Inside Tunnel［J］.Transactions of Tianjin University，2008，14（5）：358-362.

［22］Lu Guangyin，Han Xuli，Zhu Ziqiang.A Staggered-Grid High-Order Finite Difference Numerical Simulation for Tunnel Seismic Prediction Ahead of Construction［J］.Engineering Computation，2009，23（2）：225-228.

［23］牛濱華.半空間均勻各向同性單相固體彈性介質(zhì)與地震波傳播［M］.北京：地質(zhì)出版社，2005：200-205.

Niu Binhua.Single-Phase Half-Space Homogeneous Isotropic Elastic Solid Medium and Seismic Wave Propagation［M］.Beijing：Geological Publishing House，2005：200-205.

［24］Cheung Y K，Jin W G，Zienkiewicz O C.Direct Solution Procedure for Solution of Harmonic Problems Using Complete，Non-Singular，Trefftz Functions［J］.Communications in Applied Numerical Methods，1989，5（3）：159-169.

［25］宗福開(kāi).波傳播問(wèn)題中有限元分析的頻散特性及離散化準(zhǔn)則［J］.爆炸與沖擊，1984，4（4）：16-23.

Zong Fukai. Frequency-Dispersion Characteristics and Dicretization of the Finite Element Analysis in Wave Propagation Problems［J］.Explosion and Shock Waves，1984，4（4）：16-23.

［26］Kausel E.Local Transmitting Boundaries［J］.Journal of Engineering Mechanics，1988，114（6）：1011-1027.

［27］Costantino C，Miller C，Lufrano LA.Soil Structure Interaction Parameters from Finite Element Analysis［J］.Nuclear Engineering and Design，1976，38（2）：289-302.

［28］Miller T F，Schmidt F W.Use of a Pressure Weighted Interpolation Method for the Solution of the Incompressible Navier-Stokes Equations on a Nonstaggered Grid System［J］.Numerical Heat Transfer：Part A：Applications，1988，14（2）：213-233.

［29］Lysmer J，Kuhlemeyer R L.Finite Dynamic Model for Infinite Media［J］.Journal of Engineering Mechanics，1969，95（4）：859-877.

［30］Novak M，Hindy A.Seismic Analysis of Underground Tubular Structures［C］//Proceedings of the 7th World Conference on Earthquake Engineering.Istanbul：Ankara Natl Common Earthquake Eng，1980：287-294.

［31］Deeks J，Randolph M.Axisymmetrical Time-Domain Transmitting Boundaries［J］.Journal of Engineering Mechanics-Asce，1994，120（1）：25-42.

［32］Wolf J P.A Comparison of Time-Domain Transmitting Boundaries［J］.Earthquake Engineering ＆Structural Dynamics，1986，14（4）：655-673.

［33］劉晶波，谷音，杜義欣.一致粘彈性人工邊界及粘彈性邊界單 元 ［J］.巖 土 工 程 學(xué) 報(bào)，2006，28（9）：1070-1076.

Liu Jingbo，Gu Yin，Du Yixin.Consistent Viscous-Spring Artificial Boundaries and Viscous-Spring Boundary Elements［J］.Chinese Journal of Geotechnical Engineering，2006，28（9）：1070-1075.

Wavefield Modeling Based on the Finite Element Method for the Tunnel Seismic Prediction

Wang Zhaoling1，2，Liu Zhengping2，Huang Yunyan1，Huang Xianbin1

1.CollegeofCivilEngineering，SichuanAgricultureUniversity，Chengdu611830，China
2.CollegeofCivilEngineering，SouthwestJiaotongUniversity，Chengdu610031，China

In the field of exploration geophysics，the software packages that can fast and efficiently simulate 2Dand 3Dseismic full-wavefield，including simultaneously simulate P-wave，S-wave，and surface wave are not very much.In the structure mechanics，however，many sophisticated commercial software packages based on the finite element method including ANSYS and PLAXIS for multi-wave and multi-component have been employed for the solutions of the elastic constitutive equations.Compared to the complete elastic wave equations，the elastic ones include a damping term，which is a first-order differential term time.Therefore，by adjusting the parameter of the medium，that is，the Rayleigh coefficient to make the damping to zero，the elastic constitutive equations will degenerate into full elastic wave equation.The wavefields propagating in the rocks of a tunnel are free of the disturbance of low velocity resulting from the cover layer in the conventional seismic exploration.We simulated the wavefields of the tunnel seismic prediction using the ANSYS software，and researched the dispersion，damping and absorbing conditions etc.in the numerical modeling.For medium absorption，we compared time records and snapshots of the wavefields in the numerical simulation for the cases with and without the damping.In the tunnel seismic prediction，setting the damping to zero is a reasonable assumption.By comparing the relation of different grid length and dispersion in the numerical example，we find that the dispersion disappears when the mesh size is smaller than the wavelength of 1/π.By introducing a viscoelastic boundary condition deduced from the wave equation，boundary reflections can be effectively absorbed.Finally，an numerical example of the tunnel seismic prediction shows that the usage of the ANSYS software can effectively simulate the propagation of waves in the complicated geological conditions.

numerical modeling；ANSYS；FEM；tunnel seismic prediction；seismic wavefield；damping；dispersion；boundary condition

10.13278/j.cnki.jjuese.201404305

P631.4

A

王朝令，劉爭(zhēng)平，黃云艷，等.隧道地震預(yù)報(bào)波場(chǎng)的有限元數(shù)值模擬.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2014，44（4）：1369-1381.

10.13278/j.cnki.jjuese.201404305.

Wang Zhaoling，Liu Zhengping，Huang Yunyan，et al.Wavefield Modeling Based on the Finite Element Method for the Tunnel Seismic Prediction.Journal of Jilin University：Earth Science Edition，2014，44（4）：1369-1381.doi：10.13278/j.cnki.jjuese.201404305.

2013-12-12

國(guó)家自然科學(xué)基金項(xiàng)目（40874051）

王朝令（1980—，男，講師，博士，主要從事地震波場(chǎng)數(shù)值模擬、隧道超前預(yù)報(bào)預(yù)報(bào)研究，E-mail：wong81010＠gmail.com。