有限元數(shù)值模擬法基本原理及其在地質(zhì)構(gòu)造變形研究中的應(yīng)用綜述

（重慶科技學(xué)院石油與天然氣工程學(xué)院，重慶 401331）

傳統(tǒng)的構(gòu)造解析法是推斷地質(zhì)構(gòu)造變形演化過程的重要手段，但是由于構(gòu)造演化的動力學(xué)機制可以從不同的側(cè)重點進行分析，且對地質(zhì)過程的理解也不盡相同，所以傳統(tǒng)地質(zhì)分析方法在研究構(gòu)造變形機制中存在局限性。物理模型方法雖然能夠幫助人們理解構(gòu)造變形的演化過程，但此方法受到時空限制，不能有效的模擬出地質(zhì)構(gòu)造形態(tài)的復(fù)雜性與巖石物理性質(zhì)的多樣性[1]。隨著高性能計算機技術(shù)和圖形處理的快速發(fā)展，數(shù)值模擬的計算理論和計算方法已經(jīng)越來越成熟，在大陸構(gòu)造變形研究方面，從20世紀(jì)60年代的定性研究到后來的半定量甚至定量化研究，數(shù)值模擬成為了構(gòu)造變形研究中新的定量化方法，其應(yīng)用已涉及巖石圈流變學(xué)的模擬，褶皺、斷層的模擬，造山帶、俯沖帶與盆地的模擬等內(nèi)容。數(shù)值模擬方法可以綜合利用地質(zhì)、地球物理、地球化學(xué)等方法建立反映時空變化的各種地質(zhì)模型[2]，尤其適合大陸構(gòu)造變形和造山帶動力學(xué)等在空間和時間上跨度很大的領(lǐng)域的研究[3]。

數(shù)值模擬是一種計算機模擬分析方法，包括有限元法、有限差分法、邊界元法、離散元法等，其中有限元法和有限差分法廣泛應(yīng)用于構(gòu)造變形數(shù)值模擬，兩者在特定情況下得到的最終方程是一致的，但有限差分法適用于求解物理不穩(wěn)定性問題、非線性問題和大變形問題，而有限元法處理時間、步驟較少，應(yīng)對模擬線性、變形問題效率較高。因此，本文概述了有限元法數(shù)值模擬的基本原理和計算過程，列舉了構(gòu)造變形特征數(shù)值模擬的相關(guān)應(yīng)用，說明了有限元數(shù)值模擬法對構(gòu)造變形研究的優(yōu)勢，以及此方法在我國大陸構(gòu)造變形中研究進展。

1 有限元模擬的基本原理

1.1 基本原理

有限元法（FEM）是在三大守恒定律（質(zhì)量、動量、能量）的基礎(chǔ)上建立平衡方程、本構(gòu)方程和幾何方程，將連續(xù)的求解域分解為有限個單元的組合體，單元劃分越細，越接近于實際的地質(zhì)情況，計算結(jié)果就越精確。有限個單元的組合體中，單元節(jié)點量可將單元內(nèi)部結(jié)點的待求量輸入到有限元程序（如ANSYS、MARC、ABAQUS、ADINA等）中，采用最小勢能或虛位移原理建立結(jié)點平衡方程組，通過求解高階代數(shù)方程組，計算每個單元內(nèi)假設(shè)的近似值來表示工區(qū)內(nèi)所有的未知場函數(shù)。

有限元劃分出的所有單元必須滿足一個假設(shè)，有限元上的未知方程uci能被通過組合的局部節(jié)點uji近似逼近，其中uci是一組基于其節(jié)點的近似方程，并且uji必須滿足控制方程。未知方程uei和近似方程uji關(guān)系的公式如下：

式中：M—有限元的次序；—插值方程。如果上述公式成立，則以下的代數(shù)方程可以近似代替初始的控制方程：

其中系數(shù)矩陣的定義式為：

1.2 模型的建立

有限元法模型分析過程大致分為：地質(zhì)模型、力學(xué)模型、計算模型三個模型階段。

1）地質(zhì)模型

對于構(gòu)造變形演化的研究，建立地質(zhì)模型是必不可少的一步；對研究區(qū)域進行形態(tài)化，重點研究該地區(qū)地質(zhì)構(gòu)造形態(tài)、構(gòu)造組合持征；將地質(zhì)、測井、地球物理資料和各種解釋成果進行綜合分析，利用計算圖形技術(shù)繪制出三維定量特征值圖形，得到相關(guān)描述及定量數(shù)據(jù)，建立概化的地質(zhì)模型[1]。利用所建立的地質(zhì)模型和數(shù)值模擬的結(jié)果；作為分析研究地質(zhì)構(gòu)造變形的演化歷史、地質(zhì)構(gòu)造因素以及構(gòu)造成因的依據(jù)。

2）力學(xué)模型

力學(xué)模型的建立影響著數(shù)值模擬的過程及結(jié)果，所以它是構(gòu)造變形數(shù)值模擬的關(guān)鍵。在建立力學(xué)模型的時候是將地質(zhì)模型抽象化，并要考慮平面或者空間模型的選擇、材料本構(gòu)關(guān)系及參數(shù)大小的選擇、地層的劃分等[1]。

材料本構(gòu)關(guān)系及參數(shù)大小選擇要結(jié)合其他地質(zhì)研究方法，涉及巖石圈流變學(xué)的縱向和橫向分塊特性[1]，許多構(gòu)造變形數(shù)值模擬對巖石性質(zhì)、巖石圈流變學(xué)剖面參數(shù)是基于不同的構(gòu)造模式設(shè)定的。對大陸變形數(shù)值模擬中地質(zhì)材料及參數(shù)確定，以模擬構(gòu)造應(yīng)力場的需要為前提；地殼的介質(zhì)參數(shù)對構(gòu)造應(yīng)力場和位移場的變化也有至關(guān)重要的影響[4]。所以，在進行數(shù)值模擬研究構(gòu)造變形中，要選擇合乎大陸巖石層實際狀況的力學(xué)參數(shù)和經(jīng)過多種多次試驗撿驗，形成合理經(jīng)擬合的參數(shù)組合[4]。

地層應(yīng)根據(jù)可表征地質(zhì)體特點進行劃分，對所建立的地質(zhì)模型進行優(yōu)化并簡化。地質(zhì)構(gòu)造變形數(shù)值模擬涉及大變形非線性和多場耦合問題，具有大量的計算量，優(yōu)化并簡化地層，減少模型層數(shù)[1]，也可在模型邊界和斷層作適當(dāng)?shù)暮喕推交幚韀4]。

不同尺度的地質(zhì)體和不同的位置有著不同地質(zhì)演化過程，在進行力學(xué)模型建立時應(yīng)該根據(jù)具體的研究對象進行相對應(yīng)的分析。

3）計算模型

計算模型將抽象的力學(xué)模型轉(zhuǎn)化為具體方程并進行求解的過程，涉及數(shù)學(xué)模型，即方程建立和計算軟件的選擇與使用。

在軟件方面，主要是使用常用的數(shù)值模擬軟件進行模型數(shù)值的計算，通常采用算法較為可靠，計算過程和計算結(jié)果處理相對簡單的軟件，若通用軟件無法直接解決現(xiàn)有模型的計算，可自主進行算法設(shè)計并研制編寫程序，進行相應(yīng)的數(shù)值分析。由于構(gòu)造變形數(shù)值模擬大多是非線性問題，所以要選擇恰當(dāng)?shù)乃惴╗1]。

構(gòu)造變形等于連續(xù)介質(zhì)變形，連續(xù)介質(zhì)變形的能量來源于作用在該介質(zhì)上內(nèi)力和外力的平衡[3]，而這種力與變形有關(guān)的關(guān)系方程式的描述就需要用到動量方程，描述連續(xù)介質(zhì)在重力場作用下的動量守恒，作用在物體上的合外力的大小等于物體在力作用方向上的動量的變化率。計算模型建立需要利用數(shù)學(xué)動量方程作為力學(xué)模型的轉(zhuǎn)換。

1.3 邊界條件的確定

為了對方程解進行約束，有限元法還需定義初始條件和邊界條件。邊界條件是是控制大尺度構(gòu)造變形的主要因素，需要板塊理論、區(qū)域地質(zhì)研究和GPS測量結(jié)果來支撐。掌握巖體計算力學(xué)方法，深入?yún)^(qū)域地質(zhì)研究，對確定模型邊界條件最重要。構(gòu)造變形模型的邊界條件由邊界處的區(qū)域化應(yīng)力、位移變量值組成。力學(xué)邊界有應(yīng)力邊界（即給定物體表面上的面力或集中力）、位移邊界（即邊界上各點的位移分量）和混合邊界（即在一方面為已知的邊界力而在另一方面為已知的位移），應(yīng)力邊界條件和位移邊界條件尤為重要。

1.3.1 應(yīng)力邊界條件

彈性體內(nèi)部各點的約束條件在其邊界上的延續(xù)就是應(yīng)力邊界條件，主要確定方法是公式法。在空間問題中，應(yīng)力邊界條件的求解公式為：

式中：l、m、n——外法線方向分別與x、y、z軸正向夾角的余弦值，x、fyf、zf——面力分別在x、y、z軸上的投影。

平面問題中，公式可簡化為：

此外，比較法也是一種有效的方法，在求解過程中，先進行研究對象的力學(xué)分析，當(dāng)選用適當(dāng)?shù)淖鴺?biāo)系時，邊界的坐標(biāo)軸與外法線方向通常平行或垂直，這時只需將微元體“移到”邊界附近，并比較微元體的應(yīng)力分量和邊界面的面力分量fx、fy即可確定應(yīng)力邊界條件。

兩種方法各有優(yōu)點，比較法是當(dāng)邊界面與某一坐標(biāo)軸垂直或平行時的特例，求解應(yīng)力邊界條件既簡便正確性又高。當(dāng)邊界面既不平行也不垂直于坐標(biāo)軸時，只能應(yīng)用公式法，所以說公式法適用的范圍更廣泛。但是，處理同一個問題時，兩種方法得出的結(jié)果是完全相同的。

1.3.2 位移邊界條件

構(gòu)造應(yīng)力場的位移邊界條件具有一般性，為方便描述，以圖中模型的邊界為例，假設(shè)平行于東西方向的y軸兩個側(cè)面沿著y軸正方向依次為邊界1和邊界2；平行于南北方向的x軸兩個側(cè)面沿著x軸正方向依次為邊界3和邊界4[5]。

模擬地質(zhì)構(gòu)造力作用的基本因素包括以下4種[6]：自重應(yīng)力；東西向和南北向水平向均勻擠壓構(gòu)造運動；東西向和南北向垂直平面內(nèi)豎向均勻剪切變形構(gòu)造運動；水平面內(nèi)的均勻剪切變形構(gòu)造運動。

簡化計算模型[5]

不同基本因素下所對應(yīng)的位移邊界條件：

邊界1+邊界2：

邊界3+邊界4：

邊界1+邊界2+邊界3+邊界4：

式中：n為多項式的最高次數(shù)；i為對應(yīng)于位移u、v、w的第i次多項式，ai、bi、ci為位移模式系數(shù)，因為y在邊界1和邊界2上是常數(shù),所以此邊界條件不含有y，其他邊界依次類推。

3 構(gòu)造變形中有限元法的應(yīng)用

近年來，有限元法等數(shù)值模擬方法在構(gòu)造變形方面，已應(yīng)用到造山帶、俯沖帶、盆地、褶皺以及斷層的模擬等方面。

3.1 盆地的相關(guān)模擬

盆地主要是由于地殼運動形成的。通過地殼運動作用，地下的巖層受到拉伸、擠壓形成彎曲或斷裂，使部分巖石隆起，部分下降，若下降部分被隆起部分所包圍，就形成了盆地的雛形。在構(gòu)造類型中，以斷裂為邊界的張性盆地和壓性盆地發(fā)育較為廣泛，其中伸展盆地和前陸盆地與石油、天然氣、礦產(chǎn)資源和地震活動有著十分密切的關(guān)系。因此，對盆地構(gòu)造特征的模擬一直是廣大學(xué)者們研究的重點。

Bott（1997）[7]利用赫彈塑性有限元方法對上地殼中由高角度平面正斷層界定的半地塹演化機制進行了模擬，計算了演化各個階段的應(yīng)力分布、塑性破裂分布和撓曲剖面，并發(fā)現(xiàn)了塑性破裂和有效彈性厚度之間的相互關(guān)系。梁海華等（1997）[8]根據(jù)地質(zhì)模型用有限元數(shù)值計算方法模擬了古構(gòu)造應(yīng)力場和構(gòu)造變形場，發(fā)現(xiàn)擠壓作用對盆地各段邊界上的方向有細微的變化，盆地的不對稱形態(tài)以及主應(yīng)力方向隨深度的變化而變化。譚成軒等（1997）[9]運用SuperSAP有限元程序進行三維構(gòu)造應(yīng)力場的數(shù)值模擬，探索了含油氣盆地三維構(gòu)造應(yīng)力場數(shù)值模擬方法，定量化研究了排烴史和油氣運移聚集史，對油田合理注采開發(fā)，尋找殘余油等均具有重要意義。王喜雙等（1997）[10]利用有限元方法對盆地進行了應(yīng)力場模擬，找出了現(xiàn)代應(yīng)力場與油氣聚集的關(guān)系，并以比較地質(zhì)學(xué)思想為基礎(chǔ)，推出了地質(zhì)歷史演化中各時期應(yīng)力場與油氣聚集的相互關(guān)系，得出了盆地在構(gòu)造應(yīng)力場控制下的油氣聚集特征。王喜雙等（1999）[11]又建立了盆地結(jié)構(gòu)和地層資料三維地質(zhì)模型，用三維有限單元方法分析出在構(gòu)造力和重力共同作用下，盆地從淺至深各構(gòu)造層應(yīng)力變化的特征。Zhang G.B.等（2000）[12]模擬了上地殼中由高角度平面逆斷層界定的非對稱盆地演化機制，并對比了正斷層界定的半地塹演化機制，模擬了拉張半地塹和非對稱壓性盆地兩者的撓曲剖面在不同時期的演化過程。張貴賓等（2000）[13]又利用粘彈塑構(gòu)造模擬軟件包對盆地的動力學(xué)演化進行模擬，得出塑性變形使得有效彈性厚度減小，盆地寬度和極限深度受到上地殼層厚度和沉積物密度兩個因素所控制。羅曉容（2000）[14]應(yīng)用數(shù)值模型方法研究異常流體壓力機制及演化過程，得出地溫梯度的增加使得地層壓力變小，地層滲透率越低，有機質(zhì)裂解生成氣態(tài)烴的增壓效應(yīng)不一定越顯著。李鐵飛等（2013）[15]以沉積環(huán)境為基礎(chǔ)，通過有限元法建立了盆地充填過程的數(shù)值模型，研究了沉積環(huán)境對沉積物物理參數(shù)的影響。

3.2 造山帶的相關(guān)模擬

造山帶是由巖石圈劇烈構(gòu)造運動所造成的狹長構(gòu)造變形帶，往往在地表形成線狀相對隆起的山脈。造山帶是巖石圈中最活躍的地帶，也是巖石圈與軟流圈、地殼與地幔之間大規(guī)模物質(zhì)與能量的交換地帶。因此，造山帶是了解地球深部與淺表地殼構(gòu)造關(guān)系最好的場所，對造山帶的研究是現(xiàn)今構(gòu)造變形中數(shù)值模擬研究的核心問題之一。

Willett等（1994）[16]利用有限元方法進行模擬，表明巖石圈地幔并非與地殼一起變形，并認(rèn)為繼印度板塊向北碰撞以后，亞洲巖石圈地幔開始向南的俯沖并伴以向北的后退。Willett（1999）[17]又利用有限元方法對在擠壓造山環(huán)境下出現(xiàn)拉張?zhí)卣髋c不同流變學(xué)性質(zhì)之間的關(guān)系做了系統(tǒng)的研究。曾佐勛等（2001）[18]通過對非連續(xù)介質(zhì)大變形有限元數(shù)值模擬，為喜馬拉雅造山帶楔狀擠出模式、東阿爾卑斯三維管狀擠出模式、西秦嶺和東松潘-甘孜復(fù)合造山體的三維滑脫擠出模式等提供了力學(xué)依據(jù)。何建坤等（2002）[19]應(yīng)用不可壓縮非牛頓粘性流體的本構(gòu)關(guān)系和二維有限元法，在造山帶同擠壓期對下地殼流變與上地殼構(gòu)造伸展作用的動力學(xué)關(guān)系進行了研究，探討了造山帶下地殼在特定邊界條件下的黏性流變，以及對上地殼動力學(xué)演化的制約，結(jié)果表明，在側(cè)向擠壓作用下，地殼不同圈層巖石的粘性流變的分布范圍和運動學(xué)方式既受時間的影響，還受地殼厚度轉(zhuǎn)變帶形態(tài)的制約。Mclellan等（2006）[20]基于裂谷和沉陷模式進行模擬，推出變形、流體流動、熱平流全耦合的盆地過程及它們對流體流動的影響，并且與由結(jié)構(gòu)造成的流體流動效果進行了比較。王洪亮等（2001）[21]基于非牛頓流體近似的有效黏度模型數(shù)值模擬了巖石圈拆沉的過程，分析了巖石圈的黏度結(jié)構(gòu)對巖石圈拆沉作用的影響，得出下地殼對地殼與巖石圈地幔的耦合程度起控制作用，且對拆沉作用的過程和形態(tài)有很大的影響。胡秋媛等（2015）[22]運用基于有限元方法進行了三維構(gòu)造應(yīng)力場數(shù)值模擬，得出了伸展構(gòu)造演化對坳陷區(qū)油氣的聚集及隆起區(qū)金屬礦產(chǎn)的富集具有重要的控制作用。

3.3 俯沖帶的相關(guān)模擬

俯沖帶指在大洋板塊俯沖于大陸板塊之下的構(gòu)造帶，是發(fā)生在板塊匯聚邊緣主要部分，作為地球物質(zhì)發(fā)生循環(huán)的重要場所，俯沖帶的相關(guān)研究一直以來都受到國內(nèi)外地質(zhì)學(xué)家的高度重視（張繼等，2015）[23]。

在不同運動學(xué)參數(shù)和熱參數(shù)的條件下，臧紹先等（1994）[24]使用有限元法和準(zhǔn)動力學(xué)模型計算并得出了俯沖帶的溫度、密度分布，結(jié)合負(fù)浮力對板塊形成等效應(yīng)力的影響，發(fā)現(xiàn)在板塊俯沖的過程中負(fù)浮力和等效應(yīng)力均是變化的。

石耀霖等（1998）[25]通過使用有限元法模擬計算了活動海嶺俯沖的熱演化過程，分析了周圍地區(qū)地形的變化，還合理解釋了在板塊俯沖過程中所產(chǎn)生的火山活動間斷現(xiàn)象。Devaux等（2000）[26]利用橄欖石—尖晶石相變動力學(xué)的計算結(jié)果，通過使用二維彈性和粘彈性有限元法，模擬了亞穩(wěn)態(tài)橄欖石楔存在時俯沖板塊深部的各種應(yīng)力，發(fā)現(xiàn)由相變產(chǎn)生的應(yīng)力遠遠大于由密度引起的浮力所產(chǎn)生的應(yīng)力熱應(yīng)力。劉亞靜等（2002）[27]基于俯沖帶流變結(jié)構(gòu)模型基礎(chǔ)，運用粘彈性平面有限元法，主要對俯沖帶的構(gòu)造應(yīng)力場特征進行數(shù)值模擬，認(rèn)為若存在橄欖石—尖晶石相變過渡區(qū)的粘度比上部橄欖石及下部尖晶石區(qū)域粘度更低的情況，則此粘度結(jié)構(gòu)只使相變過渡區(qū)附近的最大剪應(yīng)力降低，而對應(yīng)力分布的整體圖像不會產(chǎn)生影響。姜輝（2012）[28]結(jié)合地震資料和物理觀測數(shù)據(jù)，建立了關(guān)于俯沖帶在巖石圈有無彎折的四個模型，使用有限元法模擬了俯沖帶地震粘滑失穩(wěn)過程，重點研究分析了俯沖帶幾何形態(tài)變化對其粘滑事件發(fā)生的影響。

3.4 斷層的相關(guān)模擬

斷層是一種在地殼中廣泛發(fā)育的地質(zhì)現(xiàn)象，是指當(dāng)巖體受到構(gòu)造應(yīng)力作用后形成不連續(xù)的破裂面，并沿此破裂面發(fā)生了明顯位移的構(gòu)造。關(guān)于斷層性質(zhì)，前人利用有限元法做了大量研究。

唐湘蓉等（2005）[29]從裂縫形成的古構(gòu)造應(yīng)力場出發(fā)，利用連續(xù)介質(zhì)有限元數(shù)值模擬法和巖石破裂準(zhǔn)則，對某古潛山儲層裂縫的方位、密度分布進行研究，成功地預(yù)測了該研究區(qū)內(nèi)裂逢發(fā)育帶的宏觀平面分布。沈海超等（2007）[30]使用有限元法對斷層區(qū)塊進行了模擬，研究了影響斷層附近地應(yīng)力場的敏感因素，由于斷層及其附近應(yīng)力的變化比較復(fù)雜，高應(yīng)力集中區(qū)主要在斷層端部和幾何形態(tài)的拐點處，通過模型的建立和分析，結(jié)果表示這些敏感因素對斷層附近的應(yīng)力場都有明顯的影響，只是各個因素的影響程度不同。程遠方等（2008）[31]采用有限元約束優(yōu)化反演法對斷塊的三維應(yīng)力場進行了模擬，主要探討了斷層對地應(yīng)力場的影響，研究結(jié)果表明了鄰近斷層、擾亂程度、斷層的幾何尺寸和距離斷層的遠近等是斷層附近的地應(yīng)力分布主要影響因素。孫禮健等（2009）[32]用有限元方法模擬研究斷層構(gòu)造對地應(yīng)力場分布以及斷層端部的影響，發(fā)現(xiàn)了斷層規(guī)模、幾何形態(tài)、彈性模量、巖體、斷層走向與區(qū)域最大水平主應(yīng)力的夾角、邊界應(yīng)力比等因素對斷層端部以及附近應(yīng)力場分布有著不同程度的影響。白玉柱等（2013）[33]以對逆斷層發(fā)生逆沖構(gòu)造運動時造成的變形以及相應(yīng)的應(yīng)力場分布建立了有限元數(shù)值模型，模擬得出了由斷層錯動形成的地表的變形以及相應(yīng)的應(yīng)力分量場的變化情況。袁杰等（2014）[34]為了研究斷層自發(fā)破裂的動力過程，利用有限元法對此過程進行了動態(tài)數(shù)值模擬，發(fā)現(xiàn)經(jīng)典的滑移弱化摩擦關(guān)系并不能產(chǎn)生“脈沖型”破裂模式，但如果用改進過后的摩擦關(guān)系則能產(chǎn)生此模式。

3.5 褶皺的相關(guān)模擬

褶皺是指在構(gòu)造運動的作用下，巖層受力而發(fā)生彎曲，并發(fā)生一系列波狀的彎曲變形。作為地質(zhì)界重要的地質(zhì)現(xiàn)象之一，褶皺的形成原理和產(chǎn)生機制是科學(xué)家們一直以來所重點研究的內(nèi)容。在關(guān)于褶皺的眾多研究中，為人所知的觀點主要有兩種。第一種觀點是主波長理論，即在受力作用后巖層力學(xué)性質(zhì)發(fā)生變化，并控制最終褶皺的波長，第二種觀點是褶皺最終的形態(tài)特點決定于其巖層的最初始幾何特征（龔紀(jì)文等，2002）[35]。

Zhang Y等（1996）[36]從巖層的力學(xué)性質(zhì)及初始幾何形態(tài)對褶皺進行分析研究，模擬巖層受力變形后形態(tài)的變化，以及褶皺最后的波長，并補充了對主波長的見解，認(rèn)為地層的褶皺在“無限小”初始擾動下，褶皺的最終決定形態(tài)及波長不僅取決于巖層厚度，還由巖層的力學(xué)性質(zhì)所共同決定，在褶皺“無限小”初始擾動下，褶皺的最終形態(tài)則會受其影響。Zhang Y等（2000）[37]在利用有限元法建立褶皺數(shù)值模型后，發(fā)現(xiàn)對結(jié)果產(chǎn)生較大影響的因素為應(yīng)變速率，比奧理論認(rèn)為是層間黏性差異產(chǎn)生的不穩(wěn)定現(xiàn)象即為褶皺，但該傳統(tǒng)理論方法不適用于黏性差異極小的情況。楊玉山等（2006）[38]建立了關(guān)于褶皺變形應(yīng)力場的模型，使用有限元法模擬了變形應(yīng)力場對同生斷裂的控制作用，模擬結(jié)果發(fā)現(xiàn)裂隙的延伸方向與剪應(yīng)力大小不存在直接的聯(lián)系，主要受剪切應(yīng)力方向控制。Sanz等（2008）[39]通過有限元法研究了層面變形過程中的滑動，選取具有代表性的3層地層作為研究，對不同巖層厚度進行處理，將不同巖層力學(xué)性質(zhì)組合，模擬了40種條件下的二維平面應(yīng)變，研究了巖層各種性質(zhì)對褶皺形成的影響。為了研究煤層褶皺構(gòu)造形成演化的應(yīng)力場分布，付京斌（2010）[40]利用有限元法建立了模型，模擬了在構(gòu)造作用下煤層褶皺的形成過程，模擬結(jié)果發(fā)現(xiàn)含向斜結(jié)構(gòu)煤層兩翼的應(yīng)力低于軸部應(yīng)力，含背斜煤層兩翼應(yīng)力高于軸部應(yīng)力。

4 有限元法在我國大陸構(gòu)造變形中的應(yīng)用

有限元法廣泛應(yīng)用于大尺度構(gòu)造變形數(shù)值模擬中，且多集中在構(gòu)造演化過程。由于青藏高原在構(gòu)造演化中的獨特性，以及有限元法在研究構(gòu)造變形演化機制方面的優(yōu)越性，使得有限元法廣泛應(yīng)用于青藏高原的研究中，同時青藏高原也是世界上應(yīng)用有限元法最多的區(qū)域之一。在現(xiàn)今青藏高原運動學(xué)、動力學(xué)和構(gòu)造變形機制方面，國內(nèi)專家學(xué)者應(yīng)用有限元法取得了一定的研究成果。鄭勇等（2007）[41]運用有限元法，在GPS測量結(jié)果作為邊界條件下探討了斷層的存在與否對青藏高原現(xiàn)代運動場分布的影響，研究表明斷層的滑移運動增強了青藏高原東西兩側(cè)的拉張趨勢，加大了青藏高原始物質(zhì)東移的速度，改變了柴達木盆地和塔里木的運動狀態(tài)，充分說明了在研究青藏高原動力學(xué)機制中，斷層作用是不可忽視的因素。張東寧等（1994）[42]建立了青藏高原三維彈－粘性有限單元模型，模擬了由地震震源機制解和地質(zhì)調(diào)查資料推斷的高原巖石圈的構(gòu)造應(yīng)力狀態(tài)和巖石圈物質(zhì)的運動特征。得出了青藏高原物質(zhì)目前存在被向東擠出的水平運動，可能導(dǎo)致青藏高原北部及東部地區(qū)的現(xiàn)代左旋走滑斷層。鄧宗策等（1990）[43]建立了青藏高原總體構(gòu)造走向的南北向直立剖面的有限元模型，對成分層和有限單元進行了劃分。對此模型進行彈性材料的計算模擬和分析，得出了青藏高原受印度板塊向北運動擠壓、高原北部巖石圈阻礙、軟流圈拖曳、重力及其均衡調(diào)整作用、地殼和上地幔結(jié)構(gòu)構(gòu)造的影響，其現(xiàn)今的地殼和上地慢是經(jīng)過長期地質(zhì)歷史演化的結(jié)果。傅容珊等（2000）[44]用模擬了青藏高原的擠壓隆升演化過程，剝蝕修正了數(shù)值模擬的隆升過程，結(jié)果發(fā)現(xiàn)由擠壓模型所產(chǎn)生的地形和青藏高原及其鄰區(qū)的地形的吻合性。同時也發(fā)現(xiàn)了擠壓隆升過程受巖石層的力學(xué)特性、邊界條件以及剝蝕作用等因素的制約，高原隆升都是不均勻的演化過程。李巖峰等（2004）[45]利用黏彈性有限元法，模擬研究了青藏高原西、東剖面的地殼均衡和巖石圈根拖曳的構(gòu)造應(yīng)力機制，結(jié)果表明青藏高原西部巖石圈根的向下拖曳是造山水平擠壓力的主要來源，印度板塊向北擠壓為次要因素；而青藏高原東部巖石圈根向下拖曳則不是硬上地殼中擠壓造山的主要力源，反映了高原內(nèi)部造山演化的西、東分異特征。

此外，國內(nèi)專家學(xué)者在川滇地區(qū)、大巴山弧形構(gòu)造帶、華北盆地、喜馬拉雅造山帶等區(qū)域的構(gòu)造變形數(shù)值模擬中也有相應(yīng)的研究成果。

楊光宇（1985）[46]采用平面問題的有限元法，參考了近十年四川及云南跨斷層短基線的測量資料以及查閱43個地震機制的結(jié)果，對川滇地區(qū)強震活動與構(gòu)造應(yīng)力場的關(guān)系作了初步研究，并推演出川滇地區(qū)可能發(fā)生地震的危險地段。陳連旺等（2008）[47]建立了川滇地區(qū)三維非線性有限元模型并進行分析，研究表明：龍門山斷裂帶應(yīng)力狀態(tài)在汶川地震發(fā)生后的變化特征；汶川地震余震的序列空間分布特征及其對震源機制類型的力學(xué)解釋；汶川地震引發(fā)的川滇地區(qū)加載效應(yīng)。陳祖安等（2009）[48]用三維流變非連續(xù)有限元法，在青藏高原、龍門山斷裂帶及鄰區(qū)三維構(gòu)造塊體的背景下結(jié)合GPS 資料，計算出研究區(qū)邊界斷層上表征剪應(yīng)力及法向應(yīng)力，其結(jié)果顯示了汶川大地震的孕育發(fā)生機理，并對川滇地區(qū)的地震危險性提供了有力的參考依據(jù)。

華北盆地多年來陸續(xù)發(fā)生過多次破壞性大地震，柳暢等（2012）[49]建立了華北盆地巖石圈三維黏彈性有限元模型，分析了地殼分層流變性質(zhì)和地殼結(jié)構(gòu)對地殼應(yīng)力積累的影響。結(jié)果表明，華北盆地莫霍面的隆起與黏滯系數(shù)是本區(qū)域地震孕育的重要因素，莫霍面隆起處因構(gòu)造擠壓的持續(xù)作用產(chǎn)生了明顯應(yīng)力集中現(xiàn)象，揭示了華北地殼從脆性的巖石圈上地殼逐漸變成韌性的巖石圈上地幔的分層流變結(jié)構(gòu)。

大巴山弧形構(gòu)造位于秦嶺造山帶的南緣，它的形成機制歷來是研究的熱點。在前人的應(yīng)力場反演工作成果上，王瑞瑞等（2013）[50]運用數(shù)值模擬的方法，模擬出大巴山主造山期的位移場和應(yīng)力場，結(jié)果表明大巴山弧形構(gòu)造由三個主要因素控制而形成：早期伸展背景下形成的弧形斷裂邊界；前緣兩側(cè)隆起的砥柱作用；底部滑脫作用。武紅嶺等（2009）[51]采用非線性有限元法，通過大巴山前陸構(gòu)造帶應(yīng)力場、位移場和形變場流變學(xué)特征，研究了疊加褶皺的復(fù)合形式和聯(lián)合弧的形態(tài)特征。結(jié)果表明，大巴山弧形構(gòu)造帶從西到東的早晚兩期褶皺大致有構(gòu)造橫跨疊加、斜跨疊加和共軸疊加3種基本相干形式；北東向區(qū)域應(yīng)力的驅(qū)動和東西兩個基底構(gòu)造結(jié)的存在是大巴山弧形構(gòu)造帶形成的主要因素。

在喜馬拉雅造山帶研究中，曾佐勛等（2001）[18]采用非連續(xù)介質(zhì)大變形有限元法，結(jié)合前人工作的成果，模擬了高喜馬拉雅結(jié)晶地體的垂向擠出、東阿爾卑斯三維管狀擠出和陜甘川鄰接區(qū)復(fù)合造山體的三維滑脫擠出，提出了東阿爾卑斯三維管狀擠出模式。

5 存在問題及發(fā)展趨勢

近年來，有限元法等數(shù)值模擬方法在地質(zhì)學(xué)上應(yīng)用越來越廣。但有限元法在地質(zhì)上所模擬的構(gòu)造變形仍存在一些問題：要把所研究的地質(zhì)對象進行離散化，涉及到復(fù)雜的計算和地質(zhì)分析，對地質(zhì)工作研究來說，是一個較大的難題；要進行大尺度構(gòu)造變形的模擬，常遇到非線性問題，故很容易發(fā)生收斂困難而使計算中斷；任何構(gòu)造變形的模擬試驗，其條件均不可能與自然界完全一致，會存在一定的干擾因素，因而模擬結(jié)果常出現(xiàn)異常[1]。

隨著數(shù)值模擬理論的不斷發(fā)展和計算機技術(shù)的日益成熟，有限元法能夠更加廣泛的應(yīng)用于地震預(yù)測科學(xué)和構(gòu)造應(yīng)力場分析中去：有限元法有較強的適應(yīng)性，能夠適用于地震動力過程和各種構(gòu)造地質(zhì)過程的模擬，有利于研究者進行擴展研究[52]；在地質(zhì)資料的不斷充實的情況下，有限元法能夠更加真實的反映各種地質(zhì)構(gòu)造的變化，使模擬結(jié)果的精度有一定提升；隨著計算方法的不斷完善，需要花在計算上的時間會大大減少，為地質(zhì)研究人員節(jié)省了人力、物力。

有限元法是一種綜合分析方法，和其他使用方法一樣，未來的發(fā)展會使其變得更加高效、簡單，以后在遇到大型復(fù)雜的地質(zhì)問題時，有限元法能得到更好的應(yīng)用。