胡學(xué)龍，李克慶，璩世杰

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

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

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

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

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

1 巖石彈塑性本構(gòu)模型建立

1.1 應(yīng)力應(yīng)變關(guān)系

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

應(yīng)力（增量）與應(yīng)變（增量）的之間的關(guān)系為：

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

1.2 統(tǒng)一屈服準(zhǔn)則

當(dāng)巖石的應(yīng)力狀態(tài)達到其初始屈服極限時，巖石內(nèi)部便開始產(chǎn)生塑性。本文中以統(tǒng)一強度理論作為屈服準(zhǔn)則，用三個主應(yīng)力σ1、σ2和σ3（σ1≥σ2≥σ3）表示為（以拉應(yīng)力為正）：

式中：a 為巖石的拉壓強度比；b 為反映中間主應(yīng)力對巖石破壞影響的參數(shù)，其取值范圍為0≤b≤1（外凸理論）；σt為巖石單軸抗拉強度；σc為巖石單軸抗壓強度。a 和b 取不同參數(shù)時，統(tǒng)一屈服準(zhǔn)則可以演變?yōu)橐幌盗衅渌?zhǔn)則，例如，當(dāng)a=1 且b=0 時該準(zhǔn)則變?yōu)門resca 屈服準(zhǔn)則；當(dāng)0＜a＜1 且b=0 時該準(zhǔn)則變?yōu)镸ohr-Coulomb 準(zhǔn)則；當(dāng)a=1 且b=1 時該準(zhǔn)則變?yōu)殡p剪屈服準(zhǔn)則。統(tǒng)一屈服準(zhǔn)則在偏平面上的軌跡如圖1 所示。

圖 1 統(tǒng)一屈服準(zhǔn)則函數(shù)在偏平面上的軌跡Fig. 1 The locus of unified strength theory on deviatoric plane

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

式中：c 為巖石的準(zhǔn)靜態(tài)內(nèi)聚力，φ 為巖石的內(nèi)摩擦角。

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

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

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

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

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

式中：cd為巖石動態(tài)內(nèi)聚力；fDIF為動態(tài)增長因子，動態(tài)增長因子定義為巖石的動態(tài)強度與巖石的準(zhǔn)靜態(tài)強度之比。本文中fDIF表示為：

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

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

1.3 塑性勢函數(shù)

由塑性力學(xué)可知，塑性應(yīng)變可表示為：

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

由文獻[17]可知，統(tǒng)一強度理論的塑性勢函數(shù)可表示為：

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

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

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

聯(lián)立式(15)～(16)可得：

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

聯(lián)立式(9)～(10)和(19)可得廣義塑性剪切應(yīng)變增量：

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

式中：

聯(lián)立式(4)、(15)和(21)可得巖石本構(gòu)關(guān)系可表示為：

2 模型數(shù)值實現(xiàn)

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

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

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

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

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

式中：Δσc,n+1為修正應(yīng)力，其值為：

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

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

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

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

3 模型驗證與分析

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

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

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

圖 3 巖石彈塑性本構(gòu)模型數(shù)值實現(xiàn)流程Fig. 3 Flow chart of numerical implementation of material constitutive model

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

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

圖 5 算例1 石灰?guī)r單軸壓縮應(yīng)力應(yīng)變曲線Fig. 5 Stress-strain curves of limestone uniaxial compression in example 1

3.2 算例2 巖石動態(tài)壓縮實驗

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

圖 6 算例2 內(nèi)聚力c 與廣義剪切塑性應(yīng)變γp 之間的關(guān)系Fig. 6 Relation between cohesion c and generalized shear plastic strain γp in example 2

圖 7 fDIF 與加載應(yīng)變率之間的關(guān)系Fig. 7 Relation between fDIF and Loading rate

圖10 為石灰?guī)r準(zhǔn)靜態(tài)下的單軸應(yīng)力應(yīng)變曲線，從圖中不難看出模擬結(jié)果與試驗結(jié)果高度一致，應(yīng)力應(yīng)變曲線峰前有一定的偏差是因為數(shù)值模擬中所選取的彈性模量為平均彈性模量。

圖11 為利用石灰?guī)r試樣進行SHPB 實驗得到的應(yīng)變時間信號，從圖11 中可以看到，數(shù)值模擬結(jié)果與實驗結(jié)果比較吻合，這說明本文本構(gòu)模型的正確性，能夠反映巖石在動載作用下的力學(xué)行為。

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

圖 9 巖石SHPB 數(shù)值實驗?zāi)Ｐ虵ig. 9 Numerical model of rock SHPB test

圖 10 算例2 石灰?guī)r準(zhǔn)靜態(tài)下單軸壓縮應(yīng)力應(yīng)變曲線Fig. 10 Stress-strain curves of limestone quasi-static uniaxial compression in example 2

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

4 結(jié) 論

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

（2）采用巖石單軸壓縮試驗（準(zhǔn)靜態(tài)）和巖石SHPB 試驗（動態(tài)）兩個算例對彈塑性本構(gòu)模型進行驗證，結(jié)果表明，該本構(gòu)模型能夠刻畫巖石在準(zhǔn)靜態(tài)和動態(tài)下的力學(xué)行為。