考慮損傷的軟弱夾層剪切流變模型及程序實現(xiàn)

祝 鑫 ，胡 斌 ，李 京 ，崔 凱 ，魏二劍 ，劉 楊

(1. 武漢科技大學 資源與環(huán)境工程學院，湖北 武漢 430081； 2. 冶金礦產(chǎn)資源高效利用與造塊湖北省重點實驗室，湖北 武漢 430081)

巖土體的流變性質復雜，其本構模型研究是流變力學理論研究的一大難點、熱點問題，同時也是數(shù)值試驗研究的模型基礎。國內外許多學者經(jīng)過對巖體微細觀破裂機制的不斷研究，發(fā)現(xiàn)巖體在荷載作用下的應變與破壞主要表現(xiàn)為斷裂形式與損傷形式，且越來越多的學者更傾向于應用損傷的觀點來構建相應的流變模型。趙延林等[1]將伯格斯模型和摩爾庫侖塑性損傷元件串聯(lián)建立了新的蠕變損傷模型；楊圣奇等[2]結合損傷力學與有效應力觀點，推導了巖石在兩階段中的損傷演化方程，建立了能較好描述巖石流變三階段的損傷流變模型；Ma等[3]應用Lemaitre應變等效理論引入損傷變量，提出了基于廣義Hoek-Brown準則的非線性損傷模型；胡波等[4]考慮蠕變起始的應力閾值，將一種階躍函數(shù)表示的開關與廣義K體并聯(lián)，并引入損傷變量對CVISC模型進行了改進；徐衛(wèi)亞等[5]在流變不同階段引入不同函數(shù)和損傷變量，建立了綠片巖的改進的廣義Bingham模型；伍國軍等[6]建立了適用于工程巖體的黏彈塑性損傷本構模型；陳鐵林等[7]將結構性黏土看作由結構體和結構面雙重介質組成，在堆砌體模型的基礎上，采用過應力理論，建立了結構性黏土的流變模型；張強勇等[8]將巖體流變力學參數(shù)看作是非定常的，構建了變參數(shù)的蠕變損傷本構模型；Yang等[9]根據(jù)不同圍壓下煤的短期和蠕變試驗結果，提出了新的損傷演化方程，建立了一個新的非線性蠕變損傷模型；陳衛(wèi)忠等[10]建立了鹽巖蠕變損傷的本構方程和損傷演化方程；張亮亮等[11]推導出損傷變量的演化方程，用于黏塑性模型中的黏性參數(shù)，得到了一種新的巖石損傷蠕變模型。

隨著計算機軟硬件的不斷發(fā)展，各種數(shù)值分析軟件應運而生，適用于巖土體領域的軟件也與日俱增，使得巖土體的數(shù)值研究不斷涌現(xiàn)出新成果。但是，實際工程巖體具有非均質性，其各向異性明顯，力學性質復雜，與數(shù)值模擬軟件自帶的本構模型差別較大。因此，進行自定義本構模型的開發(fā)成為了一種研究趨勢，目前已取得一定的研究成果。胡浩[12]提出一種基于分數(shù)階微積分的本構模型，基于FLAC3D軟件進行開發(fā)，并用隧道注漿加固的算例驗證了其有效性；褚衛(wèi)江等[13]給出了西原流變損傷模型在FLAC3D中實施的程序流程；徐平等[14]研制了廣義Kelvin模型的FLAC3D接口程序，并對基坑開挖進行了實例分析；楊文東等[15]基于FLAC3D開發(fā)了改進Burgers蠕變損傷模型，并驗證了黏彈性、塑性和損傷的力學性質；謝秀棟等[16]開發(fā)了軟土彈黏塑性流變模型，通過模擬側向卸荷流變試驗，并與修正劍橋模型作比較，驗證了二次開發(fā)的可行性；何利軍等[17]使用C++語言實現(xiàn)了含SMP強度準則的黏彈塑性模型的開發(fā)，為其他強度準則的開發(fā)提供了參考；蔣邦友等[18]同時考慮巖石材料的塑性及損傷共同作用，完成了基于Mogi-Coulomb準則的巖石彈塑性損傷本構模型的程序實現(xiàn)。

在自然界中，巖體多數(shù)處于擠壓、剪切的狀態(tài)，其變形和破壞的主要形式是沿著滑面的剪切破壞，且主要表現(xiàn)為剪切流變損傷特征，對巖體破壞主要考慮的因素是剪切應力[19]。本文通過對軟弱夾層不同剪切應力水平下的試驗數(shù)據(jù)進行分析，建立一個考慮流變損傷的軟弱夾層剪切流變模型，并基于FLAC3D軟件，使用C++語言進行自定義模型的代碼編寫，最后驗證開發(fā)模型的正確性。

1 剪切流變損傷本構模型建立

軟弱夾層的剪切流變一般經(jīng)歷減速、穩(wěn)定及加速流變階段。其中，加速流變階段的變形速率會越來越快，變形迅速增大，且無法收斂，最終發(fā)生破壞?，F(xiàn)有的軟弱夾層流變模型的建立主要有兩種方法，其一是在經(jīng)典流變模型的基礎上增加新的非線性元件進行組合；其二是將損傷力學等新理論、觀點引入本構模型中，以此建立非線性流變模型。

本文基于Kachanov流變損傷理論[20]，引入流變損傷變量D來描述流變過程中軟弱夾層參數(shù)的變化。損傷變量D是流變過程中巖體材料不斷劣化，強度不斷減小的量化指標，其具體表現(xiàn)形式為：

式中：t為時間。當t=tD（臨界破壞時間）時，D=1，則可得：

由式（2）、（3）可得損傷演化方程為：

圖1 為隨時間增加，不同a取值情況下流變損傷變量D的變化，可以看到，a取值不同，流變損傷變量D的變化快慢也不相同，但是在時間t趨近破壞時間tD時，其變化率趨近于無窮大，此時D趨近于1，表明已完全損傷，處于破壞狀態(tài)。

圖1 a不同取值下?lián)p傷變量D的趨勢曲線Fig. 1 Trend curve of damage variable D under different values of a

基于流變損傷變量D，提出一個能夠描述軟弱夾層加速流變階段的黏彈塑性流變模型，見圖2所示。

圖2 黏彈塑性流變模型Fig. 2 Viscoelastic plastic rheological model

軟弱夾層的剪切流變是黏、彈、塑性等多種變形共同作用的過程，要使用基本線性元件及其他非線性元件的串并聯(lián)組合形式進行描述，對于減速、穩(wěn)定兩個階段可以采用常規(guī)的伯格斯模型進行描述，對于加速流變階段需要采用能夠描述黏彈塑性特征的模型進行模擬。因此，本文將黏彈塑性流變模型與伯格斯模型串聯(lián)起來，可以同時描述軟弱夾層的減速、穩(wěn)定及加速流變3個階段的流變過程（圖3）。

圖3 非線性流變損傷模型Fig. 3 Nonlinear rheological damage model

由模型的串聯(lián)關系可知，建立的軟弱夾層非線性流變損傷組合模型滿足以下條件：

式中：τ1、τ2、τ3分 別為模型A、B及C組件剪應力（MPa）；ξ 為總應變量，ξ1、ξ2、ξ3分別為模型A、B及C組件產(chǎn)生的應變量。

根據(jù)剪切屈服應力τs與 τ 之間的關系，所建立的軟弱夾層非線性流變損傷模型的流變方程需分為以下兩種情況：

（1）τ ＜τs時 ，模型中C組件不工作，流變方程為：

式中：k1、η1、k2、η2分別為A和B組件中的彈性系數(shù)（GPa）及黏性系數(shù)（GPa·h）。

（2）τ ≥τs時，C組件開始工作，流變方程為：

式中：k3、η3分別為C組件的彈性系數(shù)（GPa）和黏性系數(shù)（GPa·h）。

式（9）可以同時描述軟弱夾層的減速流變、穩(wěn)定流變和加速流變3個階段的流變變形過程。

2 流變損傷模型的三維差分形式

為了將本文所提出的損傷本構模型導入FLAC3D，需要將模型的流變方程改寫為FLAC3D所能識別差分格式，主要包括偏應力增量和流變時間的關系及應變增量與流變時間的關系。

由模型串聯(lián)規(guī)則可以得到，總偏應變增量為A組件馬克斯韋爾體、B組件開爾文體、C組件黏彈塑性體3個偏應變增量的和。即：

式中：ij表示各個分量方向分別為馬克斯韋爾體、開爾文體、黏彈塑性體的偏應變增量。

對于A組件的馬克斯韋爾體，由彈性元件與黏性元件串聯(lián)而成，其偏應變增量為：

對于B組件開爾文體，由彈性元件與黏性元件并聯(lián)而成，其偏應力與模型總偏應力Sij相等，且與其偏應變平均值偏應變增量的關系為：

式中：k2、 η2分別為開爾文體的彈性系數(shù)（GPa）和黏性系數(shù)（GPa·h）。

又因為式中：

求解可得：

將式（11）、（13）、（16）、（18）代入式（10）并化簡可得：

其中：

模型中的球應力張量不發(fā)生塑性應變，滿足下列關系式：

最終的總應力張量可以分解為球應力張量σm和 偏應力張量Sij，即：

3 剪切流變損傷模型二次開發(fā)

流變損傷本構模型的開發(fā)在Microsoft Visual Studio 2015軟件環(huán)境下完成。由于在FLAC3D中，損傷變量D形式較為復雜，無法寫成流變時間的增量差分格式，因此，損傷變量的定義采用FLAC3D軟件自帶的FISH語言進行編寫，具體為使用FISH語言將本文推導的損傷演化方程定義為一個函數(shù)，隨著流變時間的變化，損傷變量也發(fā)生變化，進而相應的模型參數(shù)也發(fā)生變化，以此實現(xiàn)參數(shù)的損傷演化過程。

本文開發(fā)的軟弱夾層剪切流變損傷模型主要包括：（1）模型開發(fā)插件、開發(fā)平臺軟件及解決方案的生成；（2）模型的名稱、注冊ID、版本號、自定義變量k1、n1、k2、n2、k3、n3及部分中間變量等頭文件的定義與修改；（3）核心Initialize（）函數(shù)、Run（）函數(shù)流變計算指令的輸入、偏應力核心公式（20）的編寫及其公式中必要參數(shù)P、Q、m、n的準備、計算過程中的單元剪切屈服準則的代碼編入、模型剪切破壞的識別及應力修正；（4）FLAC3D軟件可以識別、調用的動態(tài)鏈接庫DLL文件。在完成主要的頭文件和源文件的編譯和修改后，將其重新生成解決方案，并儲存在Release文件夾中，命名為modelmknp.dll之后，將其拷貝到FLAC3D安裝地址中的exe64文件夾中以供使用。

4 剪切流變損傷本構模型驗證

4.1 計算模型

為驗證本文提出的非線性軟弱夾層剪切流變損傷本構模型數(shù)值模擬的正確性，選擇文獻[21]中室內剪切流變試驗不同剪切荷載的應變時間曲線來進行對比驗證，模型尺寸大小、加載條件及方式、加載時間等條件相同。數(shù)值計算的模型形狀為正方體，邊長為150 mm，共劃分為1 000個單元、1 331個節(jié)點。具體的數(shù)值計算模型及加載條件見圖4。

由圖4可得，該數(shù)值模型的下半部分剪切盒全固定，不能產(chǎn)生位移，而對上半部施加剪切荷載，剪切荷載分5級從0.10 MPa增加到0.59 MPa，且每級剪切荷載保持的時長為34 h，上表面再施加軸向荷載，軸向荷載為0.5 MPa，由文獻[21]的試驗數(shù)據(jù)可知，該條件下對應的屈服強度為0.423 MPa，臨界破壞時間為34 h。流變本構模型采用自定義的modelmknp模型，共進行5級不同剪切應力水平下的剪切流變數(shù)值模擬計算，對試件的頂部左側頂點M（0，0，0.15）進行水平切向變形監(jiān)測，數(shù)值試驗計算采用的參數(shù)見表1。

圖4 數(shù)值試驗模型Fig. 4 Numerical test model

表1 軟弱夾層力學參數(shù)Tab. 1 Mechanical parameters of weak interlayer

4.2 數(shù)值計算結果分析

根據(jù)不同剪切應力水平下的數(shù)值模擬試驗計算，獲得了該軟弱夾層在5級剪切荷載下的非線性流變損傷模型參數(shù)的最優(yōu)組合，最優(yōu)組合的確定具體為：在每一級剪切應力水平下，先通過不斷模擬調整參數(shù)取值，得到一組擬合效果較好的參數(shù)組合，在此基礎上，再進行參數(shù)的細微調整，得到10組模型參數(shù)組合，進行每組參數(shù)的模擬，將得到的位移時間曲線與室內試驗位移時間曲線進行對比，選擇曲線最接近、相關度最高的參數(shù)組合即為最優(yōu)組合（表2）。

表2 不同剪切應力下流變模型參數(shù)Tab. 2 Parameters of creep model under different shear stress levels

圖5 為數(shù)值模擬試驗結果的塑性區(qū)分布，圖中shear-p表示單元在試驗過程中受剪，tension-n表示單元直到試驗加載結束時受拉，tension-p表示單元在試驗過程中受拉的塑性區(qū)分布情況，即上半部模型右下角單元在試驗過程中受剪，而其他單元在試驗過程中既有受拉也有受剪，且直到試驗加載結束時受拉。

圖5 數(shù)值試驗結果塑性區(qū)分布Fig. 5 Plastic zone distribution of numerical test result

圖6 為5級剪切應力下的試驗曲線和非線性流變損傷本構模型中M點的位移曲線的對比。由圖6可見，在相同的加載時長下，模型每級剪應力的應變值變化都較為接近，且每級剪應力下流變變形的變化趨勢也大致相同，總體上吻合度較高，但也存在部分差異性。

圖6 數(shù)值曲線與試驗數(shù)據(jù)曲線Fig. 6 Numerical and experimental data curves

當施加的剪切荷載較小而未到達屈服應力，即在前4級剪切應力水平下時，其應變變形都只有減速流變階段和穩(wěn)定流變階段，應變的變化速率先不斷減小然后保持不變，應變變形量逐漸穩(wěn)定，趨向于一個不變的值。而當施加的水平方向的剪切應力較大，超過其屈服應力閾值即0.59 MPa時，其應變變形階段會歷經(jīng)減速過渡到穩(wěn)定最后再加速共3個階段；應變速率先減小后不變再迅速增大；應變值不斷增大，不會趨于一個定值，直至破壞。對于每級剪應力下應變的增量值，其瞬時應變都較為接近，且經(jīng)過相同時間的應力加載后，無論是穩(wěn)定后的應變值還是加速階段的應變值都較為接近。差異性主要表現(xiàn)為本文數(shù)值模擬試驗在施加的剪切荷載較小而未到達屈服應力，即在前4級剪切應力水平下時得到的位移時間曲線，其減速階段歷時較室內試驗減速階段歷時短，即更快地達到穩(wěn)定流變階段；而在施加的水平方向的剪切應力較大，超過其屈服應力閾值時得到的位移時間曲線，其加速流變階段位移變化速率較室內試驗位移變化速率大，且最終的位移比室內試驗位移稍大。因此，可以說明本文開發(fā)的非線性剪切流變損傷本構模型是可靠的。

5 結 語

（1）本文通過軟弱夾層不同剪切應力水平下的試驗數(shù)據(jù)分析，引入可以表征其流變損傷的變量D，提出了一個基于D的可以反映軟弱夾層加速流變特性的黏彈塑性非線性流變模型，將該模型與伯格斯模型串聯(lián)，構成了能全面反映3個流變階段的新的軟弱夾層剪切流變損傷模型。

（2）推導出了模型在三維情況下的差分格式，并基于FLAC3D使用C++語言對該剪切流變損傷模型進行了二次開發(fā)與試驗驗證，證明了模型的可靠性。