構(gòu)造應(yīng)力場模擬
——有限元理論、方法和研究進展①

張勝利

（1.五邑大學(xué)信息學(xué)院，廣東江門 529020；2.中科院廣州地化所，廣東廣州 510640）

構(gòu)造應(yīng)力場模擬
——有限元理論、方法和研究進展①

張勝利1，2

（1.五邑大學(xué)信息學(xué)院，廣東江門 529020；2.中科院廣州地化所，廣東廣州 510640）

采用有限元數(shù)值模擬方法對構(gòu)造地質(zhì)問題進行描述和定量化求解是當(dāng)前地質(zhì)學(xué)領(lǐng)域的研究的一個熱點，在近10年以來取得了重要進展，形成了比較完整的理論和技術(shù)體系，并在一些典型的地質(zhì)構(gòu)造帶獲得了重要的研究成果。本文以有限元數(shù)值模擬方法理論作為出發(fā)點，總結(jié)分析了國內(nèi)外有限元數(shù)值模擬方法在構(gòu)造應(yīng)力場領(lǐng)域的研究進展情況和技術(shù)方法，并討論了其目前存在的問題和未來發(fā)展方向。

構(gòu)造應(yīng)力場；數(shù)值模擬；有限單元法

Abstract：The Finite Element Method（FEM）has been used in the study of tectonic stress field for a long time，and the essence of numerical modeling has been adopted to the well－established numerical methods of multidisciplinary acknowledge including mathematics，physics and mechanics for studing characters of geological tectonics.In the last decade，great advances have been made on the numerical simulation method，not only an integrated theory has been built up，but also some significant results have been born from several typical tectonic belts.So the FEM becomes one of the most important numerical methods in the study of tectonic stress field.In this paper，taking theory of FEM as a springboard，the new progress and methods in this field at home and abroad is summarized and analyzed.Some problems and prospect of the researching on the field is also given.

Key words：Tectonic stress field；Numerical model；Finite element method

0 引言

地殼中的各種地質(zhì)構(gòu)造都是巖石受力發(fā)生變形的產(chǎn)物，它們的產(chǎn)生和發(fā)展必然也受力學(xué)規(guī)律的支配。因此研究構(gòu)造應(yīng)力場特征是闡明和解釋各種地質(zhì)構(gòu)造的產(chǎn)生、分布和演化規(guī)律的重要途徑。但由于地質(zhì)構(gòu)造演化過程的特殊性，采用傳統(tǒng)的地質(zhì)－地球物理－地球化學(xué)的方法研究應(yīng)力場特征顯得異常困難。而物理模擬雖然可以幫助人們理解構(gòu)造變形過程，但是存在著嚴重的時空尺度的局限性，并且一些因素在實驗室條件下尚不能獲得，尤其是在伴隨著高溫高壓的地球深部構(gòu)造帶。

有限單元法又稱有限元法，是上世紀50年代中期發(fā)展起來的一種數(shù)值計算方法，最初主要用于結(jié)構(gòu)和應(yīng)力分析，60年代后期被引入到地質(zhì)構(gòu)造分析中，由于其巨大優(yōu)越性而被迅速推廣，已成功的解決了很多地學(xué)的難題。我國是從70年代中期開始將有限元數(shù)值模擬技術(shù)應(yīng)用在地學(xué)領(lǐng)域的。隨著有限元數(shù)值模擬技術(shù)的出現(xiàn)，在綜合利用其它方法（如地質(zhì)、地球物理、地球化學(xué)）研究成果的基礎(chǔ)之上，采用這類數(shù)值模擬技術(shù)，可以建立不受時空限制的各種地質(zhì)模型，幫助研究人員比較直觀地認識和理解構(gòu)造應(yīng)力場的特征和演化歷史，并可以對其演化趨勢做出預(yù)測?，F(xiàn)在該技術(shù)除了在基礎(chǔ)地質(zhì)方面得到廣泛應(yīng)用外，在天文、環(huán)境、礦產(chǎn)勘探、地震預(yù)測等方面也取得了眾多的研究成果［1－3］。

隨著計算機計算速度的提高，構(gòu)造應(yīng)力場的數(shù)值模擬方法也獲得了很大的進展。國內(nèi)外的研究結(jié)果表明，目前已經(jīng)可以在構(gòu)造應(yīng)力場的模擬過程，在不同的本構(gòu)關(guān)系下，將流體、溫度、剝蝕、斷層等因素考慮進去，進行二維和三維的數(shù)值模擬［4－5］。但是也存在很多的問題，例如模擬的主觀性（邊界條件、模型等）、模擬結(jié)果的多解性和可靠性等。

1 基本理論

有限元法的基本原理是根據(jù)三大守恒定律（質(zhì)量、動量、能量）建立三個偏微分方程：平衡方程、本構(gòu)方程和幾何方程。然后將擬分析的連續(xù)體假想地分為有限個單元，離散后的單元與單元之間利用單元節(jié)點相互連接起來。單元內(nèi)部點的待求量可由單元節(jié)點量通過選定的插值函數(shù)求得。由于單元形狀簡單，易于用平衡關(guān)系或能量關(guān)系建立節(jié)點量之間的方程式。然后將各個單元方程“組裝”在一起形成總體代數(shù)方程組，計入邊界條件后即可對方程組求解。單元劃分越細，計算結(jié)果就越精確。

采用有限元法模擬構(gòu)造應(yīng)力場一般包括一下幾個步驟：

（1）在詳細分析研究目標地質(zhì)資料的基礎(chǔ)之上建立計算模型，這也是數(shù)值模擬試驗成敗的關(guān)鍵。

（2）單元剖分和選擇插值函數(shù)。根據(jù)模型的幾何特性、載荷情況及所要求的變形點，將模型離散化，單元大小應(yīng)當(dāng)可以給出有用的結(jié)果，但又要盡可能的節(jié)省計算量。

（3）確定應(yīng)變位移和應(yīng)力應(yīng)變關(guān)系。應(yīng)力與應(yīng)變之間關(guān)系由本構(gòu)方程確定，不同的本構(gòu)關(guān)系對應(yīng)不同的物理介質(zhì)。目前構(gòu)造應(yīng)力場模擬中較常用的本構(gòu)關(guān)系有4種：線彈性、彈塑性、粘彈性和彈粘塑性本構(gòu)關(guān)系［6］。

其中彈塑性本構(gòu)一般用于模擬地殼淺部的脆性層，粘彈性和彈粘塑性本構(gòu)多用于模擬地殼深部。也可以在同一個模型中采用多種本構(gòu)關(guān)系。

（1）線彈性本構(gòu)關(guān)系：

其中σij是應(yīng)力張量；eij是應(yīng)變張量；λ和μ是拉梅系數(shù)；θ是體積模量。線彈性本構(gòu)一般用于模擬小范圍的地殼淺部沉積層的應(yīng)力場。

（2）彈塑性本構(gòu)關(guān)系：

其中Dijkl為彈性系數(shù)；εkl為應(yīng)變張量；dλ是由一致性條件確定的尺度因子；F是用應(yīng)變空間表示的屈服函數(shù)。彈塑性本構(gòu)一般用于模擬地殼淺部脆性層。

（3）粘彈性本構(gòu)關(guān)系（以Kelvin體為例）：其中E為彈性模量；εij為應(yīng)變速率張量；η是粘性系數(shù)。式（3）也常被稱為線性流變方程。

（4）彈粘塑性本構(gòu)關(guān)系：

其中γ是控制塑性流動速率的流度參數(shù)，是時間和粘塑性應(yīng)變張量的不變量的函數(shù)；f是屈服函數(shù)；f0是一個使表達式無量綱化的、正的參考值。公式（2）、（3）也可以視作本式的特殊形式。

2 研究歷史及現(xiàn)狀

構(gòu)造應(yīng)力場模擬和板塊碰撞、大陸變形數(shù)值模擬緊密相關(guān)。Bott［7］對板緣力驅(qū)動下板內(nèi)應(yīng)力應(yīng)變的演化模擬可以說是這方面的開山之作。1976年Tapponnier［8］等人用理想剛塑性體所產(chǎn)生的變形來模擬亞洲大陸自新生代以來所經(jīng)歷的大尺度碰撞變形和斷層的產(chǎn)生，并為板塊內(nèi)部強烈變形提供了有力證據(jù)，從而使得應(yīng)力場的數(shù)值模擬方法引起了人們的廣泛關(guān)注。England［9］等在1982年首次將粘性薄層流變模型運用到以印度板塊和歐亞板塊的碰撞為代表的大陸變形帶，并考慮了地殼增厚造成的應(yīng)力變化，對青藏高原的擠壓隆升演化進行了數(shù)值模擬；England和Houseman等［10－11］也成功模擬了印度板塊和歐亞板塊碰撞后板塊邊緣變形與地幔對流之問的耦合關(guān)系，并指出地殼厚度變化是造成應(yīng)力場變化的主要因素，走滑斷裂只是部分地緩解了這種應(yīng)力場的變化。

但是上述研究中也存在明顯的不足之處，如考慮的只是小變形，忽略了溫度、重力、沉積、剝蝕等因素的影響等。但實際的地質(zhì)過程中，一方面地質(zhì)構(gòu)造過程屬于大變形（幾何非線性），并不滿足上述有限元法變形（小變形）的前提條件，此外所有的地質(zhì)因素是相互影響耦合在一起的，單獨考慮某一個和幾個因素必然會影響模擬的結(jié)果，甚至出現(xiàn)大的偏差。因此Houseman和England［10－11］對大變形，大旋轉(zhuǎn)的速度場和應(yīng)力場等進行了較全面的分析，并提出了新的有限元模擬方法。Zhang等用有限元的模擬手段成功地解釋了澳大利亞東部的巖石圈構(gòu)造和異常應(yīng)力場的成因［12］，并用變形一流體一熱力學(xué)一化學(xué)反應(yīng)全耦合的方法實現(xiàn)了在構(gòu)造控礦理論研究和礦產(chǎn)勘探中的應(yīng)用。

實際上以加拿大人Beaumont為核心的研究小組已經(jīng)建立了一個完整的數(shù)值模擬體系，包括理論方法，計算技術(shù)以及各類計算模［13－15］等。他們采用不同的邊界條件進行計算，對斜向匯聚碰撞和巨型造山帶的應(yīng)力場進行了模擬分析，模擬介質(zhì)包括彈塑性，剛塑性，粘彈性和粘彈塑性等，同時還考慮了溫度場、應(yīng)力場和流體場的耦合問題，并較好的解釋了新英格蘭阿巴拉契亞、喜馬拉雅、阿爾卑斯等著名造山帶的應(yīng)力異常問題［5，13－16］。此外，Beaumont C、Batt G E等人在研究擠壓構(gòu)造帶時，還全面的分析了山體抬升引起的附加重力、剝蝕過程對質(zhì)量的重新分配，以及由此造成的應(yīng)力應(yīng)變的變化。并給出了構(gòu)造變形與剝蝕過程的耦合數(shù)值求解方法，其中計算剝蝕量的地表剝蝕模型包括了長距離的河流搬運作用和短距離的巖崩、滑坡、沖刷等擴散作用的網(wǎng)格地形模型［17－18］。

國內(nèi)一些地質(zhì)學(xué)者在有限元數(shù)值模擬方面也做了很多重要工作。其中青藏高原由于其在全球構(gòu)造中獨特性，以及有限元數(shù)值模擬方法在研究構(gòu)造應(yīng)力場演化機制方面的優(yōu)越性，使得青藏高原一直是國內(nèi)專家學(xué)者的研究熱點和重點，同時也是世界上數(shù)值模擬方法應(yīng)用最多的一個區(qū)域之一。

石耀林［19］等用有限元方法研究分析青藏高原熱演化中長期變形粘性流動，是國際上最早對陸陸碰撞非彈性力學(xué)分析的論文之一。熊熊、傅容珊等［20－21］假設(shè)地幔是粘滯度為常數(shù)的二維牛頓流體，模擬了青藏高原增厚大陸巖石層熱邊界被對流地幔剝離并為軟流圈物質(zhì)替代的動力學(xué)過程，并推斷高原地殼和巖石層在巖石層增厚和剝離以前就很熱，指出擠壓隆升受多種因素的影響，而且高原隆升是不均勻演化的過程。鄭勇等［22］采用有限元方法模擬了近20萬年來青藏高原巖石圈形變演化過程，探討了印度一歐亞大陸的碰撞對中國大陸巖石層形變和應(yīng)力場的影響以及它們與強地震活動性的關(guān)系。段巖等［23］人結(jié)合深部地質(zhì)與地球物理資料，對喜馬拉雅構(gòu)造帶中的札達盆地斷裂的構(gòu)造應(yīng)力場進行了模擬，計算結(jié)果表明札達盆地的演化明顯受盆地兩側(cè)邊界斷裂的控制，札達盆地是在整體南北向擠壓應(yīng)力的作用下不同塊體差異隆升作用的結(jié)果。曹建玲等［24］根據(jù)GPS測量結(jié)果，利用三維球坐標系下的Maxwell體有限元模型模擬了青藏高原在印度板塊推擠下的變形。他們認為柔軟下地殼的存在使整個青藏高原在印度板塊的推擠下表現(xiàn)為整體抬升，高原南緣和喜馬拉雅東構(gòu)造結(jié)抬升迅速，而周邊地塊下地殼相對較硬而封閉，僅高原東南部存在高溫柔軟的通道，使得高原整體隆起達到一定高度后，下地殼和軟流層的物質(zhì)向東、東南流動，并拖曳上地殼作類似運動，形成繞喜馬拉雅東構(gòu)造結(jié)的順時針轉(zhuǎn)動。

此外，有限元數(shù)值模擬方法在地震成因以及孕震過程等方面也得到了非常廣泛的應(yīng)用，這方面的研究已有許多報道。早在1972年，Lysmer［25］就將有限元法應(yīng)用于地震數(shù)值模擬之中，Mat suura等［26］用一個簡化了的巖石圈－軟流圈耦合模型數(shù)值模擬了橫穿板塊邊界的構(gòu)造載荷情況，研究結(jié)果表明地震斷層的應(yīng)力累積一部分來源于其基底載荷（基底巖石圈的粘性阻力），一部分來源于其邊緣載荷（斷層水平邊緣上的位錯累積）；而Hashimoto等［27－28］通過綜合考慮粘彈性滑動響應(yīng)函數(shù)、隨時間與滑動變化的斷層本構(gòu)關(guān)系，以及給定一個板塊運動驅(qū)動力建立了仿真3D模型，計算模擬了板塊邊界地震發(fā)生的整個過程。在國內(nèi)，王仁等人［29］是較早開展這方面的研究人員。朱守彪［30］以2004年發(fā)生過特大地震的蘇門答臘俯沖帶為例，模擬了俯沖帶上俯沖板片與上伏板塊之間的閉鎖、解鎖、滑動到再閉鎖這一準周期性過程，成功的模擬了地震的孕育、發(fā)生過程。實際上目前有限元數(shù)值模擬方法孕震過程與前兆機理研究中的應(yīng)用已經(jīng)取得了很大進展：地質(zhì)模型從二維到三維；本構(gòu)關(guān)系從單一的彈（塑）性到粘－彈性、粘－塑性，從線性到非線性；模型中地殼（巖石圈）的分層從簡單的上、下兩層到細劃分多層［1－2，31－34］；計算方法從原來的單機串行計算到目前的CPU矩陣并行計算甚至是云計算等［35］（更多的文獻可在有限元方法在地球科學(xué)中應(yīng)用的評述性文章中查到［36］）。劉啟元、馬宗晉［36］等人曾經(jīng)指出，地震預(yù)報研究應(yīng)充分吸收和借鑒氣象預(yù)報的成功經(jīng)驗，以GPS為代表的空間對地觀測技術(shù)，巖石圈巨型高分辨率寬頻帶流動地震臺陣觀測技術(shù)以及數(shù)值模擬技術(shù)已經(jīng)為地震數(shù)值預(yù)報研究提供了前所未有的技術(shù)基礎(chǔ)。當(dāng)前地震數(shù)值預(yù)報突破的目標無疑應(yīng)首先主要集中在地震的中期或中短期預(yù)報。以地震數(shù)值預(yù)報為目標的GPS陣列地殼形變連續(xù)觀測、高分辨率地殼上地幔結(jié)構(gòu)探測、地殼動力學(xué)、地震孕育和破裂過程的理論、模擬試驗和實際觀測、數(shù)據(jù)同化和計算軟件的開發(fā)應(yīng)成為今后研究的重點。地震數(shù)值預(yù)報研究必將極大地促進中國地震科學(xué)基礎(chǔ)研究和地震預(yù)報的進一步發(fā)展。

盡管國內(nèi)的研究人員已將構(gòu)造應(yīng)力場的有限元數(shù)值模擬技術(shù)廣泛地應(yīng)用在了地震、地質(zhì)、油田開發(fā)等諸多行業(yè)，并且在計算方法、模型處理做了很多有益的嘗試，但是在模型建立、技術(shù)處理等方面與國際上還存在著一定的差距，并且目前還沒有開發(fā)出能夠獲得同行普遍認同的地學(xué)專業(yè)有限元軟件。很多研究者還是采用的線性模型，用一些通用的有限元軟件（ANSYS，ALGOR等）來處理復(fù)雜地質(zhì)問題，存在較大的誤差，致使一些模擬結(jié)果難以令人信服，此外也較少考慮時間以及地表過程（沉積、剝蝕、滑坡等因素）的影響等。

3 相關(guān)的技術(shù)方法

由于數(shù)值模擬結(jié)果存在多解性，因而建立一個合適的計算模型是模擬構(gòu)造應(yīng)力場中十分重要的一個環(huán)節(jié)，從而減少多解性。它包括建立地質(zhì)模型和數(shù)值（物理）模型兩個階段：在建立數(shù)值模型階段，需要定義模型的幾何形狀，選擇表達模型物理行為特征的數(shù)學(xué)描述方法，并確定適合該地質(zhì)模型有限單元類型；在建立地質(zhì)模型階段，需要根據(jù)觀測或推測的邊界條件、實驗參數(shù)，來建立實驗研究的初始條件和邊界條件，并在模擬過程中不斷地修正、更改實驗參數(shù)、初始條件和邊界條件進行模擬實驗，逐步使實驗結(jié)果趨近于觀察到或科學(xué)推測的地質(zhì)構(gòu)造模型，從數(shù)學(xué)理論和方法上來證實或證偽已有的地質(zhì)模型，從而建立起自己的符合區(qū)域地質(zhì)實況的地質(zhì)模式。因此，綜合分析已有資料，構(gòu)建合理的地質(zhì)模型，選擇不同的數(shù)學(xué)模型進行實驗研究，并從不同方面選擇實驗研究的初始條件、邊界條件及實驗參數(shù)，從地質(zhì)、地球物理、地球化學(xué)等方面提供約束條件，限制實驗結(jié)果的多解性，并使實驗結(jié)果趨近于地質(zhì)現(xiàn)狀，是數(shù)值模擬的主要研究內(nèi)容與方法。

此外，在構(gòu)造應(yīng)力場的模擬過程中，如果考慮地殼深部溫度場的影響，還需要引入熱動力學(xué)方程，將熱動力學(xué)方程和三大方程（平衡方程、本構(gòu)方程和幾何方程）進行聯(lián)合求解，即溫度場和應(yīng)力場的耦合處理。Fullsack［13］、Beaumont［15］、Batt［18］等對熱一力學(xué)耦合問題給出了相應(yīng)的結(jié)果。

采用有限元數(shù)值方法模擬構(gòu)造應(yīng)力場的過程中，斷層的處理是其中至關(guān)重要的一環(huán)。因為有限元法較適合分析連續(xù)介質(zhì)的力學(xué)問題，對于地質(zhì)構(gòu)造中涉及到的節(jié)理、斷層等不連續(xù)現(xiàn)象需要特殊的處理方法。一般斷層問題可以采用以下幾種辦法解決：一種是goodman單元法［37］，這種單元比較簡單，但是不適合算大的位移（變形）；一種是滑移節(jié)點法，這種方法常用在模擬地震的瞬間演化；另外是拉格朗日乘子法，這種方法可以直接求出斷層面上的接觸力，因為里面涉及一個滿矩陣的求逆，所以計算規(guī)模比較小，但它的最致命的弱點是不容易收斂。對于地質(zhì)上大的構(gòu)造變形問題（如造山運動），如果本構(gòu)關(guān)系選擇密律流體本構(gòu)，斷層的處理會相對容易些，只需要把斷層單元劃分薄一點，然后給斷層單元分配一個比較低的粘度系數(shù)就可以了。并且選擇密律流體本構(gòu)也優(yōu)于粘彈性本構(gòu)，因為大變形都涉及長時間的應(yīng)力松弛現(xiàn)象。這也是現(xiàn)在比較通行一種處理方法［17］，熊熊、曹建玲等［20－24］均采用了該方法處理構(gòu)造應(yīng)力場中的斷層。

當(dāng)要考慮地表過程造成的影響時，例如由于山體抬升引起的附加重力、剝蝕沉積過程對質(zhì)量的重新分配，以及由此造成的應(yīng)變應(yīng)力的變化等，需要采用特殊的耦合處理方法，有的還需要在剖分網(wǎng)格的時候進行局部加密。Beaumont［15－16］等給出了構(gòu)造變形與剝蝕過程的耦合數(shù)值求解方法。但是該耦合過程不包含基本變量的耦合求解，只是在每一個時間步的構(gòu)造抬升計算完成之后，根據(jù)剝蝕模型計算剝蝕量，再對地形進行修正。其中計算地殼變形抬升的模型為受擠壓作用的剛塑性或粘性平面應(yīng)變模型；計算剝蝕量的地表剝蝕模型為包括了長距離的河流搬運作用和短距離的巖崩、滑坡、沖刷等擴散作用的網(wǎng)格地形模型。

4 問題與討論

用有限元法模擬構(gòu)造應(yīng)力場不僅是當(dāng)前定量研究地質(zhì)力學(xué)的熱點和焦點，也是難點之一。雖然我國在這方面的研究工作早已經(jīng)開始，但是力量薄弱，只能說還處于探索階段，與國外同行還有不小的差距，尤其是多場耦合的三維非線性問題的模擬。此外構(gòu)造應(yīng)力場的模擬涉及的學(xué)科領(lǐng)域較多，需要考慮固體力學(xué)、流體力學(xué)和熱力學(xué)等各方面的信息，對于沉積、剝蝕等地表因素的處理也都需要用不同的方程表示，并需要將上述各方面的控制方程耦合在一起，這樣就極大的增加了求解難度。要很好的解決這些問題需要大的學(xué)科交叉和科學(xué)家的合作，尤其需要數(shù)學(xué)家、物理學(xué)家的參與，并盡早的開發(fā)出適用于地學(xué)的有限元專業(yè)軟件。

［1］ 章純.中國東部地區(qū)地震活動與構(gòu)造應(yīng)力場關(guān)系的有限元數(shù)值模擬［J］.西北地震學(xué)報，2007，29（3）：230－234.

［2］ 王愛國，袁道陽，梁明劍.蘭州盆地最大潛在地震變形數(shù)值模擬［J］.西北地震學(xué)報，2008，30（3）：232－238.

［3］ 王建，葉正仁.地幔對流對全球巖石圈應(yīng)力產(chǎn)生與分布的作用［J］.地球物理學(xué)報，2005，48（3）：584－590.

［4］ 陸遠忠，葉金鐸，蔣淳，等.中國強震前兆地震活動圖像機理的三維數(shù)值模擬研究［J］.地球物理學(xué)報，2007，50（2）：499－ 508.

［5］ Beaumont C，Jamieson R A.Himalayan tectonics explained by extrusion of a low－viscosity crustal channel coupled to focused surface denudation［J］.Nature，2001，414：738－742.

［6］ 殷有泉.計算固體力學(xué)方法［M］.北京.科學(xué)出版社，2003：121－148.

［7］ Bott M H，Dean D S.Stress systems at young continental margins［J］.Nature，1972，235：23－25.

［8］ Tappennier P，Molnar P.Slip line field theory and large－scale continental tectonics［J］.Nature，1976，264：319－324.

［9］ England P C，Mackenzie D A.Thin viscous sheet model for continental deformation［J］.Journal of Geophysics Research，1982，70：295－321.

［10］ England P，Houseman G.Finite strain calculation of continental deformation，2.Comparison with the India－Asia collision zone［J］.Journal of Geophysical Research，1986，91（B3）：3664－3676.

［11］ England P，Houseman G.Role of Lithospheric strength hetemgeneities in the tectonics of Tibet and neighboring regions［J］.Nature，1985，315：297－301.

［12］ Zhang Y，Scheibner E，OrdA.Numerical modeling of crustal stresses in the eastern Australian passive margin［J］.Austra1ian Journal of Earth Sciences，1996，43：161－175.

［13］ Fullsack P.An arbitrary Lagrangian－Eulerian formulation for creeping flows and its application in tectonic models［J］.Geophysical Journal International，1995，120：1－23.

［14］ Ellis S，Beaumont C.Models of convergent boundary tectonics：implications for the interpretation of Lithoprobe data［J］.Earth Sciences，1999，36：1711－1741.

［15］ Beaumont C，Ellis S，Hamilton J.Mechanical model for subduction－collision tectonics of Alpine－type compressional orogens［J］.Geology，1996，24：675－678.

［16］ Beaumont C，Muoz J A，Hamilton，et al..Factors controlling the Alpine evolution of the central Pyrenees inferred from a comparison of observations and geodynamical models［J］.Journal of Geophysical Research，2000，105：8121－8145.

［17］ Willett S，Beaumont C，F(xiàn)ullsack P.Mechanical model for the tectonics of doubly vergent compressional orogens［J］.Geology，1993，21：371－374.

［18］ Batt G E，Braun J.On the thermo－mechanical evolution of compressional orogens［J］.Geophysical Journal International，1997，128：364－382.

［19］ Chi－yue Wang，Yao－lin Shi.Dynamic uplift of the Himalaya［J］.Nature，1982，298：553－556.

［20］ 熊熊，傅容珊，許厚澤，等.增厚大陸巖石圈熱邊界層對流剝離的數(shù)值模擬［J］.地球物理學(xué)報，1998，41（增刊）：33－40.

［21］ 傅容珊，徐耀民，黃建華，等.青藏高原擠壓隆升過程的數(shù)值模擬［J］.地球物理學(xué)報，2000，43（3）：346－355.

［22］ 鄭勇，傅容珊，熊熊.中國大陸及周邊地區(qū)現(xiàn)代巖石圈演化動力學(xué)模擬［J］.地球物理學(xué)報，2006，49（2）：4153－427.

［23］ 段巖，孟憲剛，邵兆剛，等.西藏札達盆地控盆斷裂有限元數(shù)值模擬［J］.地質(zhì)通報，2008，27（10）：12－16.

［24］ 曹建玲，石耀霖，張懷，等.青藏高原GPS位移繞喜馬拉雅東構(gòu)造結(jié)順時針旋轉(zhuǎn)成因的數(shù)值模擬［J］.科學(xué)通報，，2009，54（2）：224－234.

［25］ Lysmer J.A Finite element method for seismology［J］.Journal of Computational Physics，1972，（11）：181－216.

［26］ Mat suura M，Sato T.Loading mechanism and scaling relation of large interplate earthquakes［J］.Tectonophysics，1997，277：189－198.

［27］ Hashimoto C，Mat suura M.3Dsimulation of earthquake generation cycles and evolution of fault constitutive properties［J］.Pure Appl Geophys，2002，159：2175－2199.

［28］ Aochi H，Mat suura M.Slip and time2dependent fault constitutive law and it s significance in earthquake generation cycles［J］.Pure.Appl.Geophys.，2002，159：2 02922 046.

［29］ 王仁，孫荀英，蔡永恩.華北地區(qū)近年地震序列的數(shù)學(xué)模擬［J］.中國科學(xué)（B輯），1982，18（8）：715－753.

［30］ 朱守彪，邢會林，謝富仁，等.地震發(fā)生過程的有限單元法模擬—以蘇門答臘俯沖帶上的大地震為例［J］.地球物理學(xué)報，2008，51（2）：460－468.

［31］ 王愛國，石玉成，柳煜.黃河黑山峽大柳樹壩址區(qū)最大潛在地震變形及地震應(yīng)力模擬預(yù)測［J］.西北地震學(xué)報，2007，29（4）：314－318.

［32］ 文彥君.延懷盆地設(shè)定地震淺析［J］.西北地震學(xué)報，2008，30（2）：159－162

［33］ 劉海明，陶夏新，孫曉丹，等.馬銜山北緣斷裂西段6.5級地震對蘭州市及周邊地區(qū)的影響［J］.西北地震學(xué)報，2008，30（3）：227－231.

［34］ 何麗君，石玉成，楊惠林，等.地震動作用下黃土邊坡穩(wěn)定性分析［J］.西北地震學(xué)報，2009，31（2）：142－147.

［35］ 張懷，吳忠良，張東寧，等.虛擬川滇—基于千萬網(wǎng)格并行有限元計算的區(qū)域強震演化過程數(shù)值模型設(shè)計和構(gòu)建［J］.中國科學(xué)（D輯），2009，39（3）：260－270.

［36］ 劉啟元，吳建春.論地震數(shù)值預(yù)報—關(guān)于我國地震預(yù)報研究發(fā)展戰(zhàn)略的思考［J］.地學(xué)前緣，2003，10（8）：217－224.

［37］ Goodman R E.Methods of geological engineering in discontinuous rocks［M］.Paul：West Publishing Co.，1976.

Modeling of Tectonic Stress Field——the Theroy，Method and Related Research Progress of the Finite Element Method

ZHANG Sheng－li1，2
（1.School of Information，Wuyi University，Guangdong Jiangmen 529020，China；2.Guangzhou Institute of Geochemistry，Chinese Academy of Sciences，Guangzhou 510640，China）

P315.12

A

1000－0844（2010）04－0405－06

2009－07－20

中國科學(xué)院知識創(chuàng)新工程項目（KZCXZ－SW－117）；廣東高校優(yōu)秀青年創(chuàng)新人才培養(yǎng)計劃項目（LYM10128）

張勝利（1979－），男（漢族），東單縣人，博士，主要從事數(shù)值模擬及數(shù)據(jù)挖掘方面研究.