1、1 三維快速拉格朗日法的基本原理1.1 概述目前在巖土力學(xué)中常用的數(shù)值計(jì)算方法有差分方法、有限元法、邊界元法等幾種，特別是后兩種方法，隨著計(jì)算機(jī)的發(fā)展其應(yīng)用尤為廣泛。但是，這幾種方法都是以連續(xù)介質(zhì)為出發(fā)點(diǎn)，而且往往囿于小變形的假定。它們雖然也可以用來(lái)解決由幾種介質(zhì)所組成的非均質(zhì)的問(wèn)題，并且對(duì)于個(gè)別的斷層或弱面，也可以用設(shè)置節(jié)理單元的辦法來(lái)解決，但是用以解決富含節(jié)理和大變形的巖土力學(xué)問(wèn)題，往往所得的結(jié)果與實(shí)際的物理圖景相差甚遠(yuǎn)。于是離散單元法和拉格朗日元法就應(yīng)運(yùn)而生。離散單元法是Cundall于上世紀(jì)70年代初所提出的。該法將為弱面所切割的巖體視為復(fù)雜的塊體的集合體，允許各個(gè)塊體可以平移或轉(zhuǎn)動(dòng)，

2、甚至相互分離。拉格朗日元法則是由Cundall所加盟的美國(guó)ITASCA咨詢集團(tuán)于1986年所開(kāi)發(fā)的。該法將流體力學(xué)中跟蹤流體運(yùn)動(dòng)的拉格朗日方法應(yīng)用于解決巖體力學(xué)的問(wèn)題獲得成功。三維快速拉格朗日法是一種基于三維顯式有限差分法的數(shù)值分析方法，它可以模擬巖土或其他材料的三維力學(xué)行為。三維快速拉格朗日分析將計(jì)算區(qū)域劃分為若干四面體單元，每個(gè)單元在給定的邊界條件下遵循指定的線性或非線性本構(gòu)關(guān)系，如果單元應(yīng)力使得材料屈服或產(chǎn)生塑性流動(dòng)，則單元網(wǎng)格可以隨著材料的變形而變形，這就是所謂的拉格朗日算法，這種算法非常適合于模擬大變形問(wèn)題。三維快速拉格朗日分析采用了顯式有限差分格式來(lái)求解場(chǎng)的控制微分方程，并應(yīng)用了混

3、合單元離散模型，可以準(zhǔn)確地模擬材料的屈服、塑性流動(dòng)、軟化直至大變形，尤其在材料的彈塑性分析、大變形分析以及模擬施工過(guò)程等領(lǐng)域有其獨(dú)到的優(yōu)點(diǎn)。1.2 三維快速拉格朗日分析的數(shù)學(xué)模型三維快速拉格朗日分析在求解中使用如下3種計(jì)算方法：(1)離散模型方法。連續(xù)介質(zhì)被離散為若干六面體單元，作用力均被集中在節(jié)點(diǎn)上。(2)有限差分方法。變量關(guān)于空間和時(shí)間的一階導(dǎo)數(shù)均用有限差分來(lái)近似。(3)動(dòng)態(tài)松馳方法。由質(zhì)點(diǎn)運(yùn)動(dòng)方程求解，通過(guò)阻尼使系統(tǒng)運(yùn)動(dòng)衰減至平衡狀態(tài)。1.2.1 空間導(dǎo)數(shù)的有限差分近似快速拉格朗日分析采用混合離散方法，將區(qū)域離散為常應(yīng)變六面體單元的集合體，又將每個(gè)六面體看作以六面體角點(diǎn)為角點(diǎn)的常應(yīng)變四面

4、體的集合體，應(yīng)力、應(yīng)變、節(jié)點(diǎn)不平衡力等變量均在四面體上進(jìn)行計(jì)算，六面體單元的應(yīng)力、應(yīng)變?nèi)≈禐槠鋬?nèi)四面體的體積加權(quán)平均。這種方法既避免了常應(yīng)變六面體單元常會(huì)遇到的位移剪切鎖死現(xiàn)象，又使得四面體單元的位移模式可以充分適應(yīng)一些本構(gòu)的要求，如不可壓縮塑性流動(dòng)等。圖1-1 四面體單元的面和節(jié)點(diǎn)如一四面體,節(jié)點(diǎn)編號(hào)為1到4，第n面表示與節(jié)點(diǎn)n相對(duì)的面，設(shè)其內(nèi)一點(diǎn)的速率分量為vi,由高斯公式得： (1-1) 其中V為四面體的體積，S為四面體的外表面，nj為外表面的單位法向向量分量。對(duì)于常應(yīng)變單元，vi為線性分布，nj在每個(gè)面上為常量。對(duì)式(1-1)積分得： （1-2）式中，上標(biāo)（f）指面f的相關(guān)變量值，指i

5、速度分量的均值。若速度呈線性變化，則： （1-3）上標(biāo)l指節(jié)點(diǎn)l的值。將上式代入式（1-2），有： （1-4）在式（1-1）中，若vi=1，應(yīng)用高斯法則可得： （1-5）所以，式（1-4）兩邊同除以V，則有： （1-6）而應(yīng)變速率張量則可由下式表示：應(yīng)變速率張量的分量形式為： (1-7)1.2.2 節(jié)點(diǎn)運(yùn)動(dòng)方程一定時(shí)域內(nèi)，靜力平衡問(wèn)題可通過(guò)以下的平衡方程求解得到： (1-8)式中：為介質(zhì)密度，bi為介質(zhì)單位質(zhì)量的體積力。根據(jù)虛功原理，作用于單個(gè)四面體上的節(jié)點(diǎn)力fl（l=（1,4）與四面體應(yīng)力和等效體力相平衡。引入節(jié)點(diǎn)虛速度vl（它在四面體中產(chǎn)生線性速度場(chǎng)v和常應(yīng)變速率），則節(jié)點(diǎn)力Fl和體力B產(chǎn)

6、生的外力功功率等于內(nèi)部應(yīng)力ij產(chǎn)生的內(nèi)力功功率。外力功功率可表示為： （1-9）而內(nèi)力功功率： （1-10）由式（1-7），對(duì)常應(yīng)變速率的四面體有： （1-11）應(yīng)力張量是對(duì)稱張量，定義矢量Tl： （1-12）則： （1-13）式（1-8）代入式（1-9），有： （1-14）Eb和EI分別為體力和慣性力所作的外力功功率。若四面體內(nèi)體力為常數(shù)，則有： （1-15） （1-16）根據(jù)有限差分近似，速度場(chǎng)在四面體內(nèi)線性變化。為描述它，引進(jìn)一個(gè)參考坐標(biāo)系（它的坐標(biāo)原點(diǎn)則四面體的中心上），則有： （1-17）式中Nn（n=1，4）為一線性函數(shù)： （1-18）其中，（n=1，4）為下述方程的解： （1-1

7、9）式中，是克羅內(nèi)克爾增量（Kronecker delta）。通過(guò)中心點(diǎn)的定義，所有形如的積分均為0，將式（1-18）、式（1-17）代入式（1-14）得： （1-20）由克雷姆定律，解式(1-19)得： （1-21）將上式代入式（1-20），有： （1-22）同理，將式（1-17）代入式（1-16）得到： （1-23）將式（1-22）和（1-23）代入式（1-14）： （1-24）對(duì)任何虛速度，外虛功率E等于內(nèi)虛功率I： （1-25）在四面體范圍內(nèi)，加速度場(chǎng)空間變化是微小的，則有： （1-26）為不變量，則上式可寫為： （1-27）用假想的節(jié)點(diǎn)質(zhì)量mn代替上式中的質(zhì)量：則，式（1-25）可寫

8、為： （1-28）對(duì)于等效體系，可以建立平衡狀態(tài)，要求在每個(gè)節(jié)點(diǎn)上靜態(tài)等效載荷之和為零。可以寫出全部節(jié)點(diǎn)上牛頓定律表達(dá)式： （1-29）式中，nn介質(zhì)中的所有的節(jié)點(diǎn)總數(shù)，節(jié)點(diǎn)質(zhì)量定義為： （1-30）不平衡力F<l>定義為： （1-31）當(dāng)介質(zhì)達(dá)到平衡時(shí)，不平衡力等于0。1.2.3 增量形式的本構(gòu)方程快速拉格朗日分析中，假定時(shí)間內(nèi)速度為常數(shù)，增量形式的本構(gòu)方程可表示為： (1-32)式中，稱為共轉(zhuǎn)（co-rotational）應(yīng)力增量，為一給定的函數(shù)。共轉(zhuǎn)（co-rotational）應(yīng)力速率張量等于給定參考系的介質(zhì)內(nèi)一點(diǎn)應(yīng)力的偏導(dǎo)數(shù)和以瞬時(shí)角速度的轉(zhuǎn)動(dòng)，數(shù)學(xué)表達(dá)式為： (1-33

9、)式中，w為轉(zhuǎn)動(dòng)速率張量。利用有限差分方程，可以得到轉(zhuǎn)動(dòng)速率張量的分量形式： (1-34)式中符合同前。1.2.4 時(shí)間導(dǎo)數(shù)的有限差分近似由本構(gòu)方程（式（1-32）和變形速率與節(jié)點(diǎn)速率之間的關(guān)系（式（1-7），式（1-26）可表示為一般的差分方程： （1-35）式中，<l>是指在計(jì)算過(guò)程中全局節(jié)點(diǎn)l節(jié)點(diǎn)速度值的子集（式（1-29）。在時(shí)間間隔中實(shí)際節(jié)點(diǎn)的速度假定是線性變化的，式(1-35)左邊導(dǎo)數(shù)用中心有限差分估算。 （1-36）類似地，節(jié)點(diǎn)的位置也用中心有限差分進(jìn)行迭代： （1-37）因此，節(jié)點(diǎn)位移也有如下關(guān)系： （1-38）1.2.5 阻尼力為使運(yùn)動(dòng)方程獲得靜態(tài)或準(zhǔn)靜態(tài)（非慣性

10、）解，快速拉格朗日分析的靜力分析中，在式(1-29)中加入非粘性阻尼力。則式(1-29)變?yōu)椋?(1-39)式中：為阻尼力，為阻尼系數(shù)，其默認(rèn)值為0 8。 (1-40)1.3 FLAC3D簡(jiǎn)介由以上原理可以看出，無(wú)論是動(dòng)態(tài)問(wèn)題，還是靜態(tài)問(wèn)題，三維快速拉格朗日分析均由運(yùn)動(dòng)方程用顯式方法進(jìn)行求解，這使得它很容易模擬動(dòng)態(tài)問(wèn)題，如振動(dòng)、失穩(wěn)、大變形等。對(duì)顯式法來(lái)說(shuō)非線性本構(gòu)關(guān)系與線性本構(gòu)關(guān)系并無(wú)算法上的差別，對(duì)于已知的應(yīng)變?cè)隽浚珊芊奖愕厍蟪鰬?yīng)力增量，并得到不平衡力，就同實(shí)際中的物理過(guò)程一樣，可以跟蹤系統(tǒng)的演化過(guò)程。此外，顯式法不形成剛度矩陣，每一步計(jì)算所需計(jì)算機(jī)內(nèi)存很小，使用較少的計(jì)算機(jī)內(nèi)存就可以模

11、擬大量的單元，特別適于在微機(jī)上操作。在求解大變形過(guò)程中，因每一時(shí)步變形很小，可采用小變形本構(gòu)關(guān)系，只需將各時(shí)步的變形疊加，即得到了大變形。這就避免了通常大變形問(wèn)題中推導(dǎo)大變形本構(gòu)關(guān)系及其應(yīng)用中所遇到的麻煩，也使它的求解過(guò)程與小變形問(wèn)題一樣。根據(jù)前述原理，美國(guó)Itasca Consulting Group開(kāi)發(fā)了三維快速拉格朗日分析程序FLAC一3D，該程序能較好地模擬地質(zhì)材料在達(dá)到強(qiáng)度極限或屈服極限時(shí)發(fā)生的破壞或塑性流動(dòng)的力學(xué)行為，特別適用于分析漸進(jìn)破壞和失穩(wěn)以及模擬大變形。它主要有如下一些特點(diǎn):(l)應(yīng)用范圍廣泛，可以模擬復(fù)雜的巖土工程或力學(xué)問(wèn)題。FLAC 3D包含了10種彈塑性材料本構(gòu)模型，

12、有靜力、動(dòng)力、蠕變、滲流、溫度五種計(jì)算模式，各種模式間可以互相藕合，以模擬各種復(fù)雜的工程力學(xué)行為。FLAC - 3D可以模擬多種結(jié)構(gòu)形式，如巖體、土體或其他材料實(shí)體，梁、錨元、樁、殼以及人工結(jié)構(gòu)如支護(hù)、襯砌、錨索、巖栓、土工織物、摩擦樁、板樁等，另外，F(xiàn)LAC3D設(shè)有界面單元，可以模擬節(jié)理、斷層或虛擬的物理邊界等；(2)FLAC3D具有強(qiáng)大的內(nèi)嵌程序語(yǔ)言FISH，使得用戶可以定義新的變量或函數(shù)，以適應(yīng)用戶的特殊需要。例如，利用FISH，用戶自己設(shè)計(jì)FLA C 3D內(nèi)部沒(méi)有的特殊單元形態(tài)；用戶可以在數(shù)值試驗(yàn)中進(jìn)行伺服控制；可以指定特殊的邊界條件，自動(dòng)進(jìn)行參數(shù)分析；可以獲得計(jì)算過(guò)程中節(jié)點(diǎn)、單元參數(shù)

13、，如坐標(biāo)、位移、速度、材料參數(shù)、應(yīng)力、應(yīng)奪、不平衡力等；(3)FLAC3D具有強(qiáng)大的前后處理功能。FLAC3D具有強(qiáng)大的自動(dòng)三維網(wǎng)格生成器，內(nèi)部定義了多種基本單元形態(tài)，可以生成非常復(fù)雜的三維網(wǎng)格。在計(jì)算過(guò)程中用戶可以用高分辨率的彩色或灰度圖或數(shù)據(jù)文件輸出結(jié)果，以對(duì)結(jié)果進(jìn)行實(shí)時(shí)分析，圖形可以表示網(wǎng)格、結(jié)構(gòu)以及有關(guān)變量的等值線圖、矢量圖、曲線圈等，可以給出計(jì)算域的任意截面上的變量等值線圖和矢量圖。FLAC3D具有如下缺陷：(1)對(duì)于線性問(wèn)題，F(xiàn)LAC3D要比相應(yīng)的有限元花費(fèi)更多的計(jì)算時(shí)間，F(xiàn)LAC3D在模擬非線性問(wèn)題、大變形問(wèn)題或動(dòng)態(tài)問(wèn)題時(shí)更有效。(2)FLAC3D的收斂速度取決于系統(tǒng)的最大固有周

14、期與最小固有周期的比值，這使得它對(duì)某些問(wèn)題的模擬效率非常低，如單元尺寸或材料彈性模量相差很大的情況。1.4 本構(gòu)模型FLAC 3D包含了10種彈塑性材料本構(gòu)模型，以下簡(jiǎn)單介紹報(bào)告中涉及到的4類本構(gòu)模型。1.4.1 空單元模型（Null Model）空單元材料用于描述從模型中刪除或開(kāi)挖掉的部分。在模擬的后續(xù)階段，空單元材料可以被轉(zhuǎn)變成不同的材料模型。用這種方式，例如，能夠模擬回填開(kāi)挖。在空單元區(qū)域內(nèi)的所有應(yīng)力被自動(dòng)設(shè)置成零。 (1-41)1.4.2 Mohr-Coulomb模型（Mohr-Coulomb Model）1.4.2.1 組合的破壞準(zhǔn)則和流動(dòng)法則在這個(gè)模型的破壞準(zhǔn)則是張拉剪切組合的Mo

15、hr-Coulomb準(zhǔn)則。這個(gè)準(zhǔn)則可以用圖(1-2)解釋。用Mohr-Coulomb破壞準(zhǔn)則描繪從點(diǎn)A到點(diǎn)B破壞包絡(luò)線，即： (1-42) 用式張拉破壞準(zhǔn)則描繪從B點(diǎn)到點(diǎn)C的包絡(luò)線： (1-43)式中，是摩擦角，C是粘聚力，t是張拉強(qiáng)度，且有： (1-44)張拉強(qiáng)度不超過(guò)3值。最大值由下式給定： （1-45）分別用兩個(gè)定義剪切塑性流動(dòng)和張拉塑性流動(dòng)的函數(shù)gs和gt描述勢(shì)函數(shù)。函數(shù)gs有如下形式： （1-46）圖1-2 Mohr-Coulomb 破壞準(zhǔn)則式中，是膨脹角。 （1-47）函數(shù)gt符合相關(guān)流動(dòng)法則，寫成： (1-48)用式將流動(dòng)法則寫成統(tǒng)一的形式： (1-50)式中，和是由下式定義的常

16、數(shù)： （1-51） （1-52）1.4.3 遍布節(jié)理模型（Ubiquitous-Joint Model）這個(gè)模型描述了Mohr-Coulomb模型中弱面的力學(xué)行為。弱面的方向給定，弱面的破壞準(zhǔn)則由組合的Mohr-Coulomb張拉剪切破壞包絡(luò)線組成。與Mohr-Coulomb模型一樣，弱面方向上的破壞準(zhǔn)則在剪切破壞時(shí)采用非相關(guān)聯(lián)的流動(dòng)法則，在拉張破壞時(shí)采用相關(guān)聯(lián)的流動(dòng)法則。非弱面方向的破壞準(zhǔn)則與Mohr-Coulomb模型一致，如前所述，以下簡(jiǎn)單說(shuō)明以下弱面上的破壞準(zhǔn)則。在以下各式中，各應(yīng)力分量皆是弱面局部坐標(biāo)系下的應(yīng)力分量，局部坐標(biāo)系規(guī)定如下：軸與弱面的法向一致，軸水平，軸符合右手螺旋法則。

17、FLAC-3D模型中弱面的破壞準(zhǔn)則是應(yīng)力()表示的組合的Mohr-Coulomb張拉剪切破壞準(zhǔn)則，如圖(1-3)。從點(diǎn)A到點(diǎn)B定義的包絡(luò)線：圖1-3 遍布節(jié)理模型的破壞準(zhǔn)則 (1-52)用式張拉破壞準(zhǔn)則描繪從B點(diǎn)到點(diǎn)C的包絡(luò)線： (1-53)式中，j，Cj，jt分別是弱面的摩擦角、粘聚力和張拉強(qiáng)度，對(duì)于非零摩擦角的弱面，張拉強(qiáng)度的最大值為： （1-54）勢(shì)函數(shù)由gs和gt兩個(gè)函數(shù)組成，這兩個(gè)函數(shù)分別定義剪切和張拉塑性流動(dòng)。函數(shù)gs有如下形式： （1-55）式中，j是弱面膨脹角。流動(dòng)法則用統(tǒng)一式表示為： (1-56)式中，和是由下式定義的常數(shù)： （1-57） （1-58）1.4.4 修正劍橋模型（Modified Cam-Clay Model）土的有效應(yīng)力彈塑性本構(gòu)關(guān)系，即應(yīng)力增量與應(yīng)變?cè)隽康年P(guān)系可用下式表達(dá)： = Cep (159)式中 Cep 為彈塑性矩陣。 （1-60）式中Ce為彈性矩陣；f 為屈服函數(shù)；g 為勢(shì)函數(shù)。FLAC3D中，采用相關(guān)聯(lián)流動(dòng)法則，故取g= f，A為表示硬化規(guī)律的參數(shù)。圖 1-4為修正劍橋模型