斷層系統摩擦動力學行為的有限元模擬分析

邢會林，郭志偉，王建超，張熔鑫，劉駿標，姚琪

1 深海圈層與地球系統前沿科學中心,海底科學與探測技術教育部重點實驗室, 中國海洋大學海洋高等研究院/海洋地球科學學院,青島 266100 2 青島海洋科學與技術國家實驗室 海洋礦產資源評價與探測技術功能實驗室,青島 266100 3 中國海洋大學海底科學與工程計算國際中心,青島 266100 4 中國地震局地震預測研究所，北京 100036

0 引言

地震預報是地球科學界的世界性難題，通過基于地球物理場觀測、巖石物理實驗和定量化的物理規律的數值模擬，研究地震孕育、破裂的物理機理，指導地震危險性的有效評估，最終實現基于物理的地震數值預報已成為人們的共識(石耀霖等,2013,2018;Barbot et al.,2012;Segall,2012;劉啟元和吳建春，2003).王仁等(1980，1982)曾嘗試利用有限元方法模擬分析華北地區地震序列及遷移規律，隨后朱岳清等(1988)和蔡永恩等(1999)曾分別就唐山地震孕震和震源動力過程數值分析進行了初步的有益嘗試.石耀霖等(2018)在總結統一的美國加州地震破裂預測系統基礎上，提出了我國地震數值預報路線圖的設想，吹響我國地震數值預報的起床號.

世界上80%以上的大地震發生在板塊邊界及板塊內部斷層上，被認為是斷層或板塊邊界上的一種摩擦失穩行為(Brace and Byerlee,1966).中國大陸具有以板內塊體為單元的構造變形特征，且強震多發生在活動塊體的邊界/斷層上.地震的孕育和發生過程也就是斷層閉鎖、解鎖和滑動過程，其中斷層面上的接觸摩擦失穩行為起著決定性的作用，為此國內外學者對此進行了大量的研究.一般而言，含斷層的有限元分析主要有兩種不同的處理方式：一種把斷層當作某種特殊材料(如弱材料)的連續介質(如王仁等,1982;陳化然等,2004;王凱英和馬瑾,2004;Zhang et al.,2009;Luo and Liu,2010;He et al.,2011;Liu et al.,2016;Li et al.,2017;Sun and Liu,2018)，另一種把斷層作為非連續的摩擦接觸界面(如Goodman,1968;Melosh and Williams,1989;Huang et al.,1997;Xing and Makinouchi,2002a,b;朱守彪和張培震,2009;Xing and Xu,2011;姚琪等,2012a,b,2018a;吳萍萍,2014;范桃園等,2014).兩者都可計算相關的應力及變形行為，但后者在理論上或可更加真實地模擬斷層的動力學行為.基于摩擦接觸的斷層有限元分析，摩擦面法向上多采用罰函數、拉格朗日乘子法或其擴展方法處理.‘Slippery node’方法(Melosh and Williams,1989)和雙節點節理模型(Goodman,1968)一樣，實質上采用了類似拉格朗日乘子法的處理都引入了新的未知變量，且其特殊處理也僅適用于二維模型計算(Huang et al.,1997).目前含斷層的有限元分析多集中在長時間變形分析或庫侖應力場計算，斷層的摩擦系數多取為常數，摩擦失穩細節基本不考慮，且多為二維單一斷層模型(Huang et al.,1997;朱守彪和張培震,2009;范桃園等,2014;吳萍萍等,2014).然而在現實中，大量的地震觀測和研究表明地震孕育、破裂并非單一斷層的孤立行為，而是以系統的形式存在.一方面，地震可能是由若干斷層級聯式破裂的，例如2016年新西蘭地震破裂了10條左右斷層.另一方面，同震破裂、震后應力調整使得周圍斷層面應力增加或減弱，激發大量余震或延遲觸發大地震.斷層系統的有限元數值模擬是了解它們的動力學和復雜的系統行為的有力工具，但是以往的研究由于觀測資料少、斷層系統精細網格難以生成等原因導致所用的模型過于簡單粗糙，加之斷層系統數值模擬技術不夠穩定和成熟等客觀條件的限制，對地震過程的數值模擬研究還十分有限.

雖然也有一些內部和商業軟件用于模擬一些與地震有關的摩擦現象，但含斷層相互作用的地殼動力學系統仍然是一個非常復雜的問題，因為沿斷層的黏著-滑移摩擦失穩具有很強的非線性.目前，絕大多數基于摩擦接觸的斷層/地震行為的研究多采用機械工程行業的商業化軟件，如ANSYS(如吳萍萍等,2014)和ABAQUS(如范桃園等,2014)來模擬斷層行為，且摩擦系數多簡單地取為常數.但是，機械工程行業常采用潤滑劑降低界面摩擦系數，避免加工面發生劃痕損傷同時降低能耗，摩擦系數一般最大為0.1左右.而地震斷層上的靜態摩擦系數一般在0.6以上，而且摩擦系數在摩擦失穩/地震發生時會急速下降，比如到0.2，而非一個常數.由此導致的高度非線性問題使得目前的商業化軟件(如ANSYS、ABAQUS、MSC-MARC和ADINA等)計算穩定性遇到了巨大的挑戰.此外，有的研究者嘗試用動態顯式法來模擬地震的破裂過程，但是計算時間步長太小，對于地震的孕育過程及整個地震周期的模擬則需要太長的計算時間，應力計算精度也多有問題.表1參照較為成熟的典型變形過程有限元計算分析，比較了不同時間積分方案的異同，其規律對地震過程的模擬也基本適用.

表1 靜態隱式、靜態顯式和動態顯式的比較Table 1 Comparison of static implicit,static explicit and dynamic explicit formulations

在目前的有限元地震孕育過程模擬分析中，靜態隱式算法曾被嘗試應用于此類非線性接觸問題.表2對此進行了部分總結.有限元軟件靜態隱式算法通常采用Newton、Newton-Raphson或者修正的Newton-Raphson迭代法，其將控制方程作為一個非線性方程(組)，計算相應的位移增量來求解，迭代到平衡為止.然而Newton-Raphson等迭代算法的一個主要缺點是迭代可能會不收斂，因為對于非線性問題每步開始前難以給出足夠精確的初始解.更何況，為了模擬地震行為，還要處理斷層上的狀態變化(如接觸與脫離、黏著與滑移)以及滑動速率、滑移及接觸壓力等相關的摩擦系數變化等新的突變性非線性問題，而傳統商業化軟件面臨這些問題往往收斂性極差.雖然采用迭代平滑技術來減少劇烈變化，但收斂仍然無法得到保證，而且這可能會導致黏著-滑移失穩發生時刻(即地震發生時刻)不能被準確地捕捉到.

表2 準靜態摩擦接觸問題有限元軟件Table 2 Survey of some FEM software for quasi-static frictional contact problems

綜上所述，雖然通過實驗室試驗和現場觀測地震斷層摩擦本構關系的研究取得了一系列成果，提出了諸如滑動位移、速度-狀態相關等不同的非線性摩擦本構模型(Dieterich，1978，1979；Kato and Tullis,2001)，但是目前鮮有商業有限元程序可用于研究與地震孕育成核和破裂相關的非線性摩擦失穩、滑移動力學行為.現行代碼/軟件大多在接觸界面上使用了預設的摩擦系數定值，而這些摩擦系數值與可變條件無關，因此無法研究變形巖石對摩擦界面行為的影響，如黏著-滑移摩擦失穩和沿接觸界面的非線性滑動摩擦行為.雖然或許可以通過用戶子程序來實現地震過程模擬所需的部分功能，但是由于如上所述非線性問題的收斂性難題可能導致不收斂而無解；如果采用迭代平滑技術來減少劇烈變化從而滿足收斂條件，會導致黏著-滑移失穩(即地震發生的時刻)等關鍵突變時刻漏失而不能被準確捕捉到.

因此，Xing等(1998)、Xing和Makinouchi(2002a，2002c)提出了一種基于R-minmium自適應控制時間步長的靜態顯式時間積分方案，兼顧了靜態隱式與動態顯式的部分優點(見表1和表2中PANDAS)，可以精確捕捉黏著-滑移失穩等關鍵時刻；利用節點-點接觸單元，將其應用于可變形體之間的三維任意形狀非線性摩擦接觸問題的有限元模擬，進一步利用速率和狀態相關的摩擦定律分析斷層的非線性摩擦失穩行為，取得了較好的效果，具體總結如下.

1 摩擦接觸有限元算法

1.1 有限元公式

為了更好地描述地震過程尤其是局部大滑移大變形等問題，我們采用更新的拉格朗日公式來描述有限變形問題.結合物體間摩擦接觸條件，其平衡方程可以采用虛擬速度原理等效表示為(Xing et al.,1998，Xing and Makinouchi,2002a,b,c)：

(1)

(2)

(3)

當材料表現為率相關，也就是黏彈性或黏彈塑性時，材料的本構方程可以表述為

(4)

1.2 摩擦接觸的統一化本構方程

1.2.1 法向接觸力

采用罰函數方法處理法向接觸不可穿入的約束關系.當兩個點在接觸面的法線方向上的相對位移小于零時，即gn<0時發生接觸，其法向接觸力為：

fn=f·n=Engn(≠0僅限于gn<0),

(5)

式中，En為法向上的罰因子，gn為法向穿透量，f為接觸力，n為接觸體表面的外法線方向.

1.2.2 摩擦力

斷層面上的摩擦接觸關系是斷層行為變化的關鍵因素，斷層面的形狀、運動狀態、摩擦系數及其不均勻分布和演化等條件都可以通過斷層的摩擦接觸行為來表現，這為使用非線性摩擦有限元方法來進行區域地震危險性分析提供了可行性.

大多數構造地震都是由沿先存斷層或板塊界面的突然滑動引起的.大地震被認為是沿著斷層的黏著-滑移摩擦失穩的結果(Brace and Byerlee，1966).此后，基于實驗結果提出了幾種摩擦本構方程，目前非線性摩擦問題可歸納為以下兩種模式(Dieterich,1978，1979；Kato and Tullis,2001;Xing et al.,2006)，一種是滑動位移依賴型，一種是速度-狀態依賴型.這兩種模式都可以用來描述斷層摩擦過程中的某種現象，但是關注點不一樣.滑動依賴模式著重于破裂過程的描述，譬如從黏著狀態轉換到滑動狀態的動態過程，但忽略了速度效應；而速度-狀態依賴模式描述了滑動過程中狀態以及滑動速度的影響，但難以描述從黏著/閉鎖狀態到產生初始滑動/破裂的行為.其實從實驗結果也可以看出，兩種模式描述了斷層摩擦行為的不同階段.在模擬斷層行為和地震過程中，人們一直在試圖將兩種模式結合起來使用(Xing et al.,2006,2007).Xing等(Xing and Makinouchi,2002a,2002c,2003;Xing et al.,2006,2007)導出了統一化的數學公式表述速度相關的摩擦接觸中黏著和滑移不同的運動狀態，并導入有限元分析軟件PANDAS，計算在內外載荷作用下接觸面上摩擦接觸狀態(黏著和滑移兩種狀態)及變化過程，即斷層上閉鎖和解鎖滑移兩種不同運動狀態的變化特征，也就是地震及其周期的孕育、發生的過程(參見表2中PANDAS).

(6)

(黏著狀態,l=1,2)

(7a)

(7b)

1.2.3 速度-狀態摩擦本構關系

Dieterich(1978，1979)和Ruina(1983)提出速率和狀態相關摩擦關系式來描述斷層的摩擦行為，已得到廣泛應用；據此，與滑動速率相關的瞬態摩擦行為可表示為：

μ=μ0+φ+aln(V/Vref),

dφ/dt=-[(V/L)(φ+bln(V/Vref))].

(8)

把式(8)代入式(7)，也可以得到摩擦的類似材料本構式(4)的統一化表達式.此外，摩擦力是一個矢量，它也會隨著接觸面法線方向的變化而變化(Xing and Makinouchi,2002a,c;Xing et al.,2006).

1.3 時間積分方法

如前所述，為了應對摩擦接觸狀態、摩擦非線性以及材料非線性所引起的收斂性問題，我們采用了顯式時間積分方案.假設在足夠小的時間步長內，只要狀態不發生劇烈變化(積分點發生彈塑性變化、摩擦界面節點發生接觸脫離或者黏著-滑移狀態變化)，公式(1)中的所有變化率在從t到t+Δt的時間增量內都可以被認為是恒定的.我們采用R-Minimum方法來限制步長，以避免在一個增量步長時間內狀態發生劇烈變化，從而精確捕捉到狀態變化的時刻(Xing and Makinouchi,2002a,2002c).

因此，式(1)中的所有變化率皆可以用增量來代替:

(9)

最后，結合上述公式(5)(7)和(8)，公式(1)有限元離散后可寫為:

(K+Kf)Δu=ΔF+ΔFf,

(10)

式中，K和ΔF分別是整個物體的標準剛度矩陣和外載荷邊界增量；當采用率相關材料模型時，它們也包括了黏性項的貢獻(如，Xing and Zhang,2009;Xing and Wang,1997)；Kf和ΔFf是所有接觸單元的接觸剛度矩陣和接觸力增量，Δu是節點位移增量(Xing and Makinouchi,2002a,c).

2 斷層系統有限元計算網格生成

網格生成是進行有限元計算分析之前的關鍵一步，它是將研究對象離散為計算用網格(單元)，以便施加邊界條件進行有限元模擬計算.目前主要使用商業化圖形網格軟件(如HyperMesh和Patran)以及ABAQUS、ANSYS等有限元分析軟件中的前處理等功能模塊進行網格生成.通過這些軟件，可以將二維幾何模型離散為三角形或四邊形網格，或者將三維幾何模型離散為四面體或六面體網格.由于這些商業化軟件都是為機械和土木工程行業設計分析而開發的，一般通過確定的線、面和體等幾何信息生成網格.這樣的方法對于處理簡單形狀的地質模型尚可，但對于含復雜斷層系統的地質體，一般在地質勘探圖像中提取的相關斷層數據并用一系列點陣來表示，而非確定的線、面和體等傳統網格生成軟件所需要的幾何描述，因此含復雜斷層系統地質體的計算用有限元網格模型生成還面臨著較大的挑戰.

目前常用的計算網格生成算法包括以下四種：Delaunay算法(Baker,1989;Borouchaki and Lo,1995;Borouchaki et al.,2000)，八叉樹(四叉樹)法(Frey et al.,1994;Schneiders,1996)，推進波前法(AFT)(L?hner,1996a,b;Lee and Hobbs,1999)以及掃描法或映射網格法(表3).為了在三維空間中生成包含由大量點數據集定義的斷層網格，可以使用Delaunay算法或推進波前法自動生成四面體網格.但傳統的Delaunay算法在應對復雜斷層系統的有限元網格生成時并不理想，因為該方法并不能保證網格節點精確地插入邊界或斷層面處.也就是說，在生成的地質體網格中，斷層面的形狀和位置會發生變化，不能保證其原有的形狀及特征.一般說來，雖然四邊形和六面體網格比三角形和四面體網格具有更好的品質(Schneiders and Bünten,1995)，更適合于非線性有限元分析，但是對于三維地質體模型，特別是對于包含復雜斷層的大規模地質體模型，目前仍難以自動生成高質量的非連續六面體網格(Lee and Lee,2002;Xing et al.,2009).

表3 表面域網格劃分的主要算法及介紹Table 3 Main algorithms and introduction of surface mesh generation

2.1 基于斷層面的非連續四面體網格生成

當斷層系統的斷層面非常復雜而曲折時，用六面體網格來對地質體模型進行網格劃分不僅費時費力，而且往往只能將曲折界面簡化為一系列具有簡單幾何形狀的平面.為此，采用基于斷層面約束的四面體網格自動生成方法，可以保證其中斷層的幾何形狀、位置和整體網格的形狀和質量，并用于有限元計算.具體的操作過程如下：首先對斷裂帶數據進行邊界提取與約束，通過Delaunay與AFT算法相結合的方法建立空間三角形來描述斷層空間曲面；其次對面網格進行優化處理，尤其是對斷層曲面的交叉部分進行特殊識別及優化；接著為研究區生成一個由三角形網格組成的封閉外殼以完全包裹研究區斷層面；最后以生成的內部斷層及外殼的三角形面網格為約束條件，生成四面體網格(Liu and Xing,2016).這種方法的優點是可快速地自動生成含復雜斷層系統的四面體有限元網格模型.然而，全自動化網格生成帶來的缺點是生成的四面體網格單元和節點數量很多，網格質量尤其是含空間曲面斷層時較難控制；同時對標準的線性四面體單元而言有限元應力計算精度較低，還在進一步優化研究之中.

2.2 非連續六面體網格生成

掃描法是非連續六面體網格生成較為常用的一種方法，通過直接或者非直接方法生成全四邊形表面網格，然后以此按照預定的方向推進生成六面體形狀的網格.同時將斷層線推出的內表面拆分成兩個單元間的分界面，形成沿斷層的不連續摩擦接觸界面.這種方法可以簡單、快速地生成六面體網格，但生成的斷層表面傾角多為90°，比較適用于走滑斷層.這種方法的難點在于作為基礎面的全四邊形表面網格的生成上(Xing et al.,2009).Liu和Xing(2013)將三角形合成和波前推進法相結合，自動生成具有復雜線約束的高質量全四邊形網格.下面，以含復雜裂隙結構材料的數字圖像(圖1a)為例，來說明利用該方法實現斷層系統非連續六面體網格的生成.具體步驟如下：(1)斷層結構提取與修復：利用PANDAS軟件將圖1a中的斷層特征提取出來，并采用蟻群算法、DBSCAN聚類算法和Canny算子的混合圖像特征修復算法(張愉玲等，2021)，恢復圖像中缺失的斷層特征(圖1b)；(2)邊界平滑：由于數字圖像的性質，提取的邊界幾何形狀通常是鋸齒狀的，并影響生成的四邊形網格的質量，采用基于相位的邊界平滑方法來平滑鋸齒狀邊界，并對其進行離散，考慮到三維有限元模型進行應力加載和邊界條件設置的邊界效應，在斷層區的四周加了一個較大的矩形邊界，用于約束網格的生成和防止計算時邊界出現突變(圖1c).為了使斷層周圍的網格更加規則和精細，可以圍繞斷層生成一條或者幾條盡量平行于斷層的線作為新的約束線(圖1d)；(3)初始三角形網格生成：利用改進的二維Delaunay三角剖分法生成具有斷層線段約束的三角形網格(圖1e)；(4)四邊形網格生成：使用Catmull-Clark細分在每個三角形的質心處插入一個節點，將每個三角形分成三個四邊形(將三角形網格轉換為四邊形網格).然后從線約束開始，基于生成的全四邊形網格，通過推進波前法生成新的四邊形單元；(5)四邊形網格優化：通過優化四邊形網格的拓撲結構以減少不規則節點的數量，并平滑生成的網格，交替執行全四邊形網格的生成，直到生成高質量的網格(圖1f)；(6)三維有限元網格生成：以四邊形表面網格為基礎種子面，沿著用戶約定的方向推進掃描生成六面體形狀的網格，然后采用潛在接觸對方案對摩擦接觸面進行拆解分離和標識(圖1g).

圖1基于圖像生成的三維有限元模型(d,e,f中紅色線段表示斷層)(a)原始圖像；(b)提取和修復后的斷層；(c)平滑處理；(d)沿橙色斷層生成藍色等值線；(e)初始三角形網格生成；(f)四邊形網格生成；(g)六面體網格生成(橙色四邊形單元表示斷層面).Fig.1 3D finite element model generated based on rock images (red line segments in d,e and f represent faults)(a)The original image;(b)Extracted and repaired faults;(c)Smoothing;(d)Blue contour lines generated along the fault in orange;(e)Initial triangular mesh generation;(f)Quadrilateral mesh generation;(g)Hexahedral mesh generation (Orange quadrilateral grids indicate fault planes).

對于斷層在垂直方向上并不一定是豎直等特殊平面的情況，難以采用上述掃描法生成六面體網格，一般可采用映射分塊法劃分斷層系統模型.本文以圖2a所示的澳大利亞南部斷層系統三維幾何模型為例進行說明，其中部分斷層面傾角非90°且并不是直線/平面.為了將該三維斷層系統生成六面體網格并設定有限元模擬計算所必需的條件(即邊界條件和關于斷層的信息)，首先，在斷層曲線/面的約束下，將整個斷層系統幾何模型分解成多個簡單的可以網格化的六面體幾何塊體，其中斷層面傾角為90°的可以用上述掃描法進行網格劃分，其余的利用映射分塊法進行網格劃分(如圖2b)；這些分組塊各自生成六面體有限元網格后，在相鄰組塊體之間除斷層接觸面以外的公共邊或曲面處共享相同的網格單元；最后，用“焊接”的方式(即刪除重復點)或網格劃分后用黏著接觸算法將上述六面體網格組裝在一起(如圖3)，同時保持跨斷層的網格不連續，在進一步的有限元分析中將其作為摩擦接觸界面處理(Xing and Mora,2006;Xing et al.,2009).這種方法適用于可以將斷層系統用映射分塊法劃分的模型，但對于復雜斷層系統網格生成費時費力，且容易出錯.一旦能夠生成合理的網格，即可進行有限元計算分析.

圖2 南澳斷層系統(a)三維斷層系統幾何模型；(b)—(d)斷層系統的不同幾何組件/塊體及其六面體網格;(e)組裝在一起后生成的總體網格(修改自Xing and Mora,2006).Fig.2 The Southern Australian fault system(a)Geometry model of 3D fault system;(b)—(d)The hexahedron mesh generated from the different components/blocks of the fault system model;(e)The total mesh generated after assembling all the components together (Modified from Xing and Mora,2006).

3 摩擦接觸有限元在斷層系統分析中的應用

在上述的基礎上，Xing等提出并采用R-Minimum自適應地控制不同非線性情況下的時間步長方案以保證計算精度和效率，將不穩定的隱式迭代計算轉變為顯式計算，可以自適應的控制時間步長以穩定高效地模擬地震孕育(時間可以為幾十年甚至更長時間)和動態破裂(時間可能小到以分秒計)等地震過程及地震周期，采用基于R-Minimum的靜態-顯式方案開發了工業強度的摩擦接觸有限元軟件PANDAS(參見表2)，徹底解決了非線性斷層系統摩擦接觸的收斂性難題 (Xing and Makinouchi 2002a,b,c,2003;Xing et al.,2004,2006,2007).

目前，該軟件已經成功模擬重現了經典的巖石摩擦實驗，如三明治斷層模型(Xing and Makinouchi,2002b)、單曲(Xing and Makinouchi,2003;Xing et al.,2004)和多曲的非平直三維斷層模型(Xing et al.,2006).計算所得的接觸面上的摩擦力、應變和應變率的演化與實驗室發現的滑動弱化摩擦現象取得了很好的匹配.

針對現實中的多曲單斷層摩擦行為模擬，該軟件已經得到了很多的應用.Xing和Makinouchi(2002b),Xing等(2008)將這個方法應用到太平洋板塊在日本海深部向下俯沖的模擬中，實現了連續曲度變化的三維斷層面摩擦行為的收斂計算.Xing等(ACES會議報告，2006 5thACES,2008 6thACES,see http:∥www.aces.org.au)和朱守彪等(2008)利用該軟件計算了多曲度蘇門答臘俯沖帶的摩擦行為，模擬了2004年蘇門答臘地震的非線性孕震過程和地震破裂現象.姚琪等(2012a,b)利用該軟件模擬計算了汶川地震起始破裂處的斷層傾角和物性參數對強震的控制作用，證實了高傾角和上下盤巨大的巖性差異都能造成逆沖斷層超長的孕震時間，并且龍日壩斷裂分擔了巴顏喀拉塊體內部大部分的走滑分量，使得汶川地震起始破裂處以逆沖活動為主.姚琪等(2018b)利用該軟件模擬了具有多曲度的喜馬拉雅主逆沖帶的斷層破裂行為，用橫向分段來解釋了2015年尼泊爾地震復雜的地震破裂傳播特征.郭婷婷等(2015)采用PANDAS對單斷層與交叉斷層兩種斷層模型分別進行數值模擬計算對比，結合中國大陸雙震或震群型地震孕育發生的構造條件對共軛斷層系統模型的孕震與發震機理進行了討論與分析.

基于上述對非連續網格劃分的方法和有限元基本理論的研究，本文以圖1a中的復雜結構斷裂系統為例進行研究.采用2.2節生成的非連續六面體網格模型作為本次模擬計算的復雜斷層系統模型，三維有限元模型共有69595個節點，49356個八節點六面體(見圖1f).這里使用速率相關摩擦定律(公式(8))來描述沿斷層的摩擦行為，式中μ0=0.6，a=0.05，b=0.025;本文選取Vref作為整個過程的參考速度,同時假設模型中的材料具有相同的性質：密度ρ=2.60 g·cm-3，楊氏模量E=44.8 GPa，泊松比γ=0.22.模擬分析采用以下邊界條件：在模型上面，沿Y軸的負方向以恒定的速度Vref/10施壓，但Z方向位移固定；在底部，沿X、Y和Z三個方向固定；其他部分自由無邊界條件約束.整個模型使用六面體單元進行三維計算，初始應力設為0.圖3展示了不同階段的X方向速度場(a—d)和總速度場(e—h)在XY平面分布變化情況.在加載開始階段，數值模擬結果顯示速度場及其X方向速度分量圖表現為連續分布，斷層兩側的速度分布顯示還未發生明顯的滑動(圖3a和3e)；隨著加載的進行，圖3b和3f中與主加載方向(Y軸方向)夾角較小的斷層，由于其斷層面上切向摩擦力較大而法向正應力分量較小從而較易滑動，因此在較短的時間內由逐步閉鎖狀態進入局部滑動狀態,并開始表現出局部斷層區速度變化(圖4a).隨著持續加載，開始出現較大的斷層破裂活動區，此時破裂區仍主要分布在低角度的斷層部分(圖3c和3g，圖4b)；與主加載方向(Y軸方向)夾角較大的斷層，在整體變形及周邊斷層滑動的影響下，也會逐漸進入滑動甚至脫離狀態，導致斷層兩側發生較大的速度突變(見圖3d和3h,圖4c)，斷層系統的破裂范圍發生進一步的擴大，整個計算區域整體明顯向X軸的正負方向擠出.

圖3 X方向速度場(a—d)和總速度場(e—f)在不同加載階段的分布(單位：m·s-1).加載時間為100R/Vref其中(a)(e)R=0.001，(b)(f)R=0.003，(c)(g)R=0.014，(d)(h)R=0.06；深藍曲線代表斷層，(d)中黑色箭頭表示塊體擠出方向.Fig.3 The velocity component in X direction (a—d)and the total velocity magnitude (e—f)at different loading stages (unit:m·s-1).Loading time is 100R/Vref for (a)(e)R=0.001,(b)(f)R=0.003,(c)(g)R=0.014,(d)(h)R=0.06Dark blue curves represent faults,and black arrows in (d)represent the movement direction of blocks.

圖4 (a)—(c)分別是圖3(b)—(d)對應時刻間X方向速度分量的變化率深藍色曲線代表斷層.Fig.4 (a)—(c)are variation of velocity component in the X direction respectively between different moments in Fig.3 (b)—(d)Dark blue curves represent faults.

此外，該軟件已經在地震相關的復雜斷裂系統摩擦行為模擬中取得了多方應用，譬如Xing和Mora(2006)模擬了澳洲南部的復雜斷層系統中(圖2)斷層的相互作用，證明了該軟件能夠用于不規整、多曲面、多斷層、多接觸的三維非線性摩擦模擬，實現了計算的收斂，獲取應力場、應變場和斷層節點摩擦力的時空演化特征，從而以此推測復雜斷層系統中，多斷層段逐次破裂的時空進程，圖5顯示了其中等效應力變化的一個特寫鏡頭(Xing and Mora,2006).類似的，Xing等(2007)將美國南加州地區800 km×700 km×10 km范圍內的9條斷層進行了簡化、建模和非線性摩擦有限元計算，模擬獲取的等效應力場演化顯示了不同規模斷層的破裂、愈合、觸發、延遲的復雜過程.針對國內川滇地區龍門山斷裂、鮮水河—安寧河—小江斷裂系、大涼山斷裂等具有強震危險性的斷裂，Xing和Xu(2011)、姚琪等(2018a)模擬了這些具有重要控制作用的邊界斷裂的破裂演化特征，模擬結果顯示了在現今的應力場作用下，鮮水河—安寧河—小江斷裂系和大涼山斷裂的反復破裂-愈合過程，以及龍門山斷裂帶在長時間的應力加載下發生大規模破裂的過程.模擬所得的川滇地區主要斷裂的破裂演化特征與該地區強震的時空分布基本一致，表明該軟件所計算模擬的斷層摩擦破裂行為能夠代表強震的孕震-發震-震后愈合過程.

圖5 變形某時刻等效應力變化率的分布快照(Xing and Mora,2006)Fig.5 A snapshot of the equivalent stress rate distribution at a certain deformation stage (Xing and Mora,2006)

4 結論與展望

地震預報是世界性難題，其重要價值不言而喻.要取得重大突破，基于物理的地震數值預報不可或缺.大地震被認為是斷層上的一種摩擦破裂失穩行為，地震孕育和發生過程是斷層系統相互作用動力學行為的一種表現，因此斷層破裂摩擦動力學行為的研究對地震數值預報具有重要意義.本文圍繞斷層破裂摩擦動力學行為的有限元模擬，(1)概述了斷層摩擦動力學有限元模擬相關的理論、軟件及其優缺點，總結并比較了其中不同的時間積分方案和摩擦接觸的計算方法，指出傳統方案所面臨的摩擦非線性的收斂性難題，即使采用迭代平滑技術來減少劇烈變化滿足了收斂條件，也會導致黏著-滑移失穩(即地震發生的時刻)等關鍵劇變時刻漏失而不能被捕捉到；(2)提出并采用R-Minimum自適應地控制時間步長方案，將不穩定的隱式迭代計算轉變為顯式計算以保證各種非線性情況下的計算穩定性、精度和效率；同時通過節點-點接觸單元，結合庫侖摩擦準則利用速率和狀態相關的摩擦系數分析斷層非線性摩擦失穩行為；通過對相關應用分析總結，所提方案可自適應地控制時間步長以穩定高效地模擬地震過程，取得了較好的效果；(3)總結了基于斷層系統有限元網格生成的不同方案，討論了主流方法在復雜斷層系統四邊形及六面體網格生成中的弊端；通過復雜圖像實例，介紹了所研發的基于圖像的斷層系統四邊形和六面體網格自動生成步驟以及生成的網格模型在斷層系統非線性有限元模擬計算中的應用.

為了提高斷層系統有限元模擬分析平臺的有效性和實用性，還需在以下幾個方面進一步研究完善：

(1)多場耦合摩擦模型及其有限元實現：斷層摩擦非線性行為十分復雜，除了受到接觸面巖石物性、運動狀態、滑動速度和位移、接觸力的影響，還受到溫度、流體孔隙壓力和物理化學反應等各種耦合因素的影響.結合實驗室試驗和現場觀測，發展斷層摩擦多場耦合非線性本構方程及其穩定高效的有限元計算動力學模型與典型實驗驗證對于進一步探明多場耦合條件下摩擦失穩機理十分必要.

(2)多場耦合地震動力學綜合研究：斷裂帶尤其是俯沖帶是地球上構造活動最復雜、最強烈的區域.大洋巖石圈通過在匯聚板塊邊界的俯沖將大量水帶入到地幔中，并對俯沖帶地震的發生起到了重要的控制作用.隨著俯沖板片深度的增加，在一定的溫壓條件下，水化地幔(即蛇紋巖)發生脫水相變，引發俯沖帶中源地震.脫出的水則由于運移的差異，既可以產生板內的“水壓致裂”，也會影響俯沖界面的耦合行為，進而導致慢滑移地震區的形成.由此可見，俯沖帶地區深海-巖石圈流體交換及其在深部的效應是一個包含化學反應-溫度-流體流動-應力變形/破壞的多物理場耦合的復雜動力學系統.俯沖帶等中深部地震，雖然地震機理復雜，但多場耦合效應十分明顯.需要從地球系統科學的角度出發，將流體運移、化學反應與傳統的固體地球地震動力學研究相結合，著眼于多場耦合動力學綜合研究，對俯沖帶地區深海-巖石圈流體交換及其效應進行多時空尺度定量化表征和分析.

(3)典型區域計算網格模型、邊界條件、材料特性等參數確定與快速創建.我們知道，地震研究涉及眾多不同的數據.一般的斷層數據是由點數據集約束而成，如何由此生成三維地質體有限元網格并保證其中斷層系的幾何形狀、位置、交叉和整體網格質量是地震數值預測有限元模擬的先決條件，同時地震斷裂模擬過程中或許還需要網格重劃分來應對局部大變形和破壞所引起的網格畸變問題.因此計算網格自適應快速生成與重劃分是地震數值模擬預測所面臨的關鍵問題之一.另外，還需要綜合利用地震觀測、地質調查、GNSS/InSAR 形變測量、應力場觀測數據及歷史地震資料，對研究區諸多數據進行整合，協同并有效地使用研究區內不同方法得到的跨時空的多物理量觀測數據，以快速而準確地確定有限元計算模型的邊界條件、材料特性等關鍵參數，為后續的研究區斷層系統的有限元模擬分析和未來數值預報奠定基礎.

致謝感謝編輯與三位匿名評審的細致審閱、批改及建議，使得文章得到很大完善與提高．第一作者在讀研期間，經導師王仲仁教授介紹認識了王仁先生，并有幸參加了王仁先生主持的材料本構關系國際會議預備會并做匯報.謹以此文紀念兩位先生．