層狀巖體各向異性彈塑性本構(gòu)模型的數(shù)值實(shí)現(xiàn)

高敏，賈善坡，劉曉東，張力偉，舒婧曦 （長江大學(xué)城市建設(shè)學(xué)院，湖北荊州434023）

層狀巖體在物理、力學(xué)特性上較均質(zhì)巖體具有更顯著的各向異性、非線性及非連續(xù)性，其變形與破壞性質(zhì)十分復(fù)雜。國內(nèi)外學(xué)者對層狀巖體的各向異性特征及本構(gòu)關(guān)系進(jìn)行了大量研究。Jeager［1］和Tien等［2，3］開展了大量的理論與試驗(yàn)研究，這些研究成果為層狀巖體本構(gòu)模型的進(jìn)一步研究奠定了基礎(chǔ)。文獻(xiàn) ［4～6］在橫觀各向同性理論與巖石材料各向異性理論的基礎(chǔ)上，通過試驗(yàn)測定了橫觀各向同性體的5個(gè)彈性參數(shù)，并建立了巖體的彈性模量與泊松比隨節(jié)理面傾角變化的函數(shù)關(guān)系，實(shí)現(xiàn)了層狀巖體等效彈性參數(shù)的各向異性特征。賈善坡等［7］在考慮結(jié)構(gòu)面的幾何特性和力學(xué)性質(zhì)的基礎(chǔ)上建立了巖體各向異性本構(gòu)模型。黃書嶺等［8］在考慮一組節(jié)理面的層狀巖體復(fù)合材料模型的基礎(chǔ)上，建立了考慮多組結(jié)構(gòu)特性的層狀巖體多節(jié)理本構(gòu)模型。雖然眾多學(xué)者做了大量研究工作，但是目前將層狀巖體的彈性與塑性的各向異性特征在有限元軟件中實(shí)現(xiàn)的研究較少。下面，筆者在ABAQUS自帶節(jié)理材料模型的基礎(chǔ)上，結(jié)合巖體各向異性彈性參數(shù)的研究成果，利用ABAQUS中USDFLD子程序接口，采用FORTRAN語言進(jìn)行開發(fā)并嵌入到ABAQUS軟件中實(shí)現(xiàn)非線性計(jì)算功能，通過實(shí)際算例分析驗(yàn)證了該模型的有效性和可靠性。

1 層狀巖體各向異性彈性參數(shù)研究

相關(guān)學(xué)者［4～6］在巖石橫觀各向同性理論和巖石材料的各向異性理論的基礎(chǔ)上，通過試驗(yàn)測定了橫觀各向同性體的5個(gè)彈性參數(shù)，并建立了巖體的彈性模量和泊松比隨節(jié)理面傾角變化的函數(shù)關(guān)系。

Amadei［9］認(rèn)為，含層理面的巖石應(yīng)力（σx，σy，σz，τyz，τzx，τxy）與應(yīng)變（εx，εy，εz，γyz，γyx，γxy）的關(guān)系如下：

由各向異性理論可知，軸向應(yīng)力σy與應(yīng)變（εx，εy，εz）存在如下關(guān)系：

式中，a12、a22、a23是由層理面傾角θ和5個(gè)彈性參數(shù)E1、μ1、E2、μ2、G2決定的函數(shù)關(guān)系式，其表達(dá)式分別為：

式中，E1、μ1分別為平行于橫觀各向同性面的彈性模量與泊松比；E2、μ2、G2分別為垂直于橫觀各向同性面的彈性模量、泊松比和剪切模量。上述5個(gè)彈性參數(shù)可以通過巖體的單軸及三軸壓縮試驗(yàn)測得。

只要確定了5個(gè)彈性參數(shù)E1、μ1、E2、μ2、G2，即可確定任一節(jié)理面傾角的彈性參數(shù)E和μ，即：

這樣就可以得到節(jié)理巖體彈性模量與泊松比隨層理面傾角θ變化的函數(shù)關(guān)系式：

2 子程序開發(fā)步驟

在ABAQUS自帶節(jié)理材料模型的基礎(chǔ)上，結(jié)合層狀巖體彈性參數(shù)的各向異性，以FORTRAN語言為平臺，開發(fā)了本構(gòu)模型子程序，以ABAQUS軟件為求解器使該模型可以進(jìn)行數(shù)值計(jì)算，進(jìn)而用于工程實(shí)踐中。

子程序的開發(fā)主要包括以下4個(gè)部分：確定層狀巖體5個(gè)等效彈性參數(shù)；調(diào)用每一步計(jì)算結(jié)果的應(yīng)力矩陣；調(diào)用每一步計(jì)算結(jié)果最大主應(yīng)力方向的余弦值；將式（7）的函數(shù)關(guān)系嵌入子程序，其中式（7）中的θ是表示在試驗(yàn)過程中改變的節(jié)理面的傾角，其余弦值與實(shí)際工程中最大主應(yīng)力的方向余弦值是對應(yīng)的，最大主應(yīng)力方向余弦值在ABAQUS軟件中通過調(diào)用SPIND函數(shù)得到。

數(shù)值計(jì)算步驟如下：利用ABAQUS自帶節(jié)理材料模型，在初始條件和初始彈性參數(shù)下進(jìn)行初始增量步計(jì)算，并利用計(jì)算結(jié)果以及調(diào)用開發(fā)的子程序計(jì)算得到新的彈性參數(shù)，再利用新的彈性參數(shù)進(jìn)行下一增量步計(jì)算，以此循環(huán)直至所有計(jì)算步驟完成（見圖1）。

圖1 程序計(jì)算步驟

3 數(shù)值算例

3.1 算例模型

選取某巖質(zhì)邊坡進(jìn)行數(shù)值試驗(yàn)（具體尺寸見圖2）。邊坡右側(cè)和坡腳下左側(cè)均為豎向滑動(dòng)水平約束邊界條件，底面為全部固定邊界。土體容重γ＝25kN／m3，粘聚力c＝0.05MPa，內(nèi)摩擦角φ＝45°。彈性模量E1＝37.79GPa，E2＝24.39GPa，泊松比μ1＝0.254，μ2＝0.180，剪切模量10.85GPa［6］。

3.2 計(jì)算工況

1）工況1 對比分析各向同性彈性模型與各向異性彈性模型計(jì)算結(jié)果的差異性，各向同性彈性模型彈性模量與泊松比分別取為E1＝37.79GPa、μ1＝0.254，各向異性彈性參數(shù)按式（7）選取。

2）工況2 對比分析ABAQUS自帶節(jié)理材料模型與各向異性彈塑性模型的計(jì)算結(jié)果，彈性參數(shù)的選取同工況1，其中節(jié)理面傾角θ從0～90°變化。

3）工況3 選取文獻(xiàn) ［6］中層狀巖體彈性參數(shù)E和μ隨節(jié)理面傾角θ變化的試驗(yàn)結(jié)果，與各向異性彈塑性本構(gòu)模型計(jì)算結(jié)果得到的彈性參數(shù)E和μ進(jìn)行對比分析。

3.3 計(jì)算結(jié)果分析

1）工況1計(jì)算結(jié)果分析 工況1計(jì)算結(jié)果如圖3所示。對比各向同性彈性計(jì)算結(jié)果與各向異性彈性結(jié)算結(jié)果可知，圖3（a）與圖3（c）的Mises應(yīng)力圖分布大致相同，但是圖3（c）中在邊坡角處存在著變化的不均勻性；圖3（b）與圖3（d）中的位移云圖分布差異較大，且圖3（d）圖的變化更能體現(xiàn)出位移變化的不均勻性，這說明各向異性彈性數(shù)值模型能夠在一定程度上反映出各向異性特征。

圖2 計(jì)算模型尺寸圖

圖3 工況1計(jì)算結(jié)果圖

2）工況2計(jì)算結(jié)果分析 圖4所示為ABAQUS自帶節(jié)理材料模型計(jì)算結(jié)果（節(jié)理面傾角θ從0～90°變化，其中選取10、30、60、80°的計(jì)算結(jié)果）。從圖4可以看出，隨著節(jié)理面傾角的變化其Mises應(yīng)力與位移變化不大，且云圖分布比較均勻；塑性破壞區(qū)域呈現(xiàn)出先增大再減小再增大的趨勢，且塑性破壞區(qū)域變化不均勻，呈現(xiàn)出一定的各向異性特征。

圖5所示為選取層狀巖體各向異性彈塑性本構(gòu)模型計(jì)算結(jié)果（節(jié)理面傾θ從0～90°變化，其中選取10、30、60、80°計(jì)算結(jié)果）。對比各向異性彈塑性模型與ABAQUS自帶節(jié)理材料模型計(jì)算結(jié)果可知，塑性破壞區(qū)與Mises應(yīng)力云圖差異不大，但各向異性彈塑性數(shù)值模型計(jì)算的結(jié)果偏大，特別是位移云圖的分布具有明顯的不均勻性。上述分析表明，各向異性彈塑性數(shù)值模型能夠在一定程度上準(zhǔn)確表現(xiàn)出巖石材料的各向異性特征。

圖4 ABAQUS自帶節(jié)理材料計(jì)算結(jié)果圖

3）工況3計(jì)算結(jié)果分析 通過數(shù)值計(jì)算得到的E和μ隨節(jié)理面傾角θ變化的結(jié)果與文獻(xiàn) ［6］中試驗(yàn)得到的結(jié)果對比如圖6所示。從圖6可以看出，通過數(shù)值計(jì)算得到的彈性模量E隨節(jié)理面傾角的變化與試驗(yàn)值比較接近且變化規(guī)律具有一致性，隨著節(jié)理面傾角θ的增加而增加。另外，雖然通過數(shù)值計(jì)算得到的泊松比μ與試驗(yàn)值有一定的誤差，但是也隨著節(jié)理面傾角θ的增加而增加。

4 結(jié)語

在ABAQUS自帶節(jié)理材料模型基礎(chǔ)上，結(jié)合層狀巖體彈性參數(shù)的各向異性研究，采用FORTRAN語言編制了本構(gòu)子程序，使層狀巖體各向異性彈塑性模型得以數(shù)值實(shí)現(xiàn)。以某巖質(zhì)邊坡為數(shù)值算例，采用多種工況對比分析，驗(yàn)證了層狀巖體各向異性彈塑性數(shù)值模型的可靠性，從而為實(shí)際工程中研究層狀巖體的各向異性特征提供數(shù)值計(jì)算基礎(chǔ)。層狀巖體各向異性彈塑性數(shù)值模型雖然能夠反映巖體的各向異性特征，但是如何考慮多組節(jié)理面存在并將其運(yùn)用到復(fù)雜模型及實(shí)際工程中有待今后進(jìn)一步的研究。

圖5 各向異性彈塑性本構(gòu)模型計(jì)算結(jié)果

圖6 工況3計(jì)算結(jié)果對比

［1］Jaeger J C.Shear failure of anisotropic rocks ［J］.Geology Magazine，1960，97（1）：65～72.

［2］Tien Y M，Kuo M C.A failure criterion for transversely isotropic rocks ［J］.International Journal of Rock Mechanics and Mining Sciences，2001，38（5）：399～412.

［3］Tien Y M，Kuo M C，Juang C.An experimental investigation of the failure mechanism of simulated transversely isotropic rocks ［J］.International Journal of Rock Mechanics and Mining Sciences，2006，43（8）：1163～1181.

［4］席道英，陳林 .巖石各向異性參數(shù)研究 ［J］.物探化探計(jì)算技術(shù)，1994，16（1）：16～21.

［5］張學(xué)民 .巖石材料各向異性特征及其對隧道圍巖穩(wěn)定性影響研究 ［D］.長沙：中南大學(xué)，2007.

［6］劉運(yùn)思，傅鶴林，伍毅敏，等 .橫觀各向同性巖石彈性參數(shù)及抗壓強(qiáng)度的試驗(yàn)研究 ［J］.中南大學(xué)學(xué)報(bào)（自然科學(xué)版），2013，44（8）：3398～3404.

［7］賈善坡，鄒臣頌，吳渤，等 .考慮巖體各向異性強(qiáng)度的隧洞圍巖穩(wěn)定性分析 ［J］.防災(zāi)減災(zāi)工程學(xué)報(bào)，2013，33（2）：179～184.

［8］黃書嶺，丁秀麗，鄔愛清，等 .層狀巖體多節(jié)理本構(gòu)模型與試驗(yàn)驗(yàn)證 ［J］.巖石力學(xué)與工程學(xué)報(bào)，2010，29（4）：743～756.

［9］Amadei B.Importance of anisotropy when estimating and measuring in situ stresses in rock ［J］.International Journal of Rock Mechanics and Mining Sciences，1996，33（3）：293～325.