FLAC3D原理..

1、2.2三維數(shù)值模擬方法及其原理2.2.1 FLAC3D工程分析軟件特點(diǎn)FLAC3D是由美國ItascaConsultingGroup,Inc.為地質(zhì)工程應(yīng)用而開發(fā)的連續(xù)介質(zhì)顯式有限差分計(jì)算機(jī)軟件。FLAC即FastLagrangianAnalysisofContinua的縮寫。該軟件主要適用于模擬計(jì)算巖土體材料的力學(xué)行為及巖土材料達(dá)到屈服極限后產(chǎn)生的塑性流動(dòng)，對(duì)大變形情況應(yīng)用效果更好。FLAC3D程序在數(shù)學(xué)上采用的是快速拉格朗日方法,基于顯式差分來獲得模型全部運(yùn)動(dòng)方程和本構(gòu)方程的步長解，其本構(gòu)方程由基本應(yīng)力應(yīng)變定義及虎克定律導(dǎo)出，運(yùn)動(dòng)平衡方程則直接應(yīng)用了柯西運(yùn)動(dòng)方程，該方程由牛頓運(yùn)動(dòng)定律導(dǎo)出。

2、計(jì)算模型一般是由若干不同形狀的三維單元體組成，也即剖分的空間單元網(wǎng)絡(luò)區(qū)，計(jì)算中又將每個(gè)單元體進(jìn)一步劃分成由四個(gè)節(jié)點(diǎn)構(gòu)成的四面體，四面體的應(yīng)力應(yīng)變只通過四個(gè)節(jié)點(diǎn)向其它四面體傳遞，進(jìn)而傳遞到其它單元體。當(dāng)對(duì)某一節(jié)點(diǎn)施加荷載后，在某一個(gè)微小的時(shí)間段內(nèi)，作用于該點(diǎn)的荷載只對(duì)周圍的若干節(jié)點(diǎn)(相鄰節(jié)點(diǎn))有影響。利用運(yùn)動(dòng)方程，根據(jù)單元節(jié)點(diǎn)的速度變化和時(shí)間，可計(jì)算出單元之間的相對(duì)位移，進(jìn)而求出單元應(yīng)變，再利用單元模型的本構(gòu)方程，可求出單元應(yīng)力。在計(jì)算應(yīng)變過程中，利用高斯積分理論，將三維問題轉(zhuǎn)化為二維問題而使其簡單化。在運(yùn)動(dòng)方程中，還充分考慮了巖土體所具有的粘滯性，將其視作阻尼附加于方程中。FLAC3D具有一個(gè)

3、功能強(qiáng)大的網(wǎng)格生成器，有12種基本形狀的單元體可供選擇，利用這12種基本單元體，幾乎可以構(gòu)成任何形狀的空間立體模型。FLAC3D主要是為地質(zhì)工程應(yīng)用而開發(fā)的巖土體力學(xué)數(shù)值評(píng)價(jià)計(jì)算程序，自身設(shè)計(jì)有九種材料本構(gòu)模型：(1)空模型(NullModel)(2)彈性各向同性材料模型(Elastic,IsotropicModel)(3)彈性各向異性材料模型(Elastic,anisotropicModel)(4)彳惠拉克-普拉格彈塑性材料模型(Drucker-PragerModel)(5)莫爾-庫倫彈塑性材料模型(Mohr-CoulombModel)(6)應(yīng)變硬化、軟化彈塑性材料模型(Strain-Har

4、dening/SofteningMohr-CoulombModel)(7)多節(jié)理裂隙材料模型(Ubiquitous-JointModel)(8)雙曲型應(yīng)變硬化、軟化多節(jié)理裂隙材料模型(BilinearStrain-Hardening/SofteningUbiquitous-JointModel)(9)修正的Cam粘土材料模型(ModifiedCam-clayModel)除上述本構(gòu)模型之外，F(xiàn)LAC3D還可進(jìn)行動(dòng)力學(xué)問題、水力學(xué)問題、熱力學(xué)問題等的數(shù)值模擬。在邊界條件及初始條件的考慮上，F(xiàn)LAC3D軟件十分靈活方便，可在數(shù)值計(jì)算過程中隨時(shí)調(diào)整邊界條件和初始條件。FLAC3D具有強(qiáng)大的后處理功能，

5、用戶可以直接在屏幕上繪制或以文件形式創(chuàng)建或輸出打印多種形式的圖形、文字，用戶還可根據(jù)各自的需要，將若干個(gè)變量合并在同一幅圖形中進(jìn)行研究分析。FLAC3D軟件還可對(duì)各種開挖工程或施加支護(hù)工程等進(jìn)行數(shù)值仿真模擬，軟件自身設(shè)計(jì)有錨桿、錨索、襯砌、支架等結(jié)構(gòu)元素，可以直接模擬這些支護(hù)于圍巖(土)體的相互作用。FLAC3D擁有可以自行設(shè)計(jì)的FISH語言，用戶可根據(jù)自身需求，自己設(shè)計(jì)材料的本構(gòu)模型、屈服準(zhǔn)則、支護(hù)方案、復(fù)雜形狀的開挖方式等工作。特別注意的是，巖石是一種脆性材料，當(dāng)外荷載達(dá)到巖石強(qiáng)度后，材料發(fā)生斷裂破壞，產(chǎn)生弱化現(xiàn)象，應(yīng)屬于彈塑性體。在FLAC3D中，一般對(duì)于彈塑性材料，判斷其破壞與否的基本

6、準(zhǔn)則有兩個(gè)，即Drucker-Prager準(zhǔn)則和Mohr-Coulomb準(zhǔn)則。根據(jù)室內(nèi)巖石力學(xué)性質(zhì)試驗(yàn)結(jié)果，其典型應(yīng)力應(yīng)變曲線反映出巖體破壞包絡(luò)線符合莫爾一庫倫屈服準(zhǔn)則，故本次建立的本構(gòu)力學(xué)模型選擇莫爾一庫倫彈塑性材料模型為宜。2.2.2FLAC3D分析計(jì)算原理計(jì)算所采用的數(shù)學(xué)模型是根據(jù)彈塑性理論的基本原理(應(yīng)變定義、運(yùn)動(dòng)定律、能量守衡定律、平衡方程及理想材料的連續(xù)性方程等)而建立的。2.2.2.1 基本約定在數(shù)學(xué)及數(shù)值模型的表達(dá)式中，符號(hào)有一定的約定含義，一般A表示張量，Aij表小張量A的（i,j）分量，a表小矢量，ai表小矢量a的i分量,表小對(duì)Xi的偏導(dǎo)數(shù)。Xi,Ui,vi和dvi/dt,

7、（i=1,3）分別表示一點(diǎn)的位置矢量分量、位移矢量分量、速度矢量分量和加速度矢量分量。2.2,2.2數(shù)學(xué)模型（一）柯西（Cauchy）應(yīng)力張量與柯西公式對(duì)于一個(gè)具有體積V的封閉曲面s的物體，在其上取一表面元素s,這個(gè)表面元素的單位外法向矢量為n,在某一時(shí)刻t,在表面元素對(duì)于連續(xù)介質(zhì)中一點(diǎn)，作用著對(duì)稱的應(yīng)力張量ij,根據(jù)s上作用有力巳則極限1sm0dPds稱為表面力。若用ti表示T的分量，則在三維直角坐標(biāo)系中可有關(guān)系式tiniij(1)這個(gè)關(guān)系式稱為柯西公式，其中,（二）應(yīng)變速率和旋轉(zhuǎn)速率ij稱為柯西應(yīng)力張量。如果介質(zhì)質(zhì)點(diǎn)具有運(yùn)動(dòng)速度矢量v,則在一個(gè)無限小的時(shí)間dt內(nèi)，介質(zhì)會(huì)產(chǎn)生一個(gè)由vidt決

8、定的無限小應(yīng)變，對(duì)應(yīng)的應(yīng)變速率分量ij為ij1產(chǎn)X）2XjXi而其旋轉(zhuǎn)速率分量ij為ijJ/上）2Xjx(2)(3)（三）運(yùn)動(dòng)及平衡方程根據(jù)牛頓運(yùn)動(dòng)定律與柯西應(yīng)力原理，如果質(zhì)點(diǎn)作用著應(yīng)力具有速度vi,則在無限小時(shí)間段dt內(nèi)，它們之間的關(guān)系為ij,dvibiXjdtij與體力bi,且（4）式中，為質(zhì)點(diǎn)密度。（4）式稱為柯西運(yùn)動(dòng)方程。當(dāng)質(zhì)點(diǎn)的加速度為零時(shí)，上式變?yōu)殪o力平衡方程bi0（5）Xj（四）本構(gòu)方程上述（4）式與（5）式組成的方程組中含有9個(gè)方程，15個(gè)未知量，其中12個(gè)是應(yīng)力與應(yīng)變速率分量，3個(gè)是速度分量。其余6個(gè)關(guān)系式則由本構(gòu)方程提供，本構(gòu)方程一般具有如下形式ijHj（0,0,）（6）式

9、中，為應(yīng)力變化速率，H表示一個(gè)特定的函數(shù)關(guān)系，為與荷載歷史有關(guān)的參數(shù)。2.2.2.3數(shù)值模型flac3D的數(shù)值剖分網(wǎng)格在計(jì)算中是按照四面體進(jìn)行的，四面體的節(jié)點(diǎn)也既是網(wǎng)格剖分的節(jié)點(diǎn)，因此，每個(gè)計(jì)算單元有4個(gè)面和4個(gè)節(jié)點(diǎn)（見圖2-1）圖2-1四面體單元示意圖（一）空間微分的有限差分逼近對(duì)于一個(gè)計(jì)算單元，若內(nèi)部各質(zhì)點(diǎn)速度為一連續(xù)的矢量場（GausS）積分原理有VidvvinjdsVXjSjv,則根據(jù)高斯式中，n為外法向單位矢量場。由于單元體的應(yīng)變速率是連續(xù)的，因此可以近似認(rèn)為速度是線性變化的,則（7）式可用下面的求和公式近似逼近V*4Vixjf1(f)n(f)q(f)njS(8)式中，v為單元體體積

10、,f=1,4,Vi為面平均速度的(9)Vixj13i4lVi1f4(f)nj1,flS(f)(10)S為單元體某一面的面積，f為單元體面數(shù),分量。由于速度場是線性的，則有-(f)1iVi-Vi311,if式中，l為單元體節(jié)點(diǎn)數(shù)，l=1,4;vi1為l節(jié)點(diǎn)的i速度分量將（9）式代入（8）式可得(11)則將（10）式兩邊除以Vixj13ViV,并將(4vln(l)S111）式代入可得(12)根據(jù)正交原理有4（f）Q（f）njSf1因此有146Vi1/i(l)i(VinjVjZ)ni)S(13)（二）運(yùn)動(dòng)方程的節(jié)點(diǎn)公式根據(jù)前面對(duì)質(zhì)點(diǎn)運(yùn)動(dòng)方程的討論，對(duì)于連續(xù)介質(zhì)，當(dāng)處于平衡狀態(tài)時(shí)，其平衡方程為一LBi

11、0（14）xjBi（bi*（15）dt介質(zhì)可由若干個(gè)作用著體力B,各自產(chǎn)生一定變形的四面體組成，設(shè)節(jié)點(diǎn)力為fn,n=1,4,利用虛功原理，假定單元體節(jié)點(diǎn)具有速度vn,內(nèi)部具有變形速率,則由節(jié)點(diǎn)力fn和體力B所做的外功率與ij所做的內(nèi)功率應(yīng)當(dāng)相外功率E可用下式表小EvnfinVvBdV(16)n1而內(nèi)功率I為IVjjdV(17)由(13)式，對(duì)于恒定應(yīng)變速率單元體可有1 4I1(vl"vjn)Sl(18)6 i1由于應(yīng)力張量是對(duì)稱的，定義矢量TlTi1ns1(19)可得14I-vlTil(20)3 l1將(15)式代入(16)式得4EvfinEbEI(21)n1式中，Eb和E1分別是體

12、力bi和慣性力產(chǎn)生的外功率，對(duì)于單元體內(nèi)恒定的體力bi,Eb可寫成EbhVvdV(22)而E1可以寫成EIvi皿dV(23)Vdt如前所述，單元體內(nèi)部速度場以線性變化，為方便描述，選取單元體質(zhì)心為原點(diǎn)，x、x2、*3為坐標(biāo)軸，可有4vv：Nn(24)n1(25)(26)式中，Nn(n=1,4)是具有如下形式的線性函數(shù)nnnnnNC0CiXiC2X2C3X3Cnnncnc0、c1、c2、c3(n=1,4)則由下面方程決止NnX1j,X2j,X3jnjnj為克羅內(nèi)克爾(Kronecker)記號(hào)，根據(jù)質(zhì)心定義可得用克萊默Ebbivinc0Vn1(Cramer's)法則解得EbEI在穩(wěn)定狀態(tài)下

13、,nC04,則nbiVVi，"4nVin1nVifinNn>biVTVn"Mv內(nèi)功率與外功率必定相等，則有fnTinbiVi34上式最后一項(xiàng)，如果單元體內(nèi)ndViNnLdV-VdtndViNnidVVdt恒定，則根據(jù)質(zhì)心定義有V(dVi)n4dt若將V/4看作是假定的節(jié)點(diǎn)質(zhì)量mn,則上式變?yōu)閚dVinLdVdVinfinbiV4n,dvi、nmE根據(jù)牛頓定律有l(wèi)1,nn式中，F(xiàn)i<l>為節(jié)點(diǎn)l所受力的分量，M<l>為節(jié)點(diǎn)l的質(zhì)量,(28)(29)(30)(31)(32)(33)(34)(35)(也)l為節(jié)點(diǎn)dtl在Fi<l>作用下產(chǎn)

14、生白加速度，nn為包含在全部連續(xù)介質(zhì)中單元節(jié)點(diǎn)總數(shù)。設(shè)質(zhì)點(diǎn)質(zhì)量Mlml,不平衡力F<>為lFil6等fil(36)34(三)時(shí)間微分的顯式有限差分對(duì)于上面(35)式，可寫成(37)dvil1l()Fil1,nndtMlFLAC3D計(jì)算中，對(duì)于上式的逼近是在一個(gè)足夠小t時(shí)間段內(nèi)，采用中心(38)差分格式，即上式變?yōu)閂i(t)Vi(t)-Fi22Ml同理，節(jié)點(diǎn)的坐標(biāo)差分公式與位移差分公式分別為Xil(tt)Uil(tt)iltXi(t)tvi(t)2uJ(t)tvj(t')(39)(40)(四)本構(gòu)方程的增量形式FLAC3D中，假定在t表小為時(shí)間段內(nèi)，速度是不變的，本構(gòu)方程的增

15、量形式可*ijHij(41)jtj則ij為t時(shí)間段內(nèi)應(yīng)變?cè)隽?。?yīng)力增量由下式確定(42)式中,ijijiC由下式定義(43)(ikkjikkj)t(44)14l6nVinjllolvjniS(45)2.2.2.4本構(gòu)模型(一)彈塑性理論的增量關(guān)系(46)一般來講，彈塑性體破壞的判斷準(zhǔn)則具有如下形式f(-i)0式中，f為一特定的函數(shù)關(guān)系，一為n維應(yīng)力矢量i分量。彈塑性體的變形為彈性變形與塑性變形p的總和，即ep_i_i_i彈性應(yīng)變與應(yīng)力之間的關(guān)系為_e_iSi_ni1,n其中，S為線性函數(shù)。塑性變形可用下式表示p_g(47)(48)(49)式中，為一常數(shù)，g為（_n）的某一函數(shù)關(guān)系。將（47）、

16、（49）式代入（48）式得g_iSi_nSi（-）（50）一n當(dāng)f（_n）為線性函數(shù)時(shí),（51）式可寫成_*_f_nf_n0其中，f.f.f0no由于f（_n）0,再將（50）式代入（52）式可得_*_*_gfSn_nf?()0-n若記則有因此有N_i_i_i1s_i_iSi_n_If_nfSn(_n)(51)(52)(53)(54)InfSng/_nf0n對(duì)于新應(yīng)力_n_n,仍有則新應(yīng)力為(二)莫爾庫倫(Mohr-Coulomb)模型本構(gòu)關(guān)系Hooke'slaw),應(yīng)力應(yīng)變的關(guān)系為(57)本次計(jì)算的本構(gòu)模型選擇的是莫爾一庫倫彈塑性材料模型，彈塑性體可產(chǎn)生彈性及塑性兩部分變形，根據(jù)虎克

17、定律(E其中，E為剛度矩陣。對(duì)于彈性變形，應(yīng)力增量可由下式確定式中,e1e2e3e2e1e1e3e3e2(58)2可由剪切模量與體積模量得出4G32G3(58)式還可寫成S2S3e1,e1,e1,e2,e2,e2,e3e3e3e1e2e3e2e1e1e3e3e2(59)(三)破壞準(zhǔn)則與流動(dòng)法則若123,則莫爾一庫倫破壞準(zhǔn)則具有形式fs13Nfs=0,2cN其中(60)對(duì)于拉張破壞，有ftft=0的形式，其中t(61)上兩式中，為內(nèi)摩擦角，C為內(nèi)聚力,N1sin()t為抗拉強(qiáng)度,1sin()對(duì)于材料的抗拉強(qiáng)度:其值不可能超過3(見圖2-2),(62)t的最大值為tmaxctg()對(duì)于塑性流動(dòng)狀態(tài)

18、，設(shè)g,與gt分別為剪切破壞和拉張破壞所對(duì)應(yīng)的函數(shù)關(guān)系,有g(shù)si3Ntg3式中，為膨脹角，(64)(65)(66)N流動(dòng)準(zhǔn)則可定義為函數(shù)關(guān)系h（1，3）0,在1、3平面內(nèi)將fs0和圖2-3莫爾一庫倫模型塑性流分區(qū)ft0以上部分分為兩個(gè)區(qū)域（見圖2-3）,其函數(shù)表達(dá)式如下.tPPh3a1其中，ap與p由下式定義N2NPtN2cN(68)從圖2-3中可見，若某點(diǎn)的應(yīng)力狀態(tài)在1區(qū)（fs0,h0）時(shí)，則發(fā)生剪切破壞，在2區(qū)（fs0,h0）時(shí)，則發(fā)生拉張破壞。（四）塑性階段應(yīng)力修正對(duì)于剪切破壞，由sg(64)式可得(69)3將（59）式中的gs/1、gs/2和gs/3代Si2NS2(70)1N再由(55)、(56)N1N2N