基于統一強度理論的巖石彈塑性本構模型及其數值實現*

胡學龍，李克慶，璩世杰

（1. 北京科技大學土木與資源工程學院，北京 100083；2. 莫納什大學土木工程學院，澳大利亞 墨爾本 3800）

統一強度理論是以一個統一物理模型為基礎，囊括了所有應力分量以及它們對材料破壞的不同影響，能夠適用于各種巖石類材料，Mohr-Coulomb 強度理論和雙剪強度理論均為其特例，并且還包含了可以比DP 準則更合理的新的計算準則以及可以描述非凸極限面試驗結果的新的非凸強度理論[1]，而且能與巖石材料的的真三軸試驗結果相吻合[2]，從而在巖土工程領域得到廣泛應用[3]。例如，廖紅建等[4]基于統一強度理論研究了巖土材料在復雜應力狀態下的強度理論及確定其動力強度參數的新方法。李杭州等[5]通過引入洛德參數得到了統一強度參數，從而建立了基于統一強度理論的軟巖損傷統計本構模型。張強等[6]基于深部巖體良好的塑性變形能力和高地應力下瞬時破壞特性，建立了適用深部巖體力學行為的彈塑脆性模型。曹雪葉等[7]以統一強度理論為屈服準則，經過推導得到了凍結壁的彈塑性應力場、彈性極限荷載及塑性極限荷載的解析解。統一強度理論由于其自身的優越性，在巖土界甚至其他領域得到了越來越多的應用，受到了越來越多學者的歡迎與重視。

以統一強度理論作為屈服準則的巖石彈塑性本構模型已經被諸多學者研究，比如，Yu 等[8]以統一強度理論為屈服準則，采用了相關聯和非相關聯流動法則建立了巖土材料的雙剪統一彈塑性模型，從而使彈塑性模型可以使用不同的屈服準則來模擬各種不同巖土材料的應力應變關系，但該模型并沒有考慮巖石的硬化/軟化規律，只屬于理想彈塑性模型。張傳慶等[9]將基于統一強度理論的彈塑性模型與FLAC3D數值分析軟件結合起來，推導出了統一彈塑性本構模型在FLAC3D中的計算格式，但該模型也只是理想彈塑性本構模型，并沒有考慮巖土的硬化/軟化規律。潘曉明等[3]把基于統一強度理論的彈塑性本構模型引入到通用有限元軟件ABAQUS 中，使得基于統一強度理論的彈塑性本構模型更能在巖土領域中得到應用，但其彈塑性本構模型只是把巖石的硬化函數用一個常數來代替，這與實際巖土的硬化特性不符。李杭州等[10]根據統一強度理論，以駝峰型曲線作為硬化函數，建立了可以考慮應變硬化和應變軟化的統一彈塑性模型，由于其沒有與通用有限元軟件結合起來，不便于其推廣和應用。

LS-DYNA 主要用于求解三維非彈性結構在高速碰撞、爆炸沖擊下的大變形動力響應問題[11]。雖然LS-DYNA 自身提供了很多巖土材料模型（例如FWHA 模型），但目前發現很少有文獻把基于統一強度理論的考慮應變率效應的彈塑性本構模型導入到LS-DYNA 中，這就很大程度上阻礙了該模型在碰撞、沖擊和爆破領域中的應用。若能把二者結合起來，一方面使得基于統一強度理論的彈塑性本構模型能運用到實際中去，另一方面也能豐富LS-DYNA 的材料庫，提高LS-DYNA 的計算能力。

本文中以統一強度理論為屈服準則，視巖石強度由摩擦強度和內聚強度兩部分組成，在摩擦強度不變的基礎上認為內聚強度是廣義剪切塑性應變的函數，并引入應變率函數，首先建立考慮應變軟化和應變率效應的巖石彈塑性本構模型，然后利用LS-DYNA 的用戶自定義材料本構程序接口(UMAT)，把該本構模型嵌入到LS-DYNA 中去，最后通過巖石單軸壓縮試驗和巖石SHPB 試驗兩個算例驗證該模型的正確性。

1 巖石彈塑性本構模型建立

1.1 應力應變關系

由彈塑性力學可知，應變（增量）由彈性應變（增量）和塑性應變（增量）組成，即：式中：ε 為應變列向量，εe為彈性應變列向量，εp為塑性應變列向量，dε 為應變增量列向量，dεe為彈性應變增量列向量，dεp為塑性應變增量列向量。

應力（增量）與應變（增量）的之間的關系為：

式中：σ 為應力列向量，dσ 為應力增量列向量，D 為彈性剛度矩陣。

1.2 統一屈服準則

當巖石的應力狀態達到其初始屈服極限時，巖石內部便開始產生塑性。本文中以統一強度理論作為屈服準則，用三個主應力σ1、σ2和σ3（σ1≥σ2≥σ3）表示為（以拉應力為正）：

式中：a 為巖石的拉壓強度比；b 為反映中間主應力對巖石破壞影響的參數，其取值范圍為0≤b≤1（外凸理論）；σt為巖石單軸抗拉強度；σc為巖石單軸抗壓強度。a 和b 取不同參數時，統一屈服準則可以演變為一系列其他屈服準則，例如，當a=1 且b=0 時該準則變為Tresca 屈服準則；當0＜a＜1 且b=0 時該準則變為Mohr-Coulomb 準則；當a=1 且b=1 時該準則變為雙剪屈服準則。統一屈服準則在偏平面上的軌跡如圖1 所示。

圖 1 統一屈服準則函數在偏平面上的軌跡Fig. 1 The locus of unified strength theory on deviatoric plane

巖石的單軸抗拉強度可以表示為：

式中：c 為巖石的準靜態內聚力，φ 為巖石的內摩擦角。

巖石強度通常被認為由兩部分組成，即內聚強度和摩擦強度[12]。巖石在載荷作用下變形過程中可以認為巖石的摩擦強度為定值，即內摩擦角保持不變；內聚強度是廣義塑性剪切應變的函數，即內聚力可以表示成廣義塑性剪切應變的函數[13-14]。根據文獻中的試驗數據進行擬合，內聚力可以表示為：

式中：A、B、C 和D 為擬合參數，可以通過繪制內聚力與廣義塑性剪切應變之間的關系圖從而進行擬合確定，具體確定過程可以參考文獻[14]；γp為廣義塑性剪切應變；廣義塑性剪切應變可以通過下式計算：

式中：dγp為廣義塑性剪切應變增量，dep為偏塑性應變增量。有：

式中：dεv為體積塑性應變增量，I 為單位矩陣。

在不同應變率荷載下巖石表現出不同的力學行為，也就是說巖石是一種應變率依賴型的地質材料。應變率對巖石力學行為最直接的影響即是使巖石的強度增加。根據文獻[15-16]可知應變率對巖石的摩擦強度影響較小，可以忽略不計；應變率對巖石的內聚強度影響較大，巖石動態內聚強度是應變率和準靜態內聚強度的函數。該函數可以表示為

式中：cd為巖石動態內聚力；fDIF為動態增長因子，動態增長因子定義為巖石的動態強度與巖石的準靜態強度之比。本文中fDIF表示為：

聯立式(7)～(8)、(11)～(12)可得巖石的動態單軸抗拉強度為：

用式(13)中的σtd替換式(5)中的σt即可得到考慮巖石應變硬化/軟化行為和應變率效應的統一屈服準則，即：

1.3 塑性勢函數

由塑性力學可知，塑性應變可表示為：

式中：dλ 為塑性乘子，dλ≥0；G 為塑性勢函數。

由文獻[17]可知，統一強度理論的塑性勢函數可表示為：

式中：a*=(1-sinψ)/(1+sinψ)，ψ 為膨脹角，若ψ=φ，則為關聯塑性流動法則，否者則為非關聯塑性流動法則。巖土材料一般采用非關聯塑性流動法則。

由彈塑性力學可知，在屈服面上必有：

對式（17）兩邊同時微分可得：

聯立式(15)～(16)可得：

式中：dεp,1、dεp,2和dεp,3分別為第一主塑性應變、第二主塑性應變和第三主塑性應變。

聯立式(9)～(10)和(19)可得廣義塑性剪切應變增量：

聯立式(4)、(15)和(20)可得塑性乘子：

式中：

聯立式(4)、(15)和(21)可得巖石本構關系可表示為：

2 模型數值實現

模型數值實現是一個在已知tn時刻的應力σn、應變增量Δεn+1、廣義塑性剪切應變γp,n求出tn+1時刻應力σn+1的過程。這里采用應力返回算法來達到求解的目的。應力返回算法主要分為以下幾步：

（1）彈性預測：σtr,n+1=σn+D?εn+1式中：σtr,n+1為試探應力。

（2）屈服判斷：將tn+1時刻的試探應力σtr, n+1代入到式(14)中，如果F≤0 說明巖石處于彈性狀態，此時有σn+1=σtr, n+1；如果F＞0 說明此時巖石已經進入到塑性屈服狀態，應該使應力狀態返回到屈服面。

（3）應力返回：為了使試探應力返回到屈服面，本文中采用割平面法(CPA)，它的幾何原理如圖2 所示。由式(17)可知在屈服面上有：

對式(24)在試探應力σtr, n+1處進行一階泰勒展開可得：

式中：Δσc,n+1為修正應力，其值為：

由式(25)～(26)可得：

如圖2 中紅線所示，當加載步較大時，應力從F(σtr, n+1)返回到屈服面（F(σn+1)=0）的過程中不可能一次完成，而是需要分k 次迭代，直至|F(σn+1)|≤Eallow，Eallow為允許誤差，本文中取Eallow=1.0×10?4。此時有：

式(29)中的(dλ)k可由下式求得：

在LS-DYNA 中，用戶能夠根據自己的需要采用Fortran 語言來編寫材料本構模型的子程序，通過其提供的用戶自定義材料本構程序接口(UMAT)來生成求解器，從而對材料的本構模型進行求解。本文中彈塑性本構模型的計算流程如圖3 所示。

3 模型驗證與分析

為了驗證巖石彈塑性模型的正確性，本文中主要通過巖石單軸壓縮試驗（準靜態）和巖石SHPB 試驗（動態）兩個算例來進行說明。

3.1 算例1：巖石單軸壓縮數值模擬試驗

圖 2 CPA 應力返回映射算法幾何示意圖Fig. 2 Geometric illustration of CPA stress return mapping algorithms

圖 3 巖石彈塑性本構模型數值實現流程Fig. 3 Flow chart of numerical implementation of material constitutive model

本算例的試驗數據來自zhang 等[18]的研究工作，巖石式樣為圓柱形石灰巖，其直徑為50 mm、高100 mm。巖石密度ρ=2 720 kg/m3、彈性模量E=44.76 GPa、泊松比ν=0.33、內摩擦角φ=50°，擬合參數A=132.7 MPa、B=?63.49、C=?117.8 MPa、D=?83.17。內聚力c 與廣義剪切塑性應變γp的關系如圖4 所示。

圖5 為通過數值模擬得到的單軸壓縮應力應變曲線與試驗結果的對比，由圖5 可以看到應力應變曲線的數值模擬結果與試驗結果比較一致。應力峰值后的軟化規律與試驗結果有一定的出入，這是由于擬合得到的內聚力c 曲線的峰值相較于試驗數據峰值向右有一定的偏移（見圖4），因此數值模擬得到的峰值后的軟化曲線整體均向右偏移一定距離，但是應力應變的峰值以及峰值后的變化趨勢與試驗結果均有很好的吻合度。

圖 5 算例1 石灰巖單軸壓縮應力應變曲線Fig. 5 Stress-strain curves of limestone uniaxial compression in example 1

3.2 算例2 巖石動態壓縮實驗

本算例的實驗數據來自Frew 等[19-20]和Liao 等[21]對石灰巖的研究工作，實驗中的巖石試樣為圓柱形石灰巖，其直徑與高均為12.7 mm。巖石彈性模量E=24 GPa、密度ρ=2 300 kg/m3、泊松比ν=0.23、內摩擦角φ=25°, 擬合參數A=22.11 MPa、B=?25.64、C=?5.095 MPa、D=?3 594，應變率參數l=0.352 7、m=0.165 6。圖6 為內聚力c 與廣義剪切塑性應變γp關系擬合圖，圖7 為動態增長因子fDIF與應變率的關系擬合圖。按照文獻[19]中所述來建立SHPB 試驗模型，SHPB 實驗裝置由三部分組成，即撞擊桿、入射桿和透射桿，SHPB 試驗裝置如圖8 所示，它們的直徑與巖石試樣直徑相同，長度分別為152、2 130、915 mm，入射桿和透射桿上各安裝一個應變計用來測量桿中的應變時間信號，它們的位置如圖8 所示。所建SHPB 試驗有限元模型如圖9 所示。撞擊桿、入射桿和透射桿均為VM350 鋼制成，其彈性模量Eb=200 GPa、泊松比νb=0.23、密度8 100 kg/m3、屈服強度為2 500 MPa。撞擊桿的初始速度為8.05 m/s。

圖 6 算例2 內聚力c 與廣義剪切塑性應變γp 之間的關系Fig. 6 Relation between cohesion c and generalized shear plastic strain γp in example 2

圖 7 fDIF 與加載應變率之間的關系Fig. 7 Relation between fDIF and Loading rate

圖10 為石灰巖準靜態下的單軸應力應變曲線，從圖中不難看出模擬結果與試驗結果高度一致，應力應變曲線峰前有一定的偏差是因為數值模擬中所選取的彈性模量為平均彈性模量。

圖11 為利用石灰巖試樣進行SHPB 實驗得到的應變時間信號，從圖11 中可以看到，數值模擬結果與實驗結果比較吻合，這說明本文本構模型的正確性，能夠反映巖石在動載作用下的力學行為。

圖 8 SHPB 實驗裝置示意圖Fig. 8 Illustration of SHPB test device

圖 9 巖石SHPB 數值實驗模型Fig. 9 Numerical model of rock SHPB test

圖 10 算例2 石灰巖準靜態下單軸壓縮應力應變曲線Fig. 10 Stress-strain curves of limestone quasi-static uniaxial compression in example 2

圖 11 利用石灰巖試樣進行SHPB 實驗的應變時程曲線Fig. 11 Strain time history curve for split Hopkinson pressure bar experiment with a limestone sample

4 結 論

（1）基于彈塑性力學理論，建立了巖石的彈塑性本構模型，該彈塑性本構模型一方面描述了巖石的硬化/軟化行為，另一方面反映了巖石的應變率效應。

（2）采用巖石單軸壓縮試驗（準靜態）和巖石SHPB 試驗（動態）兩個算例對彈塑性本構模型進行驗證，結果表明，該本構模型能夠刻畫巖石在準靜態和動態下的力學行為。