董培育，石耀霖＊

1中國(guó)科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100049

2中國(guó)科學(xué)院大學(xué)地球科學(xué)學(xué)院，北京 100049

在《地球物理學(xué)報(bào)》第55卷第1期的“用單元降剛法探索中國(guó)大陸強(qiáng)震遠(yuǎn)距離跳遷及主體活動(dòng)區(qū)域轉(zhuǎn)移”［1］一文中提出：震源區(qū)部分喪失承載能力是引起強(qiáng)震遠(yuǎn)距離跳遷及主體活動(dòng)區(qū)域轉(zhuǎn)移的重要原因，而在模擬計(jì)算中采用顯著降低模型中當(dāng)前包含最大應(yīng)力單元的一組單元的彈性模量的方法，計(jì)算模擬強(qiáng)震斷層活動(dòng)引起的應(yīng)力場(chǎng)的變化調(diào)整.調(diào)整后的應(yīng)力場(chǎng)中再次尋找包含最高應(yīng)力單元的單元組，作為后續(xù)強(qiáng)震的發(fā)震位置，再次顯著降低這組單元的彈性模量，又引起新的應(yīng)力應(yīng)變場(chǎng)的調(diào)整和應(yīng)力集中轉(zhuǎn)移.與此同時(shí)，上一次被降低了彈性模量的單元組恢復(fù)性增加其彈性模量數(shù)值，模擬斷層的愈合.

在有限元計(jì)算中，研究張裂隙破裂，或工程中的開挖造成的應(yīng)力變化，利用“殺死單元”（kill element）——即假定裂開或挖走的單元楊氏模量降為0，是一種得到廣泛應(yīng)用的方法［2－3］.但地震剪切破裂時(shí)，雖然有應(yīng)力降，但應(yīng)力并不降為0，斷層仍具有一定的抗剪強(qiáng)度，因此如何將殺死單元的思路應(yīng)用于地震斷層，是一個(gè)有意義的課題，文獻(xiàn)［1］的這種方法是一次非常有益的嘗試探索.通過文獻(xiàn)［1］“殺傷”而非“殺死”的方法，實(shí)現(xiàn)使震源區(qū)域部分喪失承載能力，產(chǎn)生應(yīng)力降［1，4－5］.但是該文這種簡(jiǎn)單的減小彈性模量的方法也還存在嚴(yán)重的問題.根據(jù)各向同性下的虎克定律，如果泊松比ν不變，彈性模量E和剪切模量G是相互聯(lián)系的，降低了斷層帶內(nèi)的彈性模量E，剪切模量G＝E／（2（1＋ν））也會(huì)同時(shí)降低；或者說降低了剪切模量G，彈性模量E也會(huì)按比例降低.用降低剪切模量G模擬斷層剪切破裂是一種看似可行的思路，但同時(shí)彈性模量E降低，這不能正確模擬斷層在正應(yīng)力下的作用.文獻(xiàn)［1］中提到“單元的泊松比也做相應(yīng)的變化，單元的剪切模量和體積模量合理.E的原值是80GPa，減弱到17.8GPa”.如果令G保持不變，ν＝E／（2G）－1，那么大幅度減弱E之后的ν變?yōu)樨?fù)值，雖然聚合物泡沫材料發(fā)現(xiàn)有泊松比可以是負(fù)值的材料［6］，但是在自然條件下的巖石介質(zhì)中負(fù)值的泊松比基本沒有.所以即便是泊松比取值在合理范圍內(nèi)相應(yīng)調(diào)整，總是會(huì)在減弱E的同時(shí)，G也減弱，或減弱G同時(shí)，E也減弱.

我們利用一個(gè)簡(jiǎn)單的二維的平面應(yīng)力問題的例子來探討這個(gè)問題：假設(shè)初始狀態(tài)沒有斷層的存在，整個(gè)區(qū)域都是均勻各向同性的.在邊界荷載的作用下，發(fā)生剪切破裂，產(chǎn)生走滑斷層，如圖1所示，假設(shè)斷層長(zhǎng)1km，厚度為2m（灰色區(qū)域?yàn)閿鄬?，斷層夸大?邊界條件為：左邊界，垂直向下的位移2m，水平向右的位移1m；右邊界，垂直向上的位移2m，水平向左的位移1m；上、下邊界的條件相同，均為垂直方向自由，水平方向從左到右位移線性變化.這樣的邊界條件既有左旋剪切的作用，也有東西向垂直于斷層的正壓應(yīng)力的作用.

圖1 計(jì)算采用的含有斷層的模型Fig.1 Calculation model which contains fault

把上述邊界條件下的均勻各向同性不含斷層的二維模型的位移和應(yīng)力作為斷裂前的初始狀態(tài)，與后面斷層發(fā)生滑動(dòng)后的狀況比較.考慮斷層發(fā)生時(shí)，斷層帶內(nèi)的物質(zhì)采用文獻(xiàn)［1］的方法，減弱其彈性模量E，計(jì)算其在上述邊界條件下的位移和應(yīng)力狀態(tài)（圖2）.計(jì)算值減去初始狀態(tài)的值得到的差值，反映的即是發(fā)生地震前后的狀態(tài)變化（圖3），這一差值可與Okata的解析解（圖4）進(jìn)行比較.

圖2中雖然反映出斷層走滑對(duì)y方向位移和剪應(yīng)力的影響，但是由于斷層的楊氏模量被人為地降低，導(dǎo)致x方向位移出現(xiàn)了不合理的狀況，斷層西面和東面的巖體均向斷層移動(dòng)（圖2a），造成了斷層?xùn)|西部廣大區(qū)域的東西向壓應(yīng)力下降（圖2d）.而走滑斷層發(fā)生應(yīng)該主要降低斷層的剪應(yīng)力，不應(yīng)該對(duì)斷層上正應(yīng)力有如此巨大影響的.

這一問題在圖3和圖4的比較中看得更為明顯.圖3是用文獻(xiàn)［1］的降剛法計(jì)算的地震前后變化，圖4是在斷層在相應(yīng)位錯(cuò)量下用Okata公式計(jì)算的地震后位移和應(yīng)力變化量.圖4d顯示Okada解地震前后正應(yīng)力σxx幾乎沒有變化，而圖3d中σxx卻有兆帕量級(jí)的降低.這主要在于文獻(xiàn)［1］降剛法使斷層?xùn)|西兩側(cè)廣大區(qū)域在x方向發(fā)生了實(shí)際不應(yīng)該有的位移（圖3a），這也導(dǎo)致剪應(yīng)力降與Okada解有顯著差別（圖3c與圖4c比較）.

我們建議采用橫向各向同性的“殺傷單元”（wounded element）來處理斷層破裂問題，即把達(dá)到破裂的斷層單元的力學(xué)性質(zhì)改為弱化了的橫向各向同性物質(zhì)，平行于斷層剪切面的剪切模量降低，從而模擬斷層發(fā)生時(shí)的錯(cuò)動(dòng)和應(yīng)力降；但垂直于斷層面的楊氏模量不變，模擬法向應(yīng)力基本不變的事實(shí).橫向各向同性物質(zhì)，如圖5所示，假設(shè)1軸的方向?yàn)閷?duì)稱軸，橫向（即2－3平面內(nèi)）為各向同性的，但是沿著1軸的法向變形與2或3軸方向可以不同（Ei代表沿著i軸的彈性模量），1－2、1－3平面內(nèi)剪切變形與2－3平面內(nèi)剪切變形也可以不同（G21與G31相同，但與G23不同），νij代表i方向上受力時(shí)，與它垂直的j方向上相應(yīng)變形的泊松比；Gij代表法向量為i的平面上沿著j軸方向的剪切模量.這時(shí)共有5個(gè)獨(dú)立參量，本構(gòu)關(guān)系見圖5［7］.這里由于我們采用二維平面應(yīng)力問題以1－2平面作為剖面，則所需的獨(dú)立彈性參數(shù)進(jìn)一步減小為4個(gè)，即E1，E2，ν21，G21.

圖2 文獻(xiàn)［1］單元降剛以后計(jì)算的位移和應(yīng)力（a）x方向（垂直于斷層）的位移Ux；（b）y方向（平行于斷層）位移Uy；（c）剪應(yīng)力σxy；（d）正應(yīng)力σxx.Fig.2 The displacement and stress distribution after reducing the Yong′s modulus of fault elements with the method used in Ref.［1］（a）Displacement Uxof the direction x（vertical to the fault）；（b）Displacement Uyof the direction y（parallel to the fault）；（c）Shear stressσxy；（d）Normal stressσxx.

在模擬斷層發(fā)生時(shí)，背景單元仍然為各向同性，我們僅僅降低斷層單元的剪切模量G21，垂直于斷層方向的楊氏模量E1并不需要變化，即E1＝E2，均取原來各向同性巖石的E的值，泊松比ν21也取原來各向同性巖石的值ν0.然后進(jìn)行同樣的計(jì)算，得到斷層的滑動(dòng)量，如圖6所示，斷層發(fā)生后應(yīng)力狀態(tài)的變化見圖7.把圖7與給定類似斷層錯(cuò)動(dòng)，根據(jù)Okada給出的解析解［8］，利用汪榮江的開源計(jì)算程序［8－9］來計(jì)算的圖4比較，可以看出位移和應(yīng)力變化量，不但空間分布圖像相似，而且各點(diǎn)計(jì)算值大小也很相近，各相應(yīng)圖件僅有微小差別.應(yīng)該指出，Okata的解適用于厚度趨于0的線狀斷層，我們的斷層有有限厚度，二者不會(huì)完全相同.另外目前我們計(jì)算中把區(qū)域簡(jiǎn)單按矩形網(wǎng)格劃分，因此對(duì)斷層尖端附近的應(yīng)力集中區(qū)域，網(wǎng)格不夠密，數(shù)值計(jì)算誤差會(huì)比較大，如果改進(jìn)單元?jiǎng)澐?，采用非結(jié)構(gòu)化網(wǎng)格對(duì)斷層端點(diǎn)附近加密，應(yīng)該會(huì)取得更佳結(jié)果.但目前的計(jì)算已經(jīng)足以說明，我們建議采用的橫向各向同性“殺傷單元”減小G21的方法，能夠較好地反映斷層滑動(dòng)時(shí)的應(yīng)力變化.

圖3 文獻(xiàn)［1］單元降剛以后計(jì)算斷層發(fā)生后與發(fā)生前位移和應(yīng)力變化量（a）ΔUx；（b）ΔUy；（c）Δσxy；（d）Δσxx.Fig.3 The change of displacement and stress before and after fracture with the method used in Ref.［1］

總之，強(qiáng)震跳遷過程的地球動(dòng)力學(xué)模擬是非常復(fù)雜困難的課題，文獻(xiàn)［1］的基本思路是有價(jià)值的，采用單元降剛，即顯著降低其彈性模量的方法，是一種新的探索，有助于了解大范圍的兆帕量級(jí)的應(yīng)力場(chǎng)的調(diào)整，解釋后續(xù)強(qiáng)震遠(yuǎn)距離跳遷的問題.但是在存在垂直斷層的正應(yīng)力情況下，簡(jiǎn)單通過減弱其彈性模量E的具體方法處理剪切斷層存在嚴(yán)重的問題.我們建議采用的橫向各向同性“殺傷單元”減小G21的方法，能夠較好地反映斷層滑動(dòng)時(shí)的應(yīng)力變化.而隨時(shí)間演變逐步增加G21，則能夠反映斷層的愈合.

圖6 用橫向各向同性介質(zhì)消減G21模擬斷層滑動(dòng)的計(jì)算結(jié)果（a）Ux；（b）Uy；（c）σxy；（d）σxx.Fig.6 The result of simulation fault slipping with the method of reducing G21 in the transverse isotropic material

致 謝 本文得到程惠紅和尹鳳玲的有益建議，特此感謝.

（

）土工基礎(chǔ)，2007，21（3）：54－56.

［1］ 楊樹新，陸遠(yuǎn)忠，陳連旺等.用單元降剛法探索中國(guó)大陸強(qiáng)震遠(yuǎn)距離跳遷及主體活動(dòng)區(qū)域轉(zhuǎn)移.地球物理學(xué)報(bào)，2012，55（1）：105－116.

Yang S X，Lu Y Z，Chen L W，et al.The mechanism of long－distance jumping and the migration of main active areas for strong earthquakes occurred in the Chinese continent.

ChineseJ.Geophys.（in Chinese），2012，55（1）：105－116.

［2］ 閻波，方云.基于ANSYS的邊坡開挖模擬.山西建筑，2008，34（1）：93－94.

Xie H G，Zhang G C，Yang C B.Achievement of dynamic simulation of excavation of slope in a hydropower station.SoilEng.andFoundation（in Chinese），2007，21（3）：54－56.

［4］ Lu Y，Yang S，Chen L，et al.Mechanism of the spatial distribution and migration of the strong earthquakes in China inferred from numerical simulation.JournalofAsianEarth Sciences，2011，40（4）：990－1001.

［5］ 陸遠(yuǎn)忠，葉金鐸，蔣淳等.中國(guó)強(qiáng)震前兆地震活動(dòng)圖像機(jī)理的三維數(shù)值模擬研究.地球物理學(xué)報(bào)，2007，50（2）：499－508.

Yan B，F(xiàn)ang Y.Simulation of slope excavation based upon ANSYS.ShanxiArchitecture（in Chinese），2008，34（1）：93－94.

［3］ 謝海國(guó)，章廣成，楊昌斌.邊坡分級(jí)開挖過程的動(dòng)態(tài)模擬.

Lu Y Z，Ye J D，Jiang C，et al.3Dnumerial simulation on the mechanism or precursory seismicity pattern before strong earthquake in China.ChineseJ.Geophys.（in Chinese），2007，50（2）：499－508.

［6］ Lakes R S.Foam structures with a negative Poisson′s ratio.Science，1987，235（4792）：1038－1040.

圖7 用橫向各向同性介質(zhì)消減G21模擬斷層滑動(dòng)的計(jì)算結(jié)果（a）ΔUx；（b）ΔUy；（c）Δσxy；（d）Δσxx.Fig.7 The result of simulation fault slipping with the method of reducing G21 in the transverse isotropic material

［7］ Slawinski M A.Waves and Rays in Elastic Continua：2nd Ed.World Scientific，2007.

［8］ Okada Y.Surface deformation due to shear and tensile faults in a half－space.BulletinoftheSeismologicalSocietyof America，1985，75（4）：1135－1154.

［9］ Wang R，Martin F L，Roth F.Computation of deformation induced by earthquakes in a multi－layered elastic crust－FORTRAN programs EDGRN／EDCMP.Computers＆Geosciences，2002，29（2）：195－207.