FLAC3D原理..

1、【精品文檔】如有侵權(quán)，請(qǐng)聯(lián)系網(wǎng)站刪除，僅供學(xué)習(xí)與交流FLAC3D原理.精品文檔.2.2 三維數(shù)值模擬方法及其原理2.2.1 FLAC3D工程分析軟件特點(diǎn)FLAC3D是由美國(guó)Itasca Consulting Group, Inc. 為地質(zhì)工程應(yīng)用而開(kāi)發(fā)的連續(xù)介質(zhì)顯式有限差分計(jì)算機(jī)軟件。FLAC即Fast Lagrangian Analysis of Continua 的縮寫(xiě)。該軟件主要適用于模擬計(jì)算巖土體材料的力學(xué)行為及巖土材料達(dá)到屈服極限后產(chǎn)生的塑性流動(dòng)，對(duì)大變形情況應(yīng)用效果更好。FLAC3D程序在數(shù)學(xué)上采用的是快速拉格朗日方法，基于顯式差分來(lái)獲得模型全部運(yùn)動(dòng)方程和本構(gòu)方程的步長(zhǎng)解，其本構(gòu)方

2、程由基本應(yīng)力應(yīng)變定義及虎克定律導(dǎo)出，運(yùn)動(dòng)平衡方程則直接應(yīng)用了柯西運(yùn)動(dòng)方程，該方程由牛頓運(yùn)動(dòng)定律導(dǎo)出。計(jì)算模型一般是由若干不同形狀的三維單元體組成，也即剖分的空間單元網(wǎng)絡(luò)區(qū)，計(jì)算中又將每個(gè)單元體進(jìn)一步劃分成由四個(gè)節(jié)點(diǎn)構(gòu)成的四面體，四面體的應(yīng)力應(yīng)變只通過(guò)四個(gè)節(jié)點(diǎn)向其它四面體傳遞，進(jìn)而傳遞到其它單元體。當(dāng)對(duì)某一節(jié)點(diǎn)施加荷載后，在某一個(gè)微小的時(shí)間段內(nèi)，作用于該點(diǎn)的荷載只對(duì)周?chē)娜舾晒?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)變過(guò)程中，利用高斯積分理論，將三維問(wèn)題轉(zhuǎn)化為二維問(wèn)題而使其簡(jiǎn)

3、單化。在運(yùn)動(dòng)方程中，還充分考慮了巖土體所具有的粘滯性，將其視作阻尼附加于方程中。FLAC3D具有一個(gè)功能強(qiáng)大的網(wǎng)格生成器，有12種基本形狀的單元體可供選擇，利用這12種基本單元體，幾乎可以構(gòu)成任何形狀的空間立體模型。FLAC3D主要是為地質(zhì)工程應(yīng)用而開(kāi)發(fā)的巖土體力學(xué)數(shù)值評(píng)價(jià)計(jì)算程序，自身設(shè)計(jì)有九種材料本構(gòu)模型：（1）空模型（Null Model）（2）彈性各向同性材料模型（Elastic, Isotropic Model）（3）彈性各向異性材料模型（Elastic, anisotropic Model）（4）德拉克-普拉格彈塑性材料模型（Drucker-Prager Model）（5）莫爾-庫(kù)

4、倫彈塑性材料模型（Mohr-Coulomb Model）（6）應(yīng)變硬化、軟化彈塑性材料模型（Strain-Hardening/Softening Mohr-Coulomb Model）（7）多節(jié)理裂隙材料模型（Ubiquitous-Joint Model）（8）雙曲型應(yīng)變硬化、軟化多節(jié)理裂隙材料模型（Bilinear Strain-Hardening/Softening Ubiquitous-Joint Model）（9）修正的Cam粘土材料模型（Modified Cam-clay Model）除上述本構(gòu)模型之外，F(xiàn)LAC3D還可進(jìn)行動(dòng)力學(xué)問(wèn)題、水力學(xué)問(wèn)題、熱力學(xué)問(wèn)題等的數(shù)值模擬。在邊界條件及

5、初始條件的考慮上，F(xiàn)LAC3D軟件十分靈活方便，可在數(shù)值計(jì)算過(guò)程中隨時(shí)調(diào)整邊界條件和初始條件。FLAC3D具有強(qiáng)大的后處理功能，用戶(hù)可以直接在屏幕上繪制或以文件形式創(chuàng)建或輸出打印多種形式的圖形、文字，用戶(hù)還可根據(jù)各自的需要，將若干個(gè)變量合并在同一幅圖形中進(jìn)行研究分析。FLAC3D軟件還可對(duì)各種開(kāi)挖工程或施加支護(hù)工程等進(jìn)行數(shù)值仿真模擬，軟件自身設(shè)計(jì)有錨桿、錨索、襯砌、支架等結(jié)構(gòu)元素，可以直接模擬這些支護(hù)于圍巖（土）體的相互作用。FLAC3D擁有可以自行設(shè)計(jì)的FISH語(yǔ)言，用戶(hù)可根據(jù)自身需求，自己設(shè)計(jì)材料的本構(gòu)模型、屈服準(zhǔn)則、支護(hù)方案、復(fù)雜形狀的開(kāi)挖方式等工作。特別注意的是，巖石是一種脆性材料，當(dāng)

6、外荷載達(dá)到巖石強(qiáng)度后，材料發(fā)生斷裂破壞，產(chǎn)生弱化現(xiàn)象，應(yīng)屬于彈塑性體。在FLAC3D中，一般對(duì)于彈塑性材料，判斷其破壞與否的基本準(zhǔn)則有兩個(gè)，即Drucker-Prager準(zhǔn)則和Mohr-Coulomb準(zhǔn)則。根據(jù)室內(nèi)巖石力學(xué)性質(zhì)試驗(yàn)結(jié)果，其典型應(yīng)力應(yīng)變曲線(xiàn)反映出巖體破壞包絡(luò)線(xiàn)符合莫爾庫(kù)倫屈服準(zhǔn)則，故本次建立的本構(gòu)力學(xué)模型選擇莫爾庫(kù)倫彈塑性材料模型為宜。2.2.2 FLAC3D分析計(jì)算原理計(jì)算所采用的數(shù)學(xué)模型是根據(jù)彈塑性理論的基本原理（應(yīng)變定義、運(yùn)動(dòng)定律、能量守衡定律、平衡方程及理想材料的連續(xù)性方程等）而建立的。2.2.2.1 基本約定在數(shù)學(xué)及數(shù)值模型的表達(dá)式中，符號(hào)有一定的約定含義，一般A表示張

7、量，Aij表示張量A的（i,j）分量，a表示矢量，ai表示矢量a的i分量，a,i表示a對(duì)xi的偏導(dǎo)數(shù)。xi，ui，vi和dvi/dt，（i=1,3）分別表示一點(diǎn)的位置矢量分量、位移矢量分量、速度矢量分量和加速度矢量分量。2.2.2.2 數(shù)學(xué)模型（一）柯西（Cauchy）應(yīng)力張量與柯西公式對(duì)于一個(gè)具有體積V的封閉曲面s的物體，在其上取一表面元素Ds，這個(gè)表面元素的單位外法向矢量為n，在某一時(shí)刻t，在表面元素對(duì)于連續(xù)介質(zhì)中一點(diǎn)，作用著對(duì)稱(chēng)的應(yīng)力張量sij，根據(jù)Ds上作用有力DP，則極限稱(chēng)為表面力。若用ti表示T的分量，則在三維直角坐標(biāo)系中可有關(guān)系式 （1）這個(gè)關(guān)系式稱(chēng)為柯西公式，其中，sij稱(chēng)為柯

8、西應(yīng)力張量。（二）應(yīng)變速率和旋轉(zhuǎn)速率如果介質(zhì)質(zhì)點(diǎn)具有運(yùn)動(dòng)速度矢量v，則在一個(gè)無(wú)限小的時(shí)間dt內(nèi)，介質(zhì)會(huì)產(chǎn)生一個(gè)由vidt決定的無(wú)限小應(yīng)變，對(duì)應(yīng)的應(yīng)變速率分量xij為 （2）而其旋轉(zhuǎn)速率分量wij為 （3）（三）運(yùn)動(dòng)及平衡方程根據(jù)牛頓運(yùn)動(dòng)定律與柯西應(yīng)力原理，如果質(zhì)點(diǎn)作用著應(yīng)力sij與體力bi，且具有速度vi，則在無(wú)限小時(shí)間段dt內(nèi)，它們之間的關(guān)系為 （4）式中，r為質(zhì)點(diǎn)密度。（4）式稱(chēng)為柯西運(yùn)動(dòng)方程。當(dāng)質(zhì)點(diǎn)的加速度為零時(shí)，上式變?yōu)殪o力平衡方程 （5）（四）本構(gòu)方程上述（4）式與（5）式組成的方程組中含有9個(gè)方程，15個(gè)未知量，其中12個(gè)是應(yīng)力與應(yīng)變速率分量，3個(gè)是速度分量。其余6個(gè)關(guān)系式則由本構(gòu)

9、方程提供，本構(gòu)方程一般具有如下形式 （6）式中，為應(yīng)力變化速率，H表示一個(gè)特定的函數(shù)關(guān)系，k為與荷載歷史有關(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)（見(jiàn)圖2-1）。單元面單元節(jié)點(diǎn)圖2-1 四面體單元示意圖（一）空間微分的有限差分逼近對(duì)于一個(gè)計(jì)算單元，若內(nèi)部各質(zhì)點(diǎn)速度為一連續(xù)的矢量場(chǎng)v，則根據(jù)高斯（Gauss）積分原理有 （7）式中，n為外法向單位矢量場(chǎng)。由于單元體的應(yīng)變速率是連續(xù)的，因此可以近似認(rèn)為速度是線(xiàn)性變化的，則（7）式可用下面的求和公式近似逼近 （8）式中，V為單元體體積

10、，為單元體某一面的面積，f為單元體面數(shù)，f=1,4，為面平均速度的i分量。由于速度場(chǎng)是線(xiàn)性的，則有 （9）式中，l為單元體節(jié)點(diǎn)數(shù)，l=1,4；vil為l節(jié)點(diǎn)的i速度分量。將（9）式代入（8）式可得 （10）根據(jù)正交原理有 （11）則將（10）式兩邊除以V，并將（11）式代入可得 （12）因此有 （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í)，其平衡方程為 （14） （15）介質(zhì)可由若干個(gè)作用著體力B，各自產(chǎn)生一定變形的四面體組成，設(shè)節(jié)點(diǎn)力為fn，n=1,4，利用虛功原理，假定單元體節(jié)點(diǎn)具有速度dvn，內(nèi)部具有變形速率dx，則由節(jié)點(diǎn)力fn和體力B所做

11、的外功率與sij所做的內(nèi)功率應(yīng)當(dāng)相等。外功率E可用下式表示 （16）而內(nèi)功率I為 （17）由（13）式，對(duì)于恒定應(yīng)變速率單元體可有 （18）由于應(yīng)力張量是對(duì)稱(chēng)的，定義矢量Tl （19）可得 （20）將（15）式代入（16）式得 （21）式中，Eb和EI分別是體力rbi和慣性力產(chǎn)生的外功率，對(duì)于單元體內(nèi)恒定的體力rbi，Eb可寫(xiě)成 （22）而EI可以寫(xiě)成 （23）如前所述，單元體內(nèi)部速度場(chǎng)以線(xiàn)性變化，為方便描述，選取單元體質(zhì)心為原點(diǎn)，、為坐標(biāo)軸，可有 （24）式中，Nn(n=1,4)是具有如下形式的線(xiàn)性函數(shù) （25）、（n=1,4）則由下面方程決定 （26）dnj為克羅內(nèi)克爾（Kronecker

12、）記號(hào)，根據(jù)質(zhì)心定義可得 （27）用克萊默（Cramers）法則解得，則 （28）而 （29）則 （30）在穩(wěn)定狀態(tài)下，內(nèi)功率與外功率必定相等，則有 （31）上式最后一項(xiàng)，如果單元體內(nèi)r恒定，則根據(jù)質(zhì)心定義有 （32）若將rV/4看作是假定的節(jié)點(diǎn)質(zhì)量mn，則上式變?yōu)?（33）故 （34）根據(jù)牛頓定律有 （35）式中，F(xiàn)i<l>為節(jié)點(diǎn)l所受力的i分量，M<l>為節(jié)點(diǎn)l的質(zhì)量，為節(jié)點(diǎn)l在Fi<l>作用下產(chǎn)生的加速度，nn為包含在全部連續(xù)介質(zhì)中單元節(jié)點(diǎn)總數(shù)。設(shè)質(zhì)點(diǎn)質(zhì)量，不平衡力F<l>為 （36）（三）時(shí)間微分的顯式有限差分對(duì)于上面（35）式，可寫(xiě)成

13、（37）FLAC3D計(jì)算中，對(duì)于上式的逼近是在一個(gè)足夠小Dt時(shí)間段內(nèi)，采用中心差分格式，即上式變?yōu)?（38）同理，節(jié)點(diǎn)的坐標(biāo)差分公式與位移差分公式分別為 （39） （40）（四）本構(gòu)方程的增量形式FLAC3D中，假定在Dt時(shí)間段內(nèi)，速度是不變的，本構(gòu)方程的增量形式可表示為 （41）記 （42）則Deij為Dt時(shí)間段內(nèi)應(yīng)變?cè)隽?。?yīng)力增量由下式確定 （43）式中，由下式定義 （44）而 （45）2.2.2.4 本構(gòu)模型（一）彈塑性理論的增量關(guān)系一般來(lái)講，彈塑性體破壞的判斷準(zhǔn)則具有如下形式 （46）式中，f為一特定的函數(shù)關(guān)系，為n維應(yīng)力矢量i分量。彈塑性體的變形為彈性變形與塑性變形的總和，即 （47

14、）彈性應(yīng)變與應(yīng)力之間的關(guān)系為 （48）其中，Si為線(xiàn)性函數(shù)。塑性變形可用下式表示 （49）式中，l為一常數(shù)，g為的某一函數(shù)關(guān)系。將（47）、（49）式代入（48）式得 （50）對(duì)于新應(yīng)力，仍有 （51）當(dāng)為線(xiàn)性函數(shù)時(shí)，（51）式可寫(xiě)成 （52）其中，。由于，再將（50）式代入（52）式可得 （53）若記則有 （54）因此有 （55）則新應(yīng)力為 （56）（二）莫爾庫(kù)倫（Mohr-Coulomb）模型本構(gòu)關(guān)系本次計(jì)算的本構(gòu)模型選擇的是莫爾庫(kù)倫彈塑性材料模型，彈塑性體可產(chǎn)生彈性及塑性?xún)刹糠肿冃?，根?jù)虎克定律（Hookes law），應(yīng)力應(yīng)變的關(guān)系為 （57）其中，E為剛度矩陣。對(duì)于彈性變形，應(yīng)力增量

15、可由下式確定 （58）式中，a1、a2可由剪切模量與體積模量得出（58）式還可寫(xiě)成 （59）（三）破壞準(zhǔn)則與流動(dòng)法則若，則莫爾庫(kù)倫破壞準(zhǔn)則具有形式f s = 0，其中 （60）對(duì)于拉張破壞，有f t = 0的形式，其中 （61）上兩式中，f為內(nèi)摩擦角，c為內(nèi)聚力，s t為抗拉強(qiáng)度， （62）對(duì)于材料的抗拉強(qiáng)度s t，其值不可能超過(guò)s3（見(jiàn)圖2-2），st的最大值為 （63）對(duì)于塑性流動(dòng)狀態(tài)，設(shè)g s與g t分別為剪切破壞和拉張破壞所對(duì)應(yīng)的函數(shù)關(guān)系，有 （64） （65）式中，y為膨脹角， （66）f t = 0s3s1= 0=As3=s1=Bf s = 0s tc / tan fC圖2-2 莫爾庫(kù)倫破壞準(zhǔn)則示意圖流動(dòng)準(zhǔn)則可定義為函數(shù)關(guān)系，在s 1、s 3平面內(nèi)將和以上部分分為兩個(gè)區(qū)域（見(jiàn)圖2-3），其函數(shù)表達(dá)式如下1區(qū)2區(qū)圖2-3 莫爾庫(kù)倫模型塑性流分區(qū) （67）其中，a p與s p由下式定義 （68）從圖2-3中可見(jiàn)，若某點(diǎn)的應(yīng)力狀態(tài)在1區(qū)時(shí)，則發(fā)生剪切破壞，在2區(qū)時(shí)