1、第九講第九講 巖石力學的數值模擬分析巖石力學的數值模擬分析 巖石力學在上個世紀的突出進展在于數值計算方法，即計算機的計算方法在巖土工程中的應用獲得了巨大的進展。近數十年來，我國已開展了大量的土木巖土工程。有些工程已屬于世界第一流的水平，如水電工程中已建立了大量的混凝土高壩和堆石壩，如世界第三高拱壩二灘拱壩、龍羊峽高壩。 現今可用于對巖體工程結構進行力學分析的數值方法多種多樣，每一種方法有其針對性和特點，對一個具體的問題用數值模擬方法進行分析時，應選擇一種最適合該問題的方法進行研究。 數值模擬方法的選擇，取決于研究對象即巖體工程結構的巖石力學性質和數值模擬的目的。 嚴格地講，巖石除具有彈性性質外

2、，還具有塑性性質和粘性性質，只不過在特定情況下，某種性質占主導地位而已。在巖體工程實踐中，硬巖及應力水平不甚高的中硬巖，其力學性質主要呈現為彈性或彈塑性；高應力環境下的軟巖，其力學性質主要呈現為塑性或粘塑性；對于服務時間較長的地下工程，巖石極軟或軟且應力水平很高，則在計算分析中不能忽視巖石的流變性質。 為了達到了解整個巖體工程系統的應力及變形規律的目的，各種數值方法均可采用，但以彈、塑性有限單元法或拉格朗日元法最為適宜。這兩種方法的單元劃分靈活，計算所需參數較少且易獲得，軟件也易于得到，成本較低。 局部工程結構的應力及變形分析，若巖石中硬以上，則各種方法均可采用；若巖石軟弱，則宜采用能進行大變

3、形分析的拉格朗日元法；若巖體可能發生非連續破壞，則宜采用離散單元法、非連續變形法或流形元法或大變形分析的拉格朗日元法。 要實現了解巖體結構的破壞特征及動態破壞過程的目的，則只能采用離散單元法、非連續變形法或流形元法，因為這些方法就是針對巖體介質的非連續性而提出的。 此外，對于一個具體問題，是進行平面分析或是進行三維分析，也需做出恰當選擇。嚴格地講，所有的問題都是三維的，但如果采用平面分析既能達到目的，計算結果誤差也不大，為了降低費用和快速方便起見，則以采用平面分析為宜。反之，則應采用三維程序代碼進行計算分析。 表1給出了幾種數值模擬方法所依據的基本原理、求解方式、離散化方法及其適用條件，可供選

4、擇模擬方法時參考。數值模擬數值模擬方法方法基本原理基本原理 求解方法求解方法離散方式離散方式適用條件適用條件有限單元有限單元法法最小勢能原理解方程組全區域劃分單元巖石中硬以上，小變形，巖體不會發生非連續性破壞如滑動、轉動、分離等邊界單元邊界單元法法Betti互等定理解方程組邊界上劃分單元同上離散單元離散單元法法牛頓運動定律顯式差分按結構弱面分布特征劃分單元巖石中硬以上，低應力水平，大變形，巖體沿弱面發生非連續性破壞非連續變非連續變形法形法最小勢能原理解方程組按主要結構弱面實際情況劃分單元大變形，巖體發生非連續性破壞數值流形數值流形法法最小勢能原理解方程組全區域劃分單元中硬以上巖體的連續或非連續

5、變形拉各朗日拉各朗日元法元法牛頓運動定律顯式差分全區域劃分單元巖石軟弱，大變形，巖體的破壞以變形為主表1 幾種數值模擬方法的對比近30年計算機發展很快，也推動了有限元的發展，并已應用于實際工程問題。但在遇到實際工程問題時，就有了困難，主要是實際工程量的DOF（自由度）太大。第一，工程對象幾何形狀復雜，如地下廠房的構造，尾水管內墩子，拱壩和地基往往也涉及很大的范圍。第二，在求解一些角度點，如鎮墩、排水孔結構等，自由度往往都很大。第三，在考慮巖體的開裂、巖體中的裂隙等，則要求更為龐大的自由度。這些都對計算能力提出了很高的要求，雖然現在的計算機的計算能力有了很大的提高，但是要完成更加復雜的計算任務，

6、為了提高計算的規模和減少計算時間，就需要進行并行計算。幾乎所有的大型計算機都是并行機。并行計算對于超大型的數值模型運算已經是一個必不可少的工具。并行計算在國外得到了很大的發展，自從美國國家航天局（NASA）的A.K.Noor在1975年發表第一篇關于有限元并行性計算的文章以來，有限元并行處理技術幾乎與并行計算機同步發展，并且由算法研究發展到了算法、軟件和硬件相結合的研究。在國內，并行計算還處于發展階段，而且硬件相對比較落后。在硬件上主要基于向量機、分布式并行機和共享存儲式并行機；在內容上，涉及范圍比較廣，但在系統性和深度上有待進一步發展，軟件的開發距離實際應用和商品化還很遠；專門針對有限元并行

7、計算的硬件研究很少。（一）（一）并行計算的基本概念并行計算的基本概念并行計算機并行計算機并行計算是指將計算任務分為n份，并將其分配給n個計算節點（計算節點數對網絡機群而言指工作站臺數；對大型并行機而言指CPU個數）同時進行計算。并行實際上包含了兩個概念：同時性（Simultaneity）和并發性（Concurrency）。同時性是指兩個或者多個事件在同一時刻發生；并發性是指兩個或者多個事件在同一時間間隔內發生。對于并行計算機，按照實際應用，可以分為如下的4類1：大規模并行機（Massively Parallel Processors，MPP）對稱多處理器系統（Symmetric MultiPr

8、ocessors，SMP）分布存儲系統（Distributed Shared Memory system，DSM）網絡工作站機群（Cluster Of Workstations，COW）在上面4種并行機中，COW和MPP具有類似的地方，實際上，當前MPP和COW之間的界限正在逐漸模糊。例如IBM SP2被視為MPP，但它卻具有機群結構。由于COW的性能價格比遠優于MPP，系統擴充方便，網絡連接形式多種多樣，所以COW是發展可擴放并行計算的主流趨勢。并行算法及并行程序開發并行算法及并行程序開發并行算法是區別于串行算法的另外一大類型的算法，它是適合于各種并行計算機上求解問題和處理數據的算法。并行程

9、序開發是對給定算法構造并行程序的活動，它要求算法設計者和計算機系統的體系結構的設計者進行廣泛的交互。并行程序的設計不僅編程困難，而且調試和分析更加困難。目前，并行程序設計在各方面都處于一個初級階段，它不支持一個成熟、穩定和通用的并行程序開發環境。并行開發環境指的是軟件環境，它對并行計算的影響比一般串行機要大得多。軟件對計算性能的影響的差別可達幾個數量級。并行開發環境包括兩部分內容：操作系統、通信平臺、編譯和調試工具以及性能測試軟件等。其中最主要的是操作系統和通信平臺。可以用于并行計算的操作系統主要有Windows（NT或者2000）和Unix，其中與Unix內核基本相同的自由軟件Linux在并

10、行計算中占優絕對的統治地位，而比較流行的通信平臺有P4、PVM、Express、PARMACS和MPI等。（二）有限元法的并行計算（二）有限元法的并行計算一般來說，有限元計算包括如下的步驟：（1）數據準備階段，也就是前處理階段，包括結構的離散化、單元信息和節點坐標的生成。（2）單元分析階段，包括計算單元剛度矩陣、單元荷載向量等。對于動力分析，還需要計算單元質量矩陣。（3）生成總體剛度矩陣和總荷載向量。（4）約束處理。（5）方程組求解，得到全部節點的位移。（6）后處理。在上面的計算中，（1）、（4）和（6）階段的計算花費的時間很短，而（2）、（3）和（5）階段則是有限元分析的關鍵所在。尤其是（5

11、），將占整個有限元分析的70的計算量。因此，并行計算主要是針對這3個步驟進行的。用于多機并行系統的有限元并行分析策略用于多機并行系統的有限元并行分析策略對于三維有限元的分析，如果在自帶存儲器的多機并行系統上，可以采用如下的兩種策略：SBS（Substruct By Substruct）策略 對要求解的結構采用某種剖分策略，比如，類似于子結構的方法，使得各個機器所承擔的部分剛度/質量矩陣可以進行“孤立”自由度（內部自由度）和“公共”自由度（界面自由度）分塊。“孤立”自由度將直接在“本地”進行“消去”，凝聚矩陣提交給“中心機”進行裝配并完成求解。然后由各個“節點機”完成消去未知量的換算。所采用的“

12、剖分策略”，一般以孤立自由度數目的總和最大為優化目標，兼顧及“同時完成任務”。采用這種策略，并行化計算是建立在子結構一級的水平上。EBE（Element By Element）策略 EBE策略的基本思想就是將一個總體矩陣的向量積轉化到一組單個矩陣的向量積的計算。在有限元計算中，總體矩陣就是由多個單元矩陣迭加而成，因此，非常適合使用EBE策略來進行求解。EBE策略是一種既能節省存儲量、又能適宜于并行計算而且容易實現的方法。靜力分析靜力分析在有限元的靜力分析中，主要是采用兩大類算法：直接并行算法和迭代并行算法。對于直接并行算法，一種是采用子結構直接并行算法，其基本思想就是SBS策略的基本思想。另外

13、一種直接解法是僅對占很大計算量的有限元方程組進行并行計算，其他計算步驟仍然采用串行算法。例如，并行求解方程組可以采用并行Gauss消去法、并行Cholesky分解等方法。周樹荃，鄧紹忠對直接并行算法進行了很深入的研究，并在YHI向量機上進行了大量的數值試驗。對于迭代并行算法，一般都采用同步運算或者基本完全異步控制，因此，迭代并行算法易于發揮并行機的向量處理功能。目前常用的算法有基于EBE策略的并行預處理共軛梯度法（Preconditioning Conjugate Gradient，PCG）、并行GMRES方法以及基于多色理論的超松弛迭代（Successive Over-Relaxation

14、Method，SOR）法等。動力分析動力分析TLDL動力響應分析 在這方面內容中，研究最多的是直接積分法，并取得了很大的成果。例如，高昀采用 分解方法，在Transputer及網絡機群環境下，實現了基于Newmark方法的的動力響應分析。結構的動力特性研究 在這一方面，相對于動力的響應分析來說，研究的相對較少。T.Hwang和I. D. Parsonsp網格法在MIMD并行計算機上求解特征值問題。國內在這一方面也做了大量的工作，主要是采用Lanczos方法對其進行求解。另外，徐甲同采用QIF分解方法對多處理機上矩陣特征值的并行計算進行了研究。發展方向發展方向采用可以達到完全并行的EBE策略，可

15、以完全突破大規模計算對串行計算機的內存和計算速度的限制。理論上而言，只要參與計算的節點足夠多，就可以對任意規模的數值計算問題進行求解。雖然采用EBE策略的預處理共軛梯度法是解決有限元并行計算很有效的方法，也是有限元并行計算發展的一個重要方向。尤其是在非線性問題中，EBE方法的效率是很高的。但是目前的各種預處理方法，不僅僅是EBE策略的預處理，就是對整體矩陣進行的預處理，都很難從理論上說明那一種更有效。如何給出一個有效的預處理方法仍然是一個需要解決的問題。另外，EBE方法在非線性問題中的收斂性也是需要進行研究的問題。 通常狀態巖體是被節理和斷層等構造切割成的不連續介質，這些不連續介質影響著巖體結

16、構的滑移和破壞，不連續介質數值模型即是為了解決這種地質問題而產生的。 非連續變形分析,即DDA(Discontinuous Deformation Analysis)，是由石根華博士提出的分析塊體系統運動和變形的一種數值方法。它以自然存在的節理面或斷層切割巖體形成不同的塊體單元，以各塊體的位移作為未知量，通過塊體間的接觸和幾何約束形成一個塊體系統。在塊體運動的過程中，嚴格滿足塊體間不侵入。將邊界條件和接觸條件等一同施加到總體平衡方程。總體平衡方程是由系統的最小勢能原理求得。求解方程組即可得到塊體當前時步的位移場、應力場、應變場及塊體間的作用力。反復形成和求解總體平衡方程式，即可得到多個時步后以

17、至最終平衡時的位移場及應力場等情況以及運動過程中各塊體的相對位置及接觸關系。因此，DDA方法可以模擬出巖石塊體的移動、轉動、張開、閉合等全部過程并據此判斷巖體的破壞程度、破壞范圍，從而對巖體的整體和局部穩定性作出正確的評價。DDA方法在三峽船閘高邊坡穩定分析以及錦屏的穩定分析中得到應用。1、近期的發展、近期的發展 塊體系統非連續變形分析是近年來發展起來的一種嶄新的數值計算理論。在滿足彈性理論的基本方程條件下能夠反映出巖體變形的不連續性，它既具有有限元理論基礎的嚴密性，又具有離散元法可計算塊體大位移的特點，是一種很有發展前途的數值計算方法。它以新穎的理論思想，嚴密科學的構思，廣泛應用前景相融為一

18、體。作者在三維理論的建立、完善的基礎上，結合三峽船閘高邊坡的穩定性分析進行了更深入細致的研究，全面發展完善和豐富了塊體系統非連續變形分析理論，取得了如下主要研究成果。（1）提出三維非連續變形分析理論，并對其理論進行完善。并從理論上論證和分析三維非連續變形分析理論的正確性，以12個變形參數來描述塊體系統的剛體運動和變形。建立總體平衡方程和各類子矩陣。（2）塊體間高效而正確的接觸判定是實現非連續變形分析計算的關鍵所在。提出利用以參考面法、參考線法來判定二維和三維塊體間的接觸。（3）參考面（線）是塊體邊界上的某些平面（線），不必重新生成。根據塊體間各平面的相對關系，建立了確定參考面的判據，參考面法不

19、僅適合對凸形塊體的接觸判定，而且對凹形塊體的接觸同樣適合。（4）進一步提出可采用形象法來尋找參考線、參考面，并給出了具體的尋找方法和步驟。利用形象法亦可迅速完成參考面的識別，該方法形象、直觀，這是一種高效率的搜尋方法。（5）建立了塊體靜動力學約束條件塊體間無侵入發生及不承受拉伸力的控制標準和位置鎖定原則。（6）實際工程中非連續面的大小是有限平面，它們在空間分布上具有很大的隨機性。針對這一情況，作者分析了目前處理節理網絡技術方面的不足，提出利用拓撲學單純同調原理解決塊體的形成問題。根據空間各有限平面，可以為有限元法、離散元法 塊體理論以及其他要求具有一般形狀的數值方法提供模型數據。該方法不僅可以

20、生成規則的塊體，而且還可以生成形狀十分復雜的多邊形、多面體，甚至是帶孔洞的塊體。（7）根據總剛矩陣 的特點，提出了一種高效三角分解方法圖解法。這種借助于圖解理論的非零元素存貯方法使得計算過程直觀而高效。通過圖形和連接矩陣的配合使用，使整個過程直觀明了，具有存貯量小、計算量少、舍入誤差小等優點。利用圖解技術，事先對塊體優化編號，確定存貯量，并使三角分解過程中產生的新元素個數最少，從而達到減少數據存貯和計算，提高解題效率的目的。K（8）根據非連續變形分析理論，編制和開發了相應的程序。并結合三峽船閘高邊坡穩定性分析實例進行了研究和分析，其計算結果與實際觀測結果大體接近。（9）從理論上解決了非連續變形

21、分析理論與有限元、邊界元及傳統極限平衡方法的耦合，并給出耦合原理及耦合的公式。（10）針對非連續變形分析法存在的問題，提出解決方案。時步的大小事關非連續變形分析理論計算結果的正確性，通過眾多的計算表明僅當時步小于0.005s時，計算結果是穩定可靠的。（11）針對非連續變形分析理論自身的特點及高邊坡滑坡位移復雜性，從理論上深刻分析、研究了邊坡滑坡時產生混沌的機理及過程，并提出防止邊坡進入混沌狀態的實際措施。（12）將不連續變形分析與AutoCAD動態技術結合，使塊體系統非連續變形分析法的計算實施過程形象、直觀，完成了塊體系統的變形或破壞全過程的模擬。2、不連續變形及其特點不連續變形及其特點 非連

22、續變形分析以各個塊體的位移作為求知量，通過塊體間的接觸和幾何約束形成一個塊體系統，塊體單元受非連續面的控制，在塊體運動的過程中單元之間可以接觸，也可以分離。在塊體運動過程中，滿足塊體間不侵入和不承受拉伸力的條件，總體平衡方程是由系統的最小勢能原理求得的。求解方程組即可得到塊體當前時步的位移場、應力場、應變場及塊體間的作用力，反復形成和求解總體平衡方程式，即可得到多外時步后塊體的變形，應力及位移情況，也可求得塊體系統最終達到平衡時應力場、位移場等情況以及運動過程中各塊體的相對位置及接觸關系。當然，非連續變形分析方法也和其它任何一種數值方法一樣，有其一定的適用范圍和局限性，這主要表現在：（1）巖體

23、的種類繁多，其性質極為復雜。在通常情況下，非連續巖體的力學表現為非均質、各向異性、非線性等復雜的特性，隨著應用工程環境的不同，使巖體數值計算有許多非確定因素。這對非連續變形分析在對巖體參數（ ， ）取值時帶來麻煩，巖體參數的取值直接影響到計算結果的正確性。目前，非連續變形分析的數值計算中，一般假定巖體是彈性體，對塑性、粘性體等尚不適用。c（2）計算時步的大小肯定會影響非連續變形分析的最終結果，時步過大將導致計算結果與實際測量值相差較大，從更深理論上研究時步的合理取值對非連續變形分析的計算具有較大的指導意義。計算時步的大小與實際問題聯系起來，合理地考慮時間因素也是今后需要研究的問題。（3）非連續

24、變形分析方法比較適合模擬硬巖，對于軟巖或軟硬相間的情況在判定塊體接觸時會遇到困難。也就是說非連續變形分析比較適合于塊體自身小變形的情況。（4）在非連續變形分析實施的過程中，有大量的數據需要處理，需要大量的計算機時間來存貯、轉換和調用，耗用了大量的計算機CPU時間，應優化和改進計算方法，縮短計算時間。3、DDA今后的發展方向今后的發展方向 DDA是近10年來迅速發展起來的一種較新的巖土數值計算方法，已于1996年在美國召開過一次國際性研討大會。之后，在日本、英國相繼舉行了研討會。若能將非連續變形分析和已發展成熟的有限元、邊界元以及極限平衡等方法相結合，發揮各自的優勢，則會使這一方法更趨于完善。對

25、非連續變形分析方法的研究，應重點改進前面提及的不足之處。除上述工作外，正在進行和將要進行的工作還應包括如下的內容：（1）繼續完善三維非連續變形分析理論以及計算程序的實現仍然是今后一個時期內需要解決的問題，優化計算機程序和設計方法是提高解題效率和推廣非連續變形分析法應用的前提。（2）從理論和計算機程序設計上解決顆粒形狀單元在不同顆粒級配條件下的非連續變形分析。對常見的幾種概率分布形式下的顆粒級配沿途進行非連續變形分析法的模擬計算。（3）滲流問題是影響巖體穩定的重要因素之一，因此有必要在非連續變形分析理論和相應的程序中加入解決滲流問題的功能，加強滲流問題與非連續變形分析法的結合具有重要的意義。（4

26、）加強流行、分形、混沌等新興學科在非連續變形分析中的應用與推廣，進一步研究圖形理論與非連續變形分析法的結合。DDA方法作為一種新型的巖土數值計算方法，有著廣闊的應用前景。鄔愛箐將DDA應用于三峽高邊坡，梁國平提出的LDDA也是當前研究的熱門課題，LDDA（Lagrange Discontinuous Deformation Analysis）方法即具有拉格朗日乘子的不連續變形方法，是模擬非線性界面的有效方法，其在界面上設置接觸單元來模擬張裂和滑動。拉格朗日乘子即是接觸面上的正壓力和切向力。由于區域分解方法的使用，使求解的效率大大提高；由于使用了DDA方法中接觸判斷準則，因此不必定義滑移面。 無

27、單元法是一種新的數值計算方法，其特點是采用滑動最小二乘法所產生的光滑函數來近似場函數，計算形函數，從而只需計算域的幾何邊界及計算點，擺脫了單元限制，大大簡化了前處理工作。由于提供了場函數的連續可導近似解，在材料分析中，使得位移及應力、應變計算結果均連續，不需進行后處理修勻。無單元法的結點生成非常容易，根本不存在其它方法所遇到的網格重構問題，因此在開裂計算中將有很好的應用前景。總之，無單元法保留了有限元的一些特點，克服了有限元的不足，適于進行巖土工程數值模擬，尤其便于跟蹤裂紋擴展，提供了巖土工程數值模擬的新途徑。無單元法的最早思想由Nayroles等于1992年提出，并被稱之為虛擬單元法(Dif

28、fuse Element Method)。但Nayroles的近似比較簡單，使用了低階積分，邊界條件的引入也不準確，而且在擬合函數的求導中忽略了一項，計算比較粗糙。Belytschko等對Nayroles的方法進行了改進，提出了無單元伽遼金法(Element-free Galerkin Method)，采用了高階高斯積分，引入了被Nayroles忽略的項，并用拉格朗日乘子引入本征邊界條件。Lu等又將無單元伽遼金法作了進一步改進，在滑動最小二乘法中采用了正交基，大大加快了計算效率和精度。周維垣等4用罰函數法引入了本征邊界條件，并對一些關鍵問題作了探討和改進，使得無單元法得到了與有限元法類似的方程

29、組形式，并提高了計算精度。 由于無單元法具有前后處理簡單、計算精度高、在計算過程當中便于增減結點等優點，已被廣泛應用于板彎曲計算，彈塑性問題分析以及線彈性開裂分析中。1、基本原理基本原理 滑動最小二乘法由Lancaster提出，但直到近年才由Nayroles等將其應用于求解邊值問題。它是用加權最小二乘法來近似場函數的一種方法，其詳細原理參考文獻。為了克服滑動最小二乘法在系數矩陣形態較差時計算精度低的缺點，在實際計算中采用了將基函數正交化的方法。由變分原理即：平衡位移使系統總勢能取駐值，并采用罰函數法引入位移邊界條件，在位移邊界上作用剛性很大的彈簧，并將已知位移作為彈簧的初始拉伸位移。無單元法計

30、算流程如圖N.7。開始生成節點并計算節點影響半徑存儲各網格所包含的節點號確定高斯點位置及其積分權檢索高斯點的影響節點計算高斯點對整體平衡方程的貢獻并將其集入整體平衡方程解平衡方程求應力、應變等所需變量并進行結果分析結束高斯點上逐個計算圖N.7 無單元法的計算流程2、方程組解法方程組解法 由無單元法平衡方程的集成過程可以看出，其剛度矩陣具有對稱正定性和稀疏性。如果節點號編制合理，無單元法還可以做到剛度矩陣非零元素的帶形分布。但有兩個問題：（1）帶寬較大因為無單元法中形函數是由最小二乘法得到，在求每個高斯點上的函數值時，為了使最小二乘法的正規方程系數矩陣非奇異，通常應使每個高斯點有較多的影響節點，

31、而不象有限元法中只有位于同一單元中的節點才相互作用。這使得無單元法即使節點號編制合理，其帶寬一般也是較大的。（2）不能發揮無單元法的優越性 無單元法最大的特點是消除了單元，而只用節點進行計算，這使得它在某些方面的應用有很大的優越性。首先，它可以任意加密節點，而不需象有限元那樣重新劃分網格；其次，在追蹤斷裂擴展時很方便，可在裂尖不斷布點，而不需重新劃分網格。但在這種情況下，要做到帶寬優化就很困難了，除非每加一個節點或每移動一個點都重新進行一次節點編號的優化，但如果實際操作起來，其計算量將是非常大的。為了克服這些缺陷，無單元法方程組的求解采用了一種迭代解法共軛梯度法，共軛梯度法是用來求解系數矩陣為

32、對稱正定稀疏矩陣的大型線性方程組的一種有效方法，其基本理論及算法參見文獻。 界面元又稱之為剛體彈簧元、剛體元、界面應力元等，它是基于1980年日本東京大學Kaiwai教授和周維垣、楊若瓊提出的剛體彈簧元模型而建立有限元計算的。它綜合了目前巖土工程數值分析中連續介質變形體模型和不連續剛體模型的特點，其建模思路可概括為三個方面：（1）假定單元的變形累積在界面層，單元本身就可視為剛體，因而可采用分片剛體位移模式構造結構的應力場，使其逼近真實的位移場。（2）基于相鄰塊體單元的相對變形和不同的材料特性（非均質、各向異性、彈塑性、非線性以及粘性），導出界面應力公式并以此表征結構的應力場。（3）根據微分方程的弱解形式加權殘量法或虛功原理或者放松界面位移連續條件的廣義變分原理，建立相應的支配方程。1、界