張勝利
(1.五邑大學信息學院,廣東江門 529020;2.中科院廣州地化所,廣東廣州 510640)
構造應力場模擬
——有限元理論、方法和研究進展①
張勝利1,2
(1.五邑大學信息學院,廣東江門 529020;2.中科院廣州地化所,廣東廣州 510640)
采用有限元數值模擬方法對構造地質問題進行描述和定量化求解是當前地質學領域的研究的一個熱點,在近10年以來取得了重要進展,形成了比較完整的理論和技術體系,并在一些典型的地質構造帶獲得了重要的研究成果。本文以有限元數值模擬方法理論作為出發(fā)點,總結分析了國內外有限元數值模擬方法在構造應力場領域的研究進展情況和技術方法,并討論了其目前存在的問題和未來發(fā)展方向。
構造應力場;數值模擬;有限單元法
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
地殼中的各種地質構造都是巖石受力發(fā)生變形的產物,它們的產生和發(fā)展必然也受力學規(guī)律的支配。因此研究構造應力場特征是闡明和解釋各種地質構造的產生、分布和演化規(guī)律的重要途徑。但由于地質構造演化過程的特殊性,采用傳統(tǒng)的地質-地球物理-地球化學的方法研究應力場特征顯得異常困難。而物理模擬雖然可以幫助人們理解構造變形過程,但是存在著嚴重的時空尺度的局限性,并且一些因素在實驗室條件下尚不能獲得,尤其是在伴隨著高溫高壓的地球深部構造帶。
有限單元法又稱有限元法,是上世紀50年代中期發(fā)展起來的一種數值計算方法,最初主要用于結構和應力分析,60年代后期被引入到地質構造分析中,由于其巨大優(yōu)越性而被迅速推廣,已成功的解決了很多地學的難題。我國是從70年代中期開始將有限元數值模擬技術應用在地學領域的。隨著有限元數值模擬技術的出現,在綜合利用其它方法(如地質、地球物理、地球化學)研究成果的基礎之上,采用這類數值模擬技術,可以建立不受時空限制的各種地質模型,幫助研究人員比較直觀地認識和理解構造應力場的特征和演化歷史,并可以對其演化趨勢做出預測?,F在該技術除了在基礎地質方面得到廣泛應用外,在天文、環(huán)境、礦產勘探、地震預測等方面也取得了眾多的研究成果[1-3]。
隨著計算機計算速度的提高,構造應力場的數值模擬方法也獲得了很大的進展。國內外的研究結果表明,目前已經可以在構造應力場的模擬過程,在不同的本構關系下,將流體、溫度、剝蝕、斷層等因素考慮進去,進行二維和三維的數值模擬[4-5]。但是也存在很多的問題,例如模擬的主觀性(邊界條件、模型等)、模擬結果的多解性和可靠性等。
有限元法的基本原理是根據三大守恒定律(質量、動量、能量)建立三個偏微分方程:平衡方程、本構方程和幾何方程。然后將擬分析的連續(xù)體假想地分為有限個單元,離散后的單元與單元之間利用單元節(jié)點相互連接起來。單元內部點的待求量可由單元節(jié)點量通過選定的插值函數求得。由于單元形狀簡單,易于用平衡關系或能量關系建立節(jié)點量之間的方程式。然后將各個單元方程“組裝”在一起形成總體代數方程組,計入邊界條件后即可對方程組求解。單元劃分越細,計算結果就越精確。
采用有限元法模擬構造應力場一般包括一下幾個步驟:
(1)在詳細分析研究目標地質資料的基礎之上建立計算模型,這也是數值模擬試驗成敗的關鍵。
(2)單元剖分和選擇插值函數。根據模型的幾何特性、載荷情況及所要求的變形點,將模型離散化,單元大小應當可以給出有用的結果,但又要盡可能的節(jié)省計算量。
(3)確定應變位移和應力應變關系。應力與應變之間關系由本構方程確定,不同的本構關系對應不同的物理介質。目前構造應力場模擬中較常用的本構關系有4種:線彈性、彈塑性、粘彈性和彈粘塑性本構關系[6]。
其中彈塑性本構一般用于模擬地殼淺部的脆性層,粘彈性和彈粘塑性本構多用于模擬地殼深部。也可以在同一個模型中采用多種本構關系。
(1)線彈性本構關系:
其中σij是應力張量;eij是應變張量;λ和μ是拉梅系數;θ是體積模量。線彈性本構一般用于模擬小范圍的地殼淺部沉積層的應力場。
(2)彈塑性本構關系:
其中Dijkl為彈性系數;εkl為應變張量;dλ是由一致性條件確定的尺度因子;F是用應變空間表示的屈服函數。彈塑性本構一般用于模擬地殼淺部脆性層。
(3)粘彈性本構關系(以Kelvin體為例):其中E為彈性模量;εij為應變速率張量;η是粘性系數。式(3)也常被稱為線性流變方程。
(4)彈粘塑性本構關系:
其中γ是控制塑性流動速率的流度參數,是時間和粘塑性應變張量的不變量的函數;f是屈服函數;f0是一個使表達式無量綱化的、正的參考值。公式(2)、(3)也可以視作本式的特殊形式。
構造應力場模擬和板塊碰撞、大陸變形數值模擬緊密相關。Bott[7]對板緣力驅動下板內應力應變的演化模擬可以說是這方面的開山之作。1976年Tapponnier[8]等人用理想剛塑性體所產生的變形來模擬亞洲大陸自新生代以來所經歷的大尺度碰撞變形和斷層的產生,并為板塊內部強烈變形提供了有力證據,從而使得應力場的數值模擬方法引起了人們的廣泛關注。England[9]等在1982年首次將粘性薄層流變模型運用到以印度板塊和歐亞板塊的碰撞為代表的大陸變形帶,并考慮了地殼增厚造成的應力變化,對青藏高原的擠壓隆升演化進行了數值模擬;England和Houseman等[10-11]也成功模擬了印度板塊和歐亞板塊碰撞后板塊邊緣變形與地幔對流之問的耦合關系,并指出地殼厚度變化是造成應力場變化的主要因素,走滑斷裂只是部分地緩解了這種應力場的變化。
但是上述研究中也存在明顯的不足之處,如考慮的只是小變形,忽略了溫度、重力、沉積、剝蝕等因素的影響等。但實際的地質過程中,一方面地質構造過程屬于大變形(幾何非線性),并不滿足上述有限元法變形(小變形)的前提條件,此外所有的地質因素是相互影響耦合在一起的,單獨考慮某一個和幾個因素必然會影響模擬的結果,甚至出現大的偏差。因此Houseman和England[10-11]對大變形,大旋轉的速度場和應力場等進行了較全面的分析,并提出了新的有限元模擬方法。Zhang等用有限元的模擬手段成功地解釋了澳大利亞東部的巖石圈構造和異常應力場的成因[12],并用變形一流體一熱力學一化學反應全耦合的方法實現了在構造控礦理論研究和礦產勘探中的應用。
實際上以加拿大人Beaumont為核心的研究小組已經建立了一個完整的數值模擬體系,包括理論方法,計算技術以及各類計算模[13-15]等。他們采用不同的邊界條件進行計算,對斜向匯聚碰撞和巨型造山帶的應力場進行了模擬分析,模擬介質包括彈塑性,剛塑性,粘彈性和粘彈塑性等,同時還考慮了溫度場、應力場和流體場的耦合問題,并較好的解釋了新英格蘭阿巴拉契亞、喜馬拉雅、阿爾卑斯等著名造山帶的應力異常問題[5,13-16]。此外,Beaumont C、Batt G E等人在研究擠壓構造帶時,還全面的分析了山體抬升引起的附加重力、剝蝕過程對質量的重新分配,以及由此造成的應力應變的變化。并給出了構造變形與剝蝕過程的耦合數值求解方法,其中計算剝蝕量的地表剝蝕模型包括了長距離的河流搬運作用和短距離的巖崩、滑坡、沖刷等擴散作用的網格地形模型[17-18]。
國內一些地質學者在有限元數值模擬方面也做了很多重要工作。其中青藏高原由于其在全球構造中獨特性,以及有限元數值模擬方法在研究構造應力場演化機制方面的優(yōu)越性,使得青藏高原一直是國內專家學者的研究熱點和重點,同時也是世界上數值模擬方法應用最多的一個區(qū)域之一。
石耀林[19]等用有限元方法研究分析青藏高原熱演化中長期變形粘性流動,是國際上最早對陸陸碰撞非彈性力學分析的論文之一。熊熊、傅容珊等[20-21]假設地幔是粘滯度為常數的二維牛頓流體,模擬了青藏高原增厚大陸巖石層熱邊界被對流地幔剝離并為軟流圈物質替代的動力學過程,并推斷高原地殼和巖石層在巖石層增厚和剝離以前就很熱,指出擠壓隆升受多種因素的影響,而且高原隆升是不均勻演化的過程。鄭勇等[22]采用有限元方法模擬了近20萬年來青藏高原巖石圈形變演化過程,探討了印度一歐亞大陸的碰撞對中國大陸巖石層形變和應力場的影響以及它們與強地震活動性的關系。段巖等[23]人結合深部地質與地球物理資料,對喜馬拉雅構造帶中的札達盆地斷裂的構造應力場進行了模擬,計算結果表明札達盆地的演化明顯受盆地兩側邊界斷裂的控制,札達盆地是在整體南北向擠壓應力的作用下不同塊體差異隆升作用的結果。曹建玲等[24]根據GPS測量結果,利用三維球坐標系下的Maxwell體有限元模型模擬了青藏高原在印度板塊推擠下的變形。他們認為柔軟下地殼的存在使整個青藏高原在印度板塊的推擠下表現為整體抬升,高原南緣和喜馬拉雅東構造結抬升迅速,而周邊地塊下地殼相對較硬而封閉,僅高原東南部存在高溫柔軟的通道,使得高原整體隆起達到一定高度后,下地殼和軟流層的物質向東、東南流動,并拖曳上地殼作類似運動,形成繞喜馬拉雅東構造結的順時針轉動。
此外,有限元數值模擬方法在地震成因以及孕震過程等方面也得到了非常廣泛的應用,這方面的研究已有許多報道。早在1972年,Lysmer[25]就將有限元法應用于地震數值模擬之中,Mat suura等[26]用一個簡化了的巖石圈-軟流圈耦合模型數值模擬了橫穿板塊邊界的構造載荷情況,研究結果表明地震斷層的應力累積一部分來源于其基底載荷(基底巖石圈的粘性阻力),一部分來源于其邊緣載荷(斷層水平邊緣上的位錯累積);而Hashimoto等[27-28]通過綜合考慮粘彈性滑動響應函數、隨時間與滑動變化的斷層本構關系,以及給定一個板塊運動驅動力建立了仿真3D模型,計算模擬了板塊邊界地震發(fā)生的整個過程。在國內,王仁等人[29]是較早開展這方面的研究人員。朱守彪[30]以2004年發(fā)生過特大地震的蘇門答臘俯沖帶為例,模擬了俯沖帶上俯沖板片與上伏板塊之間的閉鎖、解鎖、滑動到再閉鎖這一準周期性過程,成功的模擬了地震的孕育、發(fā)生過程。實際上目前有限元數值模擬方法孕震過程與前兆機理研究中的應用已經取得了很大進展:地質模型從二維到三維;本構關系從單一的彈(塑)性到粘-彈性、粘-塑性,從線性到非線性;模型中地殼(巖石圈)的分層從簡單的上、下兩層到細劃分多層[1-2,31-34];計算方法從原來的單機串行計算到目前的CPU矩陣并行計算甚至是云計算等[35](更多的文獻可在有限元方法在地球科學中應用的評述性文章中查到[36])。劉啟元、馬宗晉[36]等人曾經指出,地震預報研究應充分吸收和借鑒氣象預報的成功經驗,以GPS為代表的空間對地觀測技術,巖石圈巨型高分辨率寬頻帶流動地震臺陣觀測技術以及數值模擬技術已經為地震數值預報研究提供了前所未有的技術基礎。當前地震數值預報突破的目標無疑應首先主要集中在地震的中期或中短期預報。以地震數值預報為目標的GPS陣列地殼形變連續(xù)觀測、高分辨率地殼上地幔結構探測、地殼動力學、地震孕育和破裂過程的理論、模擬試驗和實際觀測、數據同化和計算軟件的開發(fā)應成為今后研究的重點。地震數值預報研究必將極大地促進中國地震科學基礎研究和地震預報的進一步發(fā)展。
盡管國內的研究人員已將構造應力場的有限元數值模擬技術廣泛地應用在了地震、地質、油田開發(fā)等諸多行業(yè),并且在計算方法、模型處理做了很多有益的嘗試,但是在模型建立、技術處理等方面與國際上還存在著一定的差距,并且目前還沒有開發(fā)出能夠獲得同行普遍認同的地學專業(yè)有限元軟件。很多研究者還是采用的線性模型,用一些通用的有限元軟件(ANSYS,ALGOR等)來處理復雜地質問題,存在較大的誤差,致使一些模擬結果難以令人信服,此外也較少考慮時間以及地表過程(沉積、剝蝕、滑坡等因素)的影響等。
由于數值模擬結果存在多解性,因而建立一個合適的計算模型是模擬構造應力場中十分重要的一個環(huán)節(jié),從而減少多解性。它包括建立地質模型和數值(物理)模型兩個階段:在建立數值模型階段,需要定義模型的幾何形狀,選擇表達模型物理行為特征的數學描述方法,并確定適合該地質模型有限單元類型;在建立地質模型階段,需要根據觀測或推測的邊界條件、實驗參數,來建立實驗研究的初始條件和邊界條件,并在模擬過程中不斷地修正、更改實驗參數、初始條件和邊界條件進行模擬實驗,逐步使實驗結果趨近于觀察到或科學推測的地質構造模型,從數學理論和方法上來證實或證偽已有的地質模型,從而建立起自己的符合區(qū)域地質實況的地質模式。因此,綜合分析已有資料,構建合理的地質模型,選擇不同的數學模型進行實驗研究,并從不同方面選擇實驗研究的初始條件、邊界條件及實驗參數,從地質、地球物理、地球化學等方面提供約束條件,限制實驗結果的多解性,并使實驗結果趨近于地質現狀,是數值模擬的主要研究內容與方法。
此外,在構造應力場的模擬過程中,如果考慮地殼深部溫度場的影響,還需要引入熱動力學方程,將熱動力學方程和三大方程(平衡方程、本構方程和幾何方程)進行聯合求解,即溫度場和應力場的耦合處理。Fullsack[13]、Beaumont[15]、Batt[18]等對熱一力學耦合問題給出了相應的結果。
采用有限元數值方法模擬構造應力場的過程中,斷層的處理是其中至關重要的一環(huán)。因為有限元法較適合分析連續(xù)介質的力學問題,對于地質構造中涉及到的節(jié)理、斷層等不連續(xù)現象需要特殊的處理方法。一般斷層問題可以采用以下幾種辦法解決:一種是goodman單元法[37],這種單元比較簡單,但是不適合算大的位移(變形);一種是滑移節(jié)點法,這種方法常用在模擬地震的瞬間演化;另外是拉格朗日乘子法,這種方法可以直接求出斷層面上的接觸力,因為里面涉及一個滿矩陣的求逆,所以計算規(guī)模比較小,但它的最致命的弱點是不容易收斂。對于地質上大的構造變形問題(如造山運動),如果本構關系選擇密律流體本構,斷層的處理會相對容易些,只需要把斷層單元劃分薄一點,然后給斷層單元分配一個比較低的粘度系數就可以了。并且選擇密律流體本構也優(yōu)于粘彈性本構,因為大變形都涉及長時間的應力松弛現象。這也是現在比較通行一種處理方法[17],熊熊、曹建玲等[20-24]均采用了該方法處理構造應力場中的斷層。
當要考慮地表過程造成的影響時,例如由于山體抬升引起的附加重力、剝蝕沉積過程對質量的重新分配,以及由此造成的應變應力的變化等,需要采用特殊的耦合處理方法,有的還需要在剖分網格的時候進行局部加密。Beaumont[15-16]等給出了構造變形與剝蝕過程的耦合數值求解方法。但是該耦合過程不包含基本變量的耦合求解,只是在每一個時間步的構造抬升計算完成之后,根據剝蝕模型計算剝蝕量,再對地形進行修正。其中計算地殼變形抬升的模型為受擠壓作用的剛塑性或粘性平面應變模型;計算剝蝕量的地表剝蝕模型為包括了長距離的河流搬運作用和短距離的巖崩、滑坡、沖刷等擴散作用的網格地形模型。
用有限元法模擬構造應力場不僅是當前定量研究地質力學的熱點和焦點,也是難點之一。雖然我國在這方面的研究工作早已經開始,但是力量薄弱,只能說還處于探索階段,與國外同行還有不小的差距,尤其是多場耦合的三維非線性問題的模擬。此外構造應力場的模擬涉及的學科領域較多,需要考慮固體力學、流體力學和熱力學等各方面的信息,對于沉積、剝蝕等地表因素的處理也都需要用不同的方程表示,并需要將上述各方面的控制方程耦合在一起,這樣就極大的增加了求解難度。要很好的解決這些問題需要大的學科交叉和科學家的合作,尤其需要數學家、物理學家的參與,并盡早的開發(fā)出適用于地學的有限元專業(yè)軟件。
[1] 章純.中國東部地區(qū)地震活動與構造應力場關系的有限元數值模擬[J].西北地震學報,2007,29(3):230-234.
[2] 王愛國,袁道陽,梁明劍.蘭州盆地最大潛在地震變形數值模擬[J].西北地震學報,2008,30(3):232-238.
[3] 王建,葉正仁.地幔對流對全球巖石圈應力產生與分布的作用[J].地球物理學報,2005,48(3):584-590.
[4] 陸遠忠,葉金鐸,蔣淳,等.中國強震前兆地震活動圖像機理的三維數值模擬研究[J].地球物理學報,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] 殷有泉.計算固體力學方法[M].北京.科學出版社,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,Fullsack 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] 熊熊,傅容珊,許厚澤,等.增厚大陸巖石圈熱邊界層對流剝離的數值模擬[J].地球物理學報,1998,41(增刊):33-40.
[21] 傅容珊,徐耀民,黃建華,等.青藏高原擠壓隆升過程的數值模擬[J].地球物理學報,2000,43(3):346-355.
[22] 鄭勇,傅容珊,熊熊.中國大陸及周邊地區(qū)現代巖石圈演化動力學模擬[J].地球物理學報,2006,49(2):4153-427.
[23] 段巖,孟憲剛,邵兆剛,等.西藏札達盆地控盆斷裂有限元數值模擬[J].地質通報,2008,27(10):12-16.
[24] 曹建玲,石耀霖,張懷,等.青藏高原GPS位移繞喜馬拉雅東構造結順時針旋轉成因的數值模擬[J].科學通報,,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ū)近年地震序列的數學模擬[J].中國科學(B輯),1982,18(8):715-753.
[30] 朱守彪,邢會林,謝富仁,等.地震發(fā)生過程的有限單元法模擬—以蘇門答臘俯沖帶上的大地震為例[J].地球物理學報,2008,51(2):460-468.
[31] 王愛國,石玉成,柳煜.黃河黑山峽大柳樹壩址區(qū)最大潛在地震變形及地震應力模擬預測[J].西北地震學報,2007,29(4):314-318.
[32] 文彥君.延懷盆地設定地震淺析[J].西北地震學報,2008,30(2):159-162
[33] 劉海明,陶夏新,孫曉丹,等.馬銜山北緣斷裂西段6.5級地震對蘭州市及周邊地區(qū)的影響[J].西北地震學報,2008,30(3):227-231.
[34] 何麗君,石玉成,楊惠林,等.地震動作用下黃土邊坡穩(wěn)定性分析[J].西北地震學報,2009,31(2):142-147.
[35] 張懷,吳忠良,張東寧,等.虛擬川滇—基于千萬網格并行有限元計算的區(qū)域強震演化過程數值模型設計和構建[J].中國科學(D輯),2009,39(3):260-270.
[36] 劉啟元,吳建春.論地震數值預報—關于我國地震預報研究發(fā)展戰(zhàn)略的思考[J].地學前緣,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
中國科學院知識創(chuàng)新工程項目(KZCXZ-SW-117);廣東高校優(yōu)秀青年創(chuàng)新人才培養(yǎng)計劃項目(LYM10128)
張勝利(1979-),男(漢族),東單縣人,博士,主要從事數值模擬及數據挖掘方面研究.