整理三維拉格朗日法計算原理

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

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

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

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

5、1四面體單元的面和節(jié)點式中，上標C）指面f的相關(guān)變量值，v.指i速度分量的均值。若速度呈線i性變化，則：1-3)v（f）=工v（l）i3il=1,l豐f上標l指節(jié)點l的值。將上式代入式（1-2），有：Vvi,.=1工vi3il=1工n(f)S(f).f=1,f知1-4)在式（1-1）中，若v=1，應(yīng)用高斯法則可得:i1-5)工n（f）S（f）=0.f=1所以，式（1-4）兩邊同除以V,則有:v-丄vinSi,.3Vi.1-6)l=1而應(yīng)變速率張量則可由下式表示應(yīng)變速率張量的分量形式為：g=-乙(vin(i)+vin(i)S(i)ij6Vijjii=1(1-7)節(jié)點運動方程一定時域內(nèi)，靜力平衡問

6、題可通過以下的平衡方程求解得到：ij,ji(1-8)(dv、式中：p為介質(zhì)密度，B=pb-,bi為介質(zhì)單位質(zhì)量的體積力。i根據(jù)虛功原理，作用于單個四面體上的節(jié)點力f(i=(1,4)與四面體應(yīng)力和等效體力相平衡。引入節(jié)點虛速度Cvi(它在四面體中產(chǎn)生線性速度場Cv和常應(yīng)變速率C則節(jié)點力Fi和體力B產(chǎn)生的外力功功率等于內(nèi)部應(yīng)力oij產(chǎn)生的內(nèi)力功功率。外力功功率可表示為：E=工8vvBdV1-9)而內(nèi)力功功率：1-10)Vijij由式(1-7)，對常應(yīng)變速率的四面體有：viQni+5viQniiijjjiji1-11)應(yīng)力張量是對稱張量，定義矢量Ti：Ti=QniS(i)1-12)則：I一1工svl

7、Tl3iil=1（1-13）式（1-8）代入式（1-9），有：E=5Vnfn+Eb+Eliin=11-14)Eb和Ei分別為體力pb和慣性力所作的外力功功率。若四面體內(nèi)體力pb為常數(shù),ii則有：Eb=pbJ5vdViVi1-15)dvEi=-Jp5vdVVidt（1-16）根據(jù)有限差分近似,速度場在四面體內(nèi)線性變化。為描述它,引進一個參考坐標系（它的坐標原點則四面體的中心上）,則有：5v=5vnNniin=1（1-17）式中Nn（n=1,4）為一線性函數(shù)：Nn=Cn+CnX+CX+CX01122331-18)其中,cn,cn,cn,cn(n=1,4)為下述方程的解：0123Nn(xj,xj,x

8、J=5123nj1-19)式中，5是克羅內(nèi)克爾增量（Kroneckerdelta）。通過中心點的定義，所有nj形如JxdV的積分均為0將式（1-18）、式（1-17）代入式（1-14）得:VjEb=pbSvnCnVii0n=11-20)由克雷姆定律，解式（1-19）得1Cn=o4將上式代入式（1-20），有:1-21)Eb=工5vn化i4n=11-22)同理，將式（1-17）代入式（1-16）得到:EI=XVnjNn4dViVdtn=11-23)將式（1-22）和（1-23）代入式（1-14）:Vnfn+in=1PLJNn巴dVi4Vdt1-24)對任何虛速度，外虛功率E等于內(nèi)虛功率I:1-2

9、5)-fn=Tn+巴TNn竺dVi34Vdt在四面體范圍內(nèi)，加速度場空間變化是微小的，則有:1-26)JpNndV=VdtP為不變量，則上式可寫為:JpNnVdvpV(dv)ndV=1idt4jdt丿1-27)用假想的節(jié)點質(zhì)量mn代替上式中的質(zhì)量PV/4:則，式（1-25）可寫為:TnpbV(dv)ni+imnJ34jdt丿fni1-28)對于等效體系，可以建立平衡狀態(tài)，要求在每個節(jié)點上靜態(tài)等效載荷之和為零?？梢詫懗鋈抗?jié)點上牛頓定律表達式：Fi=Mii(dvJ.dt丿l=1,n1-29)式中，nn介質(zhì)中的所有的節(jié)點總數(shù)，節(jié)點質(zhì)量定義為:1-30)Mi=dmlli不平衡力F定義為:F=扌+罟p

10、+P1-31)當(dāng)介質(zhì)達到平衡時，不平衡力等于0。1.2.3增量形式的本構(gòu)方程快速拉格朗日分析中，假定時間At內(nèi)速度為常數(shù)，增量形式的本構(gòu)方程可表示為:Ab=H*(c,At)ijijijij(1-32)式中，Ab稱為共轉(zhuǎn)(co-rotational)應(yīng)力增量，H*為一給定的函數(shù)。共轉(zhuǎn)ijij(co-rotational)應(yīng)力速率張量等于給定參考系的介質(zhì)內(nèi)一點應(yīng)力的偏導(dǎo)數(shù)和ij以瞬時角速度Q的轉(zhuǎn)動，數(shù)學(xué)表達式為：dbQ=jwb+bwijdtijijijkj(1-33)式中，w為轉(zhuǎn)動速率張量。利用有限差分方程，可以得到轉(zhuǎn)動速率張量的分量形式w=ij16VFCin(i)vin(i)(i)ijjii=1

11、(1-34)式中符合同前。1.2.4時間導(dǎo)數(shù)的有限差分近似由本構(gòu)方程（式（1-32）和變形速率與節(jié)點速率之間的關(guān)系（式（1-7），式（1-26）可表示為一般的差分方程：1-35)dv1i=Ft,V,v,vv-,ki=1,ndtMiiiiin式中，是指在計算過程中全局節(jié)點l節(jié)點速度值的子集（式（1-29）。在時間間隔At中實際節(jié)點的速度假定是線性變化的，式（1-35）左邊導(dǎo)數(shù)用中心有限差分估算。v-)=vt,V,VVL,/(1-36)i2i2Miiiii類似地，節(jié)點的位置也用中心有限差分進行迭代：1-37)tTOC o 1-5 h zx(t+At)=x(t)+Atv(t+)iii2因此，節(jié)點位移

12、也有如下關(guān)系：1-38)u(t+At)=u(t)+Atv(t+生)iii21.2.5阻尼力為使運動方程獲得靜態(tài)或準靜態(tài)（非慣性）解，快速拉格朗日分析的靜力分析中，在式（1-29）中加入非粘性阻尼力。則式（1-29）變?yōu)椋菏街校篿i為阻尼力，i=iisign(y)=1+1，-1,0ifififFi+二MiiiFisigni(1-39)C),idv丫iIdt丿l二1,n為阻尼系數(shù)，其默認值為08。(1-40)1.3FLAC3D簡介由以上原理可以看出，無論是動態(tài)問題，還是靜態(tài)問題，三維快速拉格朗日分析均由運動方程用顯式方法進行求解，這使得它很容易模擬動態(tài)問題，如振動、失穩(wěn)、大變形等。對顯式法來說非線

13、性本構(gòu)關(guān)系與線性本構(gòu)關(guān)系并無算法上的差別，對于已知的應(yīng)變增量，可很方便地求出應(yīng)力增量，并得到不平衡力，就同實際中的物理過程一樣，可以跟蹤系統(tǒng)的演化過程。此外，顯式法不形成剛度矩陣，每一步計算所需計算機內(nèi)存很小，使用較少的計算機內(nèi)存就可以模擬大量的單元，特別適于在微機上操作。在求解大變形過程中，因每一時步變形很小，可采用小變形本構(gòu)關(guān)系，只需將各時步的變形疊加，即得到了大變形。這就避免了通常大變形問題中推導(dǎo)大變形本構(gòu)關(guān)系及其應(yīng)用中所遇到的麻煩，也使它的求解過程與小變形問題一樣。根據(jù)前述原理，美國ItascaConsultingGroup開發(fā)了三維快速拉格朗日分析程序FLAC3D,該程序能較好地模擬

14、地質(zhì)材料在達到強度極限或屈服極限時發(fā)生的破壞或塑性流動的力學(xué)行為，特別適用于分析漸進破壞和失穩(wěn)以及模擬大變形。它主要有如下一些特點:(1)應(yīng)用范圍廣泛，可以模擬復(fù)雜的巖土工程或力學(xué)問題。FLAC3D包含了10種彈塑性材料本構(gòu)模型,有靜力、動力、蠕變、滲流、溫度五種計算模式,各種模式間可以互相藕合，以模擬各種復(fù)雜的工程力學(xué)行為。FLAC-3D可以模擬多種結(jié)構(gòu)形式，如巖體、土體或其他材料實體，梁、錨元、樁、殼以及人工結(jié)構(gòu)如支護、襯砌、錨索、巖栓、土工織物、摩擦樁、板樁等，另外，F(xiàn)LAC3D設(shè)有界面單元，可以模擬節(jié)理、斷層或虛擬的物理邊界等；FLAC3D具有強大的內(nèi)嵌程序語言FISH，使得用戶可以定

15、義新的變量或函數(shù)，以適應(yīng)用戶的特殊需要。例如，利用FISH，用戶自己設(shè)計FLAC3D內(nèi)部沒有的特殊單元形態(tài)；用戶可以在數(shù)值試驗中進行伺服控制；可以指定特殊的邊界條件，自動進行參數(shù)分析；可以獲得計算過程中節(jié)點、單元參數(shù)，如坐標、位移、速度、材料參數(shù)、應(yīng)力、應(yīng)奪、不平衡力等；FLAC3D具有強大的前后處理功能oFLAC3D具有強大的自動三維網(wǎng)格生成器，內(nèi)部定義了多種基本單元形態(tài)，可以生成非常復(fù)雜的三維網(wǎng)格。在計算過程中用戶可以用高分辨率的彩色或灰度圖或數(shù)據(jù)文件輸出結(jié)果，以對結(jié)果進行實時分析，圖形可以表示網(wǎng)格、結(jié)構(gòu)以及有關(guān)變量的等值線圖、矢量圖、曲線圈等，可以給出計算域的任意截面上的變量等值線圖和矢

16、量圖。FLAC3D具有如下缺陷：(1)對于線性問題，F(xiàn)LAC3D要比相應(yīng)的有限元花費更多的計算時間，F(xiàn)LAC3D在模擬非線性問題、大變形問題或動態(tài)問題時更有效。(2)FLAC3D的收斂速度取決于系統(tǒng)的最大固有周期與最小固有周期的比值，這使得它對某些問題的模擬效率非常低，如單元尺寸或材料彈性模量相差很大的情況。1.4本構(gòu)模型FLAC3D包含了10種彈塑性材料本構(gòu)模型，以下簡單介紹報告中涉及到的4類本構(gòu)模型?？諉卧Ｐ?NullModel)空單元材料用于描述從模型中刪除或開挖掉的部分。在模擬的后續(xù)階段，空單元材料可以被轉(zhuǎn)變成不同的材料模型。用這種方式，例如，能夠模擬回填開挖。在空單元區(qū)域內(nèi)的所有應(yīng)

17、力被自動設(shè)置成零。bN=0ij(1-41)Mohr-Coulomb模型(Mohr-CoulombModel)1.4.2.1組合的破壞準則和流動法則在這個模型的破壞準則是張拉剪切組合的Mohr-Coulomb準則。這個準則可以用圖(1-2)解釋。用Mohr-Coulomb破壞準則描繪從點A到點B破壞包絡(luò)線fs=0，即：fs=b廣3N+2c.兀(1-42)用式ft=0張拉破壞準則描繪從B點到點C的包絡(luò)線：f=bbtt3(1-43)式中，0是摩擦角，C是粘聚力，是張拉強度，且有：N1+sin|)0No1sin)o(1-44)張拉強度不超過值。最大值由下式給定：xtaij)1-45)分別用兩個定義剪切

18、塑性流動和張拉塑性流動的函數(shù)gs和gt描述勢函數(shù)。函數(shù)gs有如下形式:-aN3屮1-46)式中，/是膨脹角。N=1+血(屮)屮1一sin(V)1-47)函數(shù)gt符合相關(guān)流動法則，寫成:(1-48)用式h(a,a)二0將流動13法則寫成統(tǒng)一的形式:h=aa+ap(aap)3t1(1-50)式中，aP和aP是由下式定義的常數(shù):1-51)ap=atN2cjNj1-52)遍布節(jié)理模型(UbiqUtous-JointModel)這個模型描述了Mohr-Coulomb模型中弱面的力學(xué)行為。弱面的方向給定弱面的破壞準則由組合的Mohr-Coulomb張拉剪切破壞包絡(luò)線組成。與Mohr-Coulomb模型一樣

19、，弱面方向上的破壞準則在剪切破壞時采用非相關(guān)聯(lián)的流動法則，在拉張破壞時采用相關(guān)聯(lián)的流動法則。非弱面方向的破壞準則與Mohr-Coulomb模型一致，如前所述，以下簡單說明以下弱面上的破壞準則。在以下各式中，各應(yīng)力分量皆是弱面局部坐標系下的應(yīng)力分量，局部坐標系規(guī)定如下：x軸與弱面的法向一致，y軸水平，z軸符合右手螺旋法則。FLAC-3D模型中弱面的破壞準則是應(yīng)力（,t）表示的組合的33Mohr-Coulomb張拉剪切破壞準則，如圖（1-3）。從點A到點B定義的包絡(luò)線fsQ,t）二0:33fs=t+qtan-c（1-52）33jj從B點到點C的包絡(luò)線:=-Qtt33j(1-53)ACC1Cj/tan(pj圖1-3遍布節(jié)理模型的破壞準則用式ft二0張拉破壞準則描繪于非零摩擦角的弱面，張拉強度的最大值為：cQt=jjmaxtanj式中，匕，C號分別是弱面的摩擦角、粘聚力和張拉強度，對1-54)勢函數(shù)由gs和gt兩個函數(shù)組成，這兩個函數(shù)分別定義剪切和張拉塑性流動。函數(shù)gs有如下形式：gs=t+qta甲33j（1-55）式中，竹是弱面膨脹角。流動法則用統(tǒng)一式表示為：h=t-tpapQ-Qt)jj33j(1-