關于“用單元降剛法探索中國大陸強震遠距離跳遷及主體活動區(qū)域轉(zhuǎn)移”的討論——橫向各向同性“殺傷單元”才是更好的途徑

董培育，石耀霖＊

1中國科學院計算地球動力學重點實驗室，北京 100049

2中國科學院大學地球科學學院，北京 100049

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

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

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

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

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

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

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

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

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

圖3 文獻［1］單元降剛以后計算斷層發(fā)生后與發(fā)生前位移和應力變化量（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］

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

圖6 用橫向各向同性介質(zhì)消減G21模擬斷層滑動的計算結果（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

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

（

）土工基礎，2007，21（3）：54－56.

［1］ 楊樹新，陸遠忠，陳連旺等.用單元降剛法探索中國大陸強震遠距離跳遷及主體活動區(qū)域轉(zhuǎn)移.地球物理學報，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］ 陸遠忠，葉金鐸，蔣淳等.中國強震前兆地震活動圖像機理的三維數(shù)值模擬研究.地球物理學報，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］ 謝海國，章廣成，楊昌斌.邊坡分級開挖過程的動態(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模擬斷層滑動的計算結果（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.