基于FLAC3D黃土的濕載結(jié)構(gòu)性本構(gòu)模型二次開發(fā)

羅愛忠， 方 娟

(貴州工程應(yīng)用技術(shù)學(xué)院 土木建筑工程學(xué)院， 貴州 畢節(jié) 551700)

21世紀(jì)已經(jīng)成為土體本構(gòu)模型特別是結(jié)構(gòu)性土體本構(gòu)模型研究的世紀(jì)[1-5]，針對(duì)各種土特性的土體本構(gòu)模型已經(jīng)被提出，但是由于沒有和計(jì)算程序的有機(jī)結(jié)合，限制了這些模型在實(shí)際工程實(shí)踐中的應(yīng)用。近年來，為了使現(xiàn)有已開發(fā)的本構(gòu)模型得到更好的應(yīng)用和被工程界所接受，眾多的學(xué)者開展了現(xiàn)有本構(gòu)模型與現(xiàn)有通用軟件的有機(jī)結(jié)合，這其中，又以有限差分法類軟件為代表。王春波等[1]推導(dǎo)了硬化土本構(gòu)模型的有限差分格式，與現(xiàn)有有限差分通用商業(yè)軟件FLAC相結(jié)合進(jìn)行了二次開發(fā)，并將其應(yīng)用于深基坑穩(wěn)定性分析。李英杰等[6]考慮了變形模量的劣化，建立了相應(yīng)的應(yīng)變軟化模型，結(jié)合FLAC進(jìn)行隧道開挖穩(wěn)定性分析。楊文東等[7]改進(jìn)了Burgers蠕變損傷模型，并在FLAC軟件中實(shí)現(xiàn)了二次開發(fā)。左雙英等[8]在研究現(xiàn)有彈塑性本構(gòu)模型的基礎(chǔ)上，考慮了Zienkiewicz-Pan強(qiáng)度準(zhǔn)則，最后基于有限差分法進(jìn)行了大型地下洞室的開挖穩(wěn)定性分析。蘭航等[9]為了分析節(jié)理巖體的采動(dòng)損傷問題，建立了節(jié)理巖體采動(dòng)損傷本構(gòu)模型，在FLAC3D中分析了井工開采對(duì)上覆巖體的損傷問題。陶慧等[10]將堆石料三維邊界面模型嵌入FLAC3D中，進(jìn)行了堆石料穩(wěn)定性分析。何利君等[11]將黏彈塑性模型與SMP準(zhǔn)則結(jié)合，并在FLAC3D中進(jìn)行了二次開發(fā)。陳育民等[12]將鄧肯-張模型引入FLAC3D中，并進(jìn)行了一系列三軸試驗(yàn)驗(yàn)證。綜合已有研究，有限差分法由于其概念相對(duì)直觀，再加上其數(shù)學(xué)表達(dá)式描述得相對(duì)簡(jiǎn)單，在巖土計(jì)算領(lǐng)域得到了很大的運(yùn)用[13-17]。有限差分法直接將所要求解的計(jì)算域劃分為有限差分網(wǎng)格代替復(fù)雜的連續(xù)求解域，用網(wǎng)格節(jié)點(diǎn)上的函數(shù)的差商代替了偏微分方程中復(fù)雜的偏導(dǎo)數(shù)，從而建立以節(jié)點(diǎn)函數(shù)值為未知量的代數(shù)方程組。濕陷性黃土是一種典型的非飽和結(jié)構(gòu)性黃土，結(jié)構(gòu)性的存在使得其在工程實(shí)際中表現(xiàn)出了不同于正常固結(jié)土的特殊性質(zhì)，為了描述黃土結(jié)構(gòu)性的存在對(duì)應(yīng)力應(yīng)變關(guān)系及工程特性的影響，驗(yàn)證證明了該模型在描述黃土應(yīng)力應(yīng)變關(guān)系方面的準(zhǔn)確性和優(yōu)越性及應(yīng)用于工程實(shí)際分析，有必要將該本構(gòu)模型與通用計(jì)算程序相結(jié)合。本文以黃土的濕載結(jié)構(gòu)性本構(gòu)模型的數(shù)值化為切入點(diǎn)，以黃土壓剪濕力學(xué)作用變化過程分析為目的，驗(yàn)證壓黃土的濕載結(jié)構(gòu)性本構(gòu)模型的合理性和準(zhǔn)確性及本構(gòu)模型與大型數(shù)值計(jì)算軟件FLAC3D結(jié)合的可行性與有效性，進(jìn)而為模型的實(shí)際應(yīng)用奠定基礎(chǔ)。

1 FLAC3D分析軟件的數(shù)值格式及二次開發(fā)

1.1 FLAC3D拉格朗日顯示差分法

差分法最早由著名學(xué)者Wilkins在1963年提出[18]，這一差分格式主要是依據(jù)偏導(dǎo)數(shù)積分的定義而得到。

在FLAC3D中，單元格式為空間單元，因此其差分格式遵從空間的差分格式，其應(yīng)變率分量的差分格式可以表述為式(1)形式，F(xiàn)LAC3D中的顯式計(jì)算中，時(shí)間增量步的確定可以通過周期及假設(shè)的節(jié)點(diǎn)質(zhì)量確定。

(1)

1.2 FLAC3D分析軟件的數(shù)值格式

FLAC3D是美國(guó)Itasca Consulting Group公司開發(fā)的一款用于巖土工程數(shù)值模擬的軟件。該軟件建立在拉格朗日算法的基礎(chǔ)上，采用有限差分顯式算法來獲得模型全部運(yùn)動(dòng)方程的時(shí)間步長(zhǎng)解，該技術(shù)和混合離散劃分技術(shù)的結(jié)合，保證了非常精確地模擬塑性破壞和流變，從而追蹤材料的漸進(jìn)破壞，它適合于模擬計(jì)算巖土材料力學(xué)行為，特別適合模擬材料的高度非線性、不可逆剪切破壞和壓密等問題。

FLAC3D中，所有的本構(gòu)模型都遵從同樣的算法。這種算法是在給定t時(shí)刻的應(yīng)力狀態(tài)和Δt時(shí)間步的應(yīng)變?cè)隽浚瑥亩蟮胻+Δt時(shí)刻的應(yīng)力增量和應(yīng)力狀態(tài)。假設(shè)發(fā)生的變形全部是彈性變形，根據(jù)胡克定律直接求出應(yīng)力增量。如果塑性變形發(fā)生，根據(jù)彈塑性一般原理，應(yīng)力增量根據(jù)總應(yīng)變?cè)隽康膹椥圆糠执_定。即，在計(jì)算中在每一時(shí)步都必須對(duì)假設(shè)得出的應(yīng)力增量進(jìn)行塑性修正，從而得出新的應(yīng)力增量和應(yīng)力狀態(tài)。

(1) 破壞準(zhǔn)則。破壞準(zhǔn)則滿足下式：

f(σn)=0

(2)

式中：f為屈服函數(shù)，即當(dāng)塑性流動(dòng)發(fā)生時(shí)的應(yīng)力組合關(guān)系，該函數(shù)在應(yīng)力空間中是一個(gè)曲面，所有在曲面之內(nèi)的應(yīng)力點(diǎn)保持彈性行為；σn表示σ1，σ2，…，σn。

(2) 應(yīng)變一般分解為彈性和塑性兩部分

(3)

式中：i分別取1、2、3；上標(biāo)e和p分別表示彈性和塑性。

(3) 應(yīng)變?cè)隽颗c應(yīng)力增量之間的彈性關(guān)系

(4)

(4) 流動(dòng)法則。如果采用關(guān)聯(lián)流動(dòng)法則，及塑性勢(shì)函數(shù)和屈服函數(shù)采用相同的形式，則可以得到流動(dòng)法則滿足如下的關(guān)系式：

(5)

式中：λ為常數(shù)；g為塑性勢(shì)函數(shù)；σ為應(yīng)力。

(5) 屈服函數(shù)的塑性修正。經(jīng)過塑性修正后，新的屈服函數(shù)滿足如下的關(guān)系：

f(σn+Δσn)=0

(6)

新的塑性應(yīng)變?cè)隽糠匠瘫挥脕碓u(píng)估塑性應(yīng)變?cè)隽康拇笮　?/p>

同時(shí)考慮Si的線性關(guān)系，得：

(7)

通過流動(dòng)法則，更進(jìn)一步應(yīng)變?cè)隽靠梢员硎緸椋?/p>

(8)

式中：Si與應(yīng)變?cè)隽恐g的線性關(guān)系同樣被采用。

硬化常數(shù)λ表示為：

(9)

1.3 FLAC3D自定義本構(gòu)模型的二次開發(fā)環(huán)境及模型文件編寫

在FLAC3D中，可以通過兩種方法自定義本構(gòu)模型，一種方法是采用FISH語言，另外一種方法是采用面向?qū)ο蟮恼Z言C++編寫。FISH語言繁瑣而且不通用，因而C++就成了FLAC3D中自定義本構(gòu)模型的主流方法。模型嵌入的方法是將C++編譯的模型文件在現(xiàn)有的編譯平臺(tái)上編譯成動(dòng)態(tài)鏈接庫文件(DLL文件)，在使用時(shí)通過FLAC3D的可執(zhí)行文件載入。自定義的本構(gòu)模型的主要功能是給定應(yīng)變?cè)隽?，獲得新的應(yīng)力，輔助功能還包括提供模型的名稱、版本等基本信息及完成讀寫的基本操作。

目前FLAC3D中自定義本構(gòu)模型需要至少在Visual Studio 2005以上的版本進(jìn)行創(chuàng)建，在FLAC3D中所有的本構(gòu)模型的加載幾乎都是通過動(dòng)態(tài)鏈接庫的形式提供的，采用這種鏈接方式主要具有以下兩個(gè)優(yōu)點(diǎn)：使得用戶自定義的本構(gòu)模型與軟件自帶的本構(gòu)模型兩者在執(zhí)行時(shí)效率處于同一水平；自定義的本構(gòu)模型在更高版本的FLAC3D系列軟件中得到通用。

2 黃土濕載結(jié)構(gòu)性本構(gòu)模型在FLAC3D中的計(jì)算格式推導(dǎo)

文獻(xiàn)[19]提出的考慮壓剪濕綜合影響的結(jié)構(gòu)性本構(gòu)模型是基于一定的應(yīng)力增量得到塑性應(yīng)變?cè)隽浚窃贔LAC3D中，其計(jì)算格式是基于一定的應(yīng)變?cè)隽康玫较鄳?yīng)的應(yīng)力增量。根據(jù)彈塑性理論的基本知識(shí)，彈性應(yīng)變?cè)隽亢退苄詰?yīng)變?cè)隽恐g的關(guān)系可以表示為：

(10)

(11)

依據(jù)相關(guān)聯(lián)流動(dòng)法則，可以得到塑性應(yīng)變?cè)隽勘硎緸椋?/p>

(12)

由文獻(xiàn)[19]的屈服函數(shù)，可以得到屈服函數(shù)中球應(yīng)力和偏應(yīng)力的偏導(dǎo)如下：

(13)

(14)

根據(jù)Hooke定律，應(yīng)力增量與應(yīng)變?cè)隽康年P(guān)系可以表示為：

(15)

(16)

根據(jù)相關(guān)聯(lián)流動(dòng)法則，新的應(yīng)力狀態(tài)可以表示為：

pn=po+Δp

(17)

qn=qo+Δq

(18)

式(18)中，上標(biāo)n和o分別代表新的應(yīng)力和塑性修正前的應(yīng)力。將式(16)和式(17)分別代入式(17)和式(18)得：

pn=pI-Λ*KCa

(19)

qn=qI-Λ*3GCb

(20)

其中，pI和qI分別為：

pI=po+Kdεv

(21)

qI=qo+3Gdεs

(22)

把經(jīng)過塑性修正后的應(yīng)力代入屈服面方程，并進(jìn)一步整理則有：

aΛ2+bΛ+c=0

(23)

其中：

進(jìn)一步對(duì)式(23)求解，可得方程的兩個(gè)根，若兩根都大于零，去最小的一個(gè)根，若兩根一正一負(fù)，取正根。從而，新的應(yīng)力可以表示為：

(24)

對(duì)于濕載條件下的結(jié)構(gòu)性邊界面本構(gòu)模型，由于引入下加載屈服面，不能按常規(guī)方程退化求出式(24)，結(jié)合相關(guān)聯(lián)流動(dòng)法則，硬化參量改寫為：

(25)

進(jìn)一步，在主應(yīng)力空間表示為：

(26)

依據(jù)文獻(xiàn)[19]濕載條件下的結(jié)構(gòu)性邊界面本構(gòu)模型的屈服函數(shù)，有如下的關(guān)系式成立：

(27)

式(27)中，各分項(xiàng)滿足：

其中，δij為克羅內(nèi)爾系數(shù)。結(jié)合式(27)，則有：

式(10)～式(27)中，p，q，σ，ε為土的球應(yīng)力、偏應(yīng)力、一般應(yīng)力和一般應(yīng)變，下標(biāo)i和j表示變量順序；δij為克羅內(nèi)爾系數(shù)；psx為p-q平面上結(jié)構(gòu)屈服面與p軸正半軸交點(diǎn)坐標(biāo)；pt為p-q平面上結(jié)構(gòu)屈服面與p軸負(fù)半軸交點(diǎn)坐標(biāo)；M為p-q平面上剪切屈服面斜率。

3 黃土的濕載結(jié)構(gòu)性本構(gòu)模型的實(shí)現(xiàn)流程及加載運(yùn)行

FLAC3D中自定義本構(gòu)模型的工作空間格式為UDM.dsw，在該工作平臺(tái)下，以C++編程語言為手段，通過編輯頭文件，張量文件及運(yùn)行文件的塑性指示及流動(dòng)法則，實(shí)現(xiàn)應(yīng)力應(yīng)變的計(jì)算和應(yīng)力狀態(tài)變量的更新。基于上一節(jié)的公式推導(dǎo)，結(jié)合FLAC3D特點(diǎn)，得到黃土的濕載結(jié)構(gòu)性本構(gòu)模型的數(shù)值實(shí)現(xiàn)流程如圖1所示。

圖1程序?qū)崿F(xiàn)流程圖

模型加載運(yùn)行過程中，參數(shù)的輸入必須與程序中一致，見圖2。將編好的頭文件loess.h和源文件loess.cpp導(dǎo)入到工程文件中，經(jīng)過編譯鏈接之后，形成了動(dòng)態(tài)鏈接庫文件loesscam.dll，將loesscam.dll拷貝到安裝文件目錄下，如圖3所示。在程序調(diào)用時(shí)，通過config cppudm進(jìn)行模塊聲明，通過model load loesscam.dll語句進(jìn)行模型加載。判斷模型是否被加載成功，可以通過print model 語句判斷是否加載成功，如果加載成功將會(huì)命令窗口中找到該本構(gòu)模型及其編號(hào)，如圖4所示。這樣，程序在計(jì)算時(shí)就可以識(shí)別自定義的模型及屬性。需要注意的是，在使用RESTORE命令恢復(fù)一個(gè)含有自定義模型的文件時(shí)，也必須使用config cppudm命令和model load loesscam.dll語句。

圖2 C++程序中的模型參數(shù)

圖3 FLAC3D生成動(dòng)態(tài)鏈接庫文件

圖4自定義本構(gòu)模型加載

4 自編本構(gòu)模型合理性驗(yàn)證

模型的調(diào)試與驗(yàn)證實(shí)際上包含了兩個(gè)最主要的方面：一方面是程序編制過程中出現(xiàn)的各種語法及算法的合理性，對(duì)于這一方面的問題，參照程序語言的編程規(guī)定即可實(shí)現(xiàn)；另一方面，需要對(duì)自定義本構(gòu)模型的計(jì)算結(jié)果進(jìn)行驗(yàn)證和完善，這個(gè)問題的解決方法一般是將模型的計(jì)算結(jié)果與理論值進(jìn)行比較，以證明自定義本構(gòu)模型的合理性和正確性。為了驗(yàn)證自定義本構(gòu)模型內(nèi)嵌FLAC3D平臺(tái)的可行性及合理性，本文采應(yīng)變式加載三軸試驗(yàn)的方法驗(yàn)證模型開發(fā)的正確性。通過FLAC3D外部程序語言Fish語言編寫了軸向應(yīng)力，軸向應(yīng)變及體應(yīng)變、剪應(yīng)變求解的外部函數(shù)。模擬的關(guān)鍵在于與黃土常規(guī)三軸試驗(yàn)排水條件的相似性。采用的模擬參數(shù)見表1，該模型參數(shù)通過常規(guī)三軸試驗(yàn)求得，鑒于篇幅關(guān)系，可參看文獻(xiàn)[19]。圖5～圖8分別給出了含水率10%、固結(jié)壓力50 kPa和含水率18%、固結(jié)壓力200 kPa的試驗(yàn)結(jié)果及自定義模型模擬結(jié)果的對(duì)比關(guān)系。從圖5～圖8可以看出，自定義本構(gòu)模型與試驗(yàn)得到的試驗(yàn)值比較吻合，變化趨勢(shì)一致，大小相差不大，誤差在5%以內(nèi)，這進(jìn)一步證明了筆者所建立的黃土的濕載結(jié)構(gòu)性本構(gòu)模型的合理性及其在FLAC3D平臺(tái)上實(shí)現(xiàn)的可行性和合理性。

表1 計(jì)算參數(shù)

圖5 含水率10%，固結(jié)壓力50 kPa試驗(yàn)結(jié)果

圖6 含水率10%，固結(jié)壓力50 kPa自定義本構(gòu)模型模擬結(jié)果

圖7含水率18%，固結(jié)壓力200 kPa試驗(yàn)結(jié)果

參考文獻(xiàn)：

[1] 謝定義,姚仰平,黨發(fā)寧.高等土力學(xué)[M].北京：高等教育出版社,2008.

[2] Terzaghi K. Theoretical Soil Mechanics[M]. New York: J. Wiley and Sons, 1943．

圖8含水率18%，固結(jié)壓力200 kPa自定義本構(gòu)模型模擬結(jié)果

[3] Soga K. Mechanical behavior and constitutive modeling of natural structured soils[D]. Berkeley: Univ. of California, Calif. 1994.

[4] Martin D. Liu, Buddhima N. Indraratna, M., General strength criterion for geomaterials including anisotropic effect[J]. International Journal of Geomechanics, 2011,11(3):251-262.

[5] 金福喜，錢 毅，袁權(quán)威.某城市排土場(chǎng)邊坡三維穩(wěn)定性分析[J].水資源與水工程學(xué)報(bào)，2017,28(2):222-226.

[6] 王春波,丁文其,喬亞飛.硬化土本構(gòu)模型在FLAC～(3D)中的開發(fā)及應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào),2014,33(1):199-208.

[7] 李英杰,張頂立,劉保國(guó).考慮變形模量劣化的應(yīng)變軟化模型在FLAC～(3D)中的開發(fā)與驗(yàn)證[J].巖土力學(xué),2011,32(S2):647-652，659．

[8] 楊文東,張強(qiáng)勇,張建國(guó),等.基于FLAC～(3D)的改進(jìn)Burgers蠕變損傷模型的二次開發(fā)研究[J].巖土力學(xué),2010,31(6):1956-1964.

[9] 左雙英,肖 明,陳俊濤.基于Zienkiewicz-Pande屈服準(zhǔn)則的彈塑性本構(gòu)模型在FLAC～(3D)中的二次開發(fā)及應(yīng)用[J].巖土力學(xué),2011,32(11):3515-3520.

[10] 藍(lán) 航,姚建國(guó),張華興,等.基于FLAC～(3D)的節(jié)理巖體采動(dòng)損傷本構(gòu)模型的開發(fā)及應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào),2008,27(3):572-579.

[11] 何利軍,吳文軍,孔令偉.基于FLAC～(3D)含SMP強(qiáng)度準(zhǔn)則黏彈塑性模型的二次開發(fā)[J].巖土力學(xué),2012,33(5):1549-1556.

[12] 陶 惠,陳育民,肖 楊,等.堆石料三維邊界面模型在FLAC～(3D)中的開發(fā)與驗(yàn)證[J].巖土力學(xué),2014,35(6):1801-1808.

[13] 陳育民,劉漢龍.鄧肯-張本構(gòu)模型在FLAC～(3D)中的開發(fā)與實(shí)現(xiàn)[J].巖土力學(xué),2007,28(10):2123-2126.

[14] 鄒 博,姚仰平,路德春.變換應(yīng)力三維化方法在清華模型中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào),2005(23):4303-4307.

[15] 楊慶光，李毛毛，王國(guó)鋒，等．考慮土體結(jié)構(gòu)性損傷的靜壓樁承載力及時(shí)效性研究[J]．湖南工業(yè)大學(xué)學(xué)報(bào)，2017,31(6):13-19.

[16] 李忠友，劉元雪，李永毅，等．結(jié)構(gòu)性土在能量轉(zhuǎn)化過程中的損傷演化規(guī)律[J]．土木建筑與環(huán)境工程，2017,49(5):40-48.

[17] 胡小榮，樊曉梅，童肖龍，等．飽和粘性土的結(jié)構(gòu)性本構(gòu)模型[J]．南昌大學(xué)學(xué)報(bào)(工學(xué)版)，2017,39(3):246-253.

[18] Alder B J, Femback S, Rotenberg M. Fundamental methods in dydrodynamics[J]. Methods in Computational Physics Advances in Research & Applications, 1964:211-263.

[19] 羅愛忠,邵生俊,陳昌祿,等.黃土的濕載結(jié)構(gòu)性本構(gòu)模型研究[J].巖土力學(xué),2015,36(8):2209-2215.