禹海濤，胡曉錕，李天斌

（1. 成都理工大學地質災害防治與地質環(huán)境保護國家重點實驗室，四川 成都 610059；2. 同濟大學土木工程防災國家重點實驗室，上海 200092；3. 同濟大學土木工程學院，上海 200092）

巖石作為一種天然介質體，具有顯著的彈塑性變形特征和斷裂力學行為，如何構建一個能夠合理描述巖石彈塑性變形特征和連續(xù)-非連續(xù)力學行為的數(shù)值模型，對于深入研究巖石在復雜荷載條件下的物理力學行為至關重要。

近年來學者們提出了不同類型的巖石本構模型，其中Hoek-Brown（H-B）準則是應用最為廣泛、影響最大的巖石強度準則［1］，主要通過對大量的巖石試驗數(shù)據進行統(tǒng)計分析提出的。相比于Mohr-Coulomb（M-C）準則和Drucker-Prager（D-P）準則，H-B準則能夠綜合考慮巖塊的結構面強度以及巖體結構的影響，Hoek 等［2］建立了與地質強度指數(shù)（GSI）相關的參數(shù)確定方法，因此可以合理地反映巖石強度的非線性破壞特征。Cundall等［3］建立了基于H-B 準則的巖石彈塑性本構模型，根據應力狀態(tài)和巖石變形的破壞特點建立了不同的流動法則，并將其應用于FLAC軟件中進行實現(xiàn)。由于H-B準則的屈服面存在棱線和尖點，在數(shù)值求解中需要進行適當?shù)奶幚恚热鏗-B 準則屈服函數(shù)尖點處的平滑處理［4］，Clausen 等［5］采用主應力空間的返回映射方法并結合有限元進行彈塑性模擬分析，解決了棱線和尖點難以收斂的問題。

然而，目前H-B 準則的研究主要基于傳統(tǒng)連續(xù)介質力學框架，比如有限元法、有限差分法，對于巖石斷裂破壞等不連續(xù)問題的研究較為欠缺。近場動力學（PD）［6］作為一種新的非局部的理論，通過空間積分方程代替?zhèn)鹘y(tǒng)連續(xù)介質力學的偏微分方程，避免了模擬不連續(xù)問題時的奇異性，因此非常適合描述巖石等準脆性材料的斷裂損傷力學行為。Yu等［7］建立了廣義微極近場動力學模型來模擬黏彈性問題，并提出了相應的漸近斷裂準則，能夠合理表征準脆性材料在不同加載速率下的非線性黏彈性行為與斷裂模式。Zhou 等［8］建立了基于Drucker-Prager 準則的常規(guī)態(tài)型近場動力學模型來研究巖土材料的塑性變形行為，并模擬了巖石介質的裂紋萌生和擴展問題［9］。雖然近場動力學已經應用于彈塑性斷裂問題的模擬，但是缺少對巖石彈塑性力學行為的準確描述。此外，由于鍵型［6］和常規(guī)態(tài)型近場動力學［10］均缺少與材料對應關系的描述，故難以實現(xiàn)材料彈塑性本構模型的通用表述，而非常規(guī)態(tài)型近場動力學模型［11］中包含了與連續(xù)介質力學相對應的變形梯度張量，便于應用巖石相關的彈塑性本構關系，因此，有必要結合廣泛應用的H-B 準則和非常規(guī)態(tài)型近場動力學模型建立巖石的非局部彈塑性斷裂力學模型，以期合理表征巖石彈塑性連續(xù)變形特征以及斷裂非連續(xù)力學行為。

本文提出一種基于H-B強度準則的非常規(guī)態(tài)型近場動力學彈塑性模型，旨在描述巖石的彈塑性變形連續(xù)場特征以及斷裂力學不連續(xù)場力學行為。通過非局部理論框架下主應力空間的返回映射方法，模擬巖石的彈塑性變形特征，并基于建立的等效塑性應變斷裂準則實現(xiàn)復雜荷載條件下巖石彈塑性斷裂問題的數(shù)值模擬。基于數(shù)值模擬結果與有限元和試驗數(shù)據的對比分析，驗證該模型的正確性。

1 非常規(guī)態(tài)型近場動力學理論框架

近場動力學基于非局部理論的思想，將模型離散為物質點的形式，并認為一個物質點不僅與它鄰近點發(fā)生相互作用，還會受到整個近場非局部作用區(qū)域Hx內其他物質點的影響，且物質點之間相互作用以長程力的形式表征，通過對鄰域內積分可得近場動力學的運動方程為［6］

式中：D為彈性剛度矩陣。進而可得第一類Piola-Kirchhoff（P-K）應力張量和柯西應力張量

式中：J為雅可比矩陣行列式的值。根據第一類P-K應力與應變能和變形梯度的關系，得到近場動力學的力態(tài)表達式

式中：E為彈性模量；υ為泊松比；h為二維問題的厚度；h1為一維問題的截面面積。

2 基于Hoek-Brown 準則的近場動力學彈塑性模型

2.1 屈服方程

Hoek 和Brown［15］在1980 年通過對大量試驗進行分析得到了H-B準則，1992年Hoek等［16］對H-B準則進行改進，提出了適用更廣泛的廣義H-B準則，為

式中：σ1、σ3分別為最大和最小主應力（以拉應力為正）；σc為巖石的單軸抗壓強度；a為經驗參數(shù)，反映材料的非線性程度；mb為經驗參數(shù)，是巖石的材料常數(shù)，表示巖石軟硬程度；s為經驗參數(shù)，與巖體的完整性有關，當巖體完整性較好時，s取1.0。Hoek等［2］還建立了這3個經驗參數(shù)與地質強度指標（GSI）之間的聯(lián)系

式中：mi為不同巖體的經驗參數(shù)；D為爆破影響和應力釋放對節(jié)理巖體擾動程度的參數(shù)，取值范圍0～1，當巖體沒有受到外界擾動影響時，D=0。

考慮主應力之間大小的關系，基于廣義H-B 準則的屈服函數(shù)在主應力空間可以表示為

為了在NOSBPD 中建立基于H-B 準則的彈塑性本構模型，通過式（6）求解柯西應力張量，然后基于柯西應力張量的特征矩陣將三維空間下的應力張量轉化為主應力張量，從而實現(xiàn)NOSBPD模型中主應力空間的H-B準則的描述。

2.2 主應力空間的塑性流動法則

通過主應力空間能夠判斷出某點應力狀態(tài)在屈服面上的位置關系，如單一屈服面、雙屈服面相交的棱線或者多屈服面相交的尖點處，如圖1所示，若應力狀態(tài)在σ1=σ2棱線上時，塑性應變的流動需要考慮左右2個屈服面的流動方向ra和rb。

圖1 屈服面在π平面的投影Fig.1 Projection of yield surface on π plane

采用相關聯(lián)流動法則，塑性勢函數(shù)的形式與屈服函數(shù)一致，則塑性應變張量的增量表達式為

式中：n為式（13）中屈服函數(shù)大于零的個數(shù)，當應力在屈服面、棱線和尖點時，分別需要1個、2個或多個Δλ，其中λ為塑性因子；gi為屈服函數(shù)；σ為主應力張量。當σ1＞σ2＞σ3時，主應變和主應力的塑性增量為

2.3 應力更新

當巖石進入屈服產生塑性應變時，采用主應力空間的返回映射方法計算塑性應變增量，使得更新后的應力落在屈服面上，可得主應力與塑性應變增量表達的非線性方程組為

由于在返回映射的過程中剪應力始終為零，且主應力方向不變，因此可依據之前求解的主應力特征矩陣將滿足屈服函數(shù)的主應力張量轉化為三維空間的應力張量σe，進而得到基于H-B 強度準則的非常規(guī)態(tài)型力態(tài)表達式為

將力態(tài)方程（19）代入式（1）即可求解運動方程。

2.4 斷裂準則

為了考慮巖石塑性損傷的影響，采用基于等效塑性應變的斷裂準則。等效塑性應變的表達式為

2.5 數(shù)值求解方法

為了求解NOSBPD運動方程，將模型離散為物質點的形式，位移邊界條件通過施加在邊界處的虛擬物質點進行實現(xiàn)，即通過邊界的虛擬物質點帶動內部物質點進行傳遞，力邊界條件可直接將外力轉化為體力施加在邊界處的物質點上。

基于物質點的控制方程，采用顯式的動態(tài)松弛方法，即通過增加阻尼項來求解準靜態(tài)問題［17］，如式（23）：

式中：C為人工阻尼。

然后，采用中心差分方法求解位移場u，具體的計算流程如圖2所示。

圖2 數(shù)值求解流程Fig.2 Computational flowchart of the model proposed

3 算例驗證

3.1 巖石試樣的壓縮模擬

為了驗證提出的巖石彈塑性非局部力學模型，基于平面應變假設，分別對邊長為1 m 的巖塊和含孔洞巖塊進行模擬分析，并結合有限元模擬結果進行對比驗證。假設巖石楊氏模量為28 GPa，泊松比為0.2，單軸抗壓強度為100 MPa，廣義Hoek-Brown準則中的參數(shù)s為1，mb為0.5，a為0.5。對不含孔洞的巖石模型施加位移邊界條件，對含圓孔的巖石模型施加力邊界條件進行模擬，具體模型及邊界條件如圖3所示。

圖3 模型及邊界條件Fig.3 Model and boundary conditions

對于不含孔洞的模型，將位移加載分為1 000個時間步，對模型上下表面均施加0.002 m 的位移壓縮荷載。選擇圖3中（-0.3，-0.3）為觀測點，該點豎直方向的壓縮應力與壓縮應變響應關系曲線如圖4所示。從圖中可以看出，當應變接近0.003 5時，巖石進入塑性屈服階段，本文提出的NOSBPD彈塑性模型計算結果與有限元結果完全一致。

圖4 觀測點在垂直方向的應力-應變關系Fig.4 Stress-strain relationship of observation points in vertical direction

對于含圓孔巖石的數(shù)值模型，為了精確刻畫孔洞形狀，模型離散為40 090個物質點，物質點間距為0.005 m，對模型四周直接施加80 MPa壓力進行計算。圖5給出均勻受壓荷載作用下含中心圓孔巖石模型的主應力云圖和等效塑性應變云圖。由圖可見，平面模型在四周均布壓力作用下，圓孔周邊主應力方向為沿徑向和環(huán)向的圓環(huán)形分布，本文NOSBPD模型計算得到的等效塑性應變分布與有限元結果一致，進一步驗證了所提出模型的準確性。

圖5 主應力及等效塑性應變云圖Fig.5 Contours of principal stress and equivalent plastic strain

3.2 含預制裂紋巖石的壓縮破壞模擬

巖體中預先存在的缺陷是影響巖石結構穩(wěn)定性的重要因素，本算例以文獻［18］中的受壓傾斜缺陷板試驗為對象，模擬內部裂紋的萌生與擴展過程。依據文獻［18］，試驗對象為高0.15 m、寬0.075 m的大理巖模型，其中心處的預制裂紋長0.012 7 m、傾角為30°，模型上下邊界分別施加50 MPa 的壓縮荷載。大理巖模型的楊氏模量為49 GPa，泊松比為0.19，單軸抗壓強度為84.67 MPa，GSI取95，擾動參數(shù)D為零，mi為9，通過式（12）可求得廣義H-B 準則中的參數(shù)s為0.573 8，mb為7.528 2，a為0.500 1，臨界等效塑性應變取7.5×10-6。

基于本文NOSBPD 彈塑性模型的計算結果與試驗數(shù)據對比如圖6所示?？梢姳疚哪M的裂紋擴展模式與文獻［18］中的試驗觀測基本一致，即初始的裂紋擴展路徑與預制缺陷呈90°夾角，然后隨著荷載的增加，產生翼型裂紋。此外，分別采用了11 250個、45 000個物質點離散的數(shù)值模型來模擬此問題，得到的裂紋擴展模式均與20 000個物質點離散的數(shù)值模型結果圖6a一致，驗證了本文模型處理彈塑性斷裂問題的收斂性。說明本文模型可以有效描述含裂隙巖石在受壓荷載下的翼型裂紋擴展問題。

圖6 PD預測的裂紋擴展模式與試驗結果的對比Fig.6 Comparison of PD-predicted crack growth modes with experimental results

3.3 拱形洞室的模型試驗對比分析

Vu 等［19］開展了不同埋深條件下軟巖拱形洞室模型的變形破壞試驗研究，如圖7，其中整體模型尺寸為2.00 m×2.00 m×0.40 m，洞室模型尺寸為0.37 m×0.22 m。軟巖模型的密度為2 100 kg·m-3，楊氏模量為0.3 GPa，泊松比為0.32，單軸抗壓強度為0.95 MPa，GSI為30，擾動參數(shù)D為零，mi為12，mb為0.985，s為4.2×10-4，a為0.522，臨界等效塑性應變取8×10-3［20］。限制模型前后表面的移動以模擬平面應變假設，模型底部固定，頂部和兩側壓力由單獨的液壓千斤頂系統(tǒng)控制，通過調整施加的豎向和水平應力可以模擬不同埋深和側向壓力系數(shù)。限于篇幅，本文選取側壓力系數(shù)為1 的試驗數(shù)據進行對比分析，即施加于試驗模型頂部和兩側壓力均相等，且分別考慮3 種圍壓條件，分別為280 kPa、525 kPa 和700 kPa。

圖7 拱形洞室試驗示意（單位：m）Fig.7 Sketch of vaulted cave test（unit:m）

圖8 為280 kPa 的圍壓條件下采用本文NOSBPD 模型得到的塑性區(qū)分布結果與試驗受損情況的對比分析，可見本文模型的塑性預測結果與試驗結果吻合較好。圖9給出不同圍壓條件下計算得到的模型結構損傷云圖和試驗觀測結果，可見，圍壓525 kPa條件下洞室模型基本處于完好狀態(tài)，而當圍壓增加至700 kPa 時洞室模型的底部和拱頂發(fā)生損傷破壞，整體破壞形態(tài)趨向圓形，這與圖9c的試驗觀測現(xiàn)象基本一致，從而進一步說明本文非局部模型用于模擬地下洞室彈塑性損傷分析的有效性。

圖8 塑性區(qū)域對比Fig.8 Comparison of plastic regions

圖9 不同圍壓條件下的洞室損傷Fig.9 Damage of cavern under different confining pressure conditions

4 結論

提出了一種基于Hoek-Brown（H-B）強度準則的非常規(guī)態(tài)型近場動力學（NOSBPD）彈塑性模型，可以同時描述巖石的彈塑性變形連續(xù)場特征以及斷裂力學的非連續(xù)場行為，即該模型基于H-B 準則在主應力空間的返回映射方法，能夠準確地描述巖石材料的彈塑性變形特征；同時，結合建立的等效塑性應變斷裂準則，還可以實現(xiàn)巖石在復雜荷載作用條件下的斷裂行為全過程力學分析?；谒岢龅腘OSBPD彈塑性模型，模擬分析了壓縮荷載作用下完整巖石試樣與含孔洞巖石試樣的彈塑性變形行為，計算結果與有限元結果基本一致，從而驗證了本文模型的準確性。通過與試驗觀測結果對比分析表明，本文模型還可以準確地描述在單軸壓縮條件下含預制斜裂紋的巖石試樣沿裂尖的翼型裂紋的萌生與擴展過程。此外，還將該模型應用于拱形洞室的彈塑性破壞數(shù)值模擬，發(fā)現(xiàn)當圍壓較大時，洞室的拱頂和底部首先發(fā)生破壞，分析結果與室內模型試驗觀測現(xiàn)象基本一致，進一步說明了本文模型用于巖石彈塑性斷裂問題模擬的有效性，為實際工程巖石彈塑性破壞的連續(xù)-非連續(xù)問題研究提供了可靠的分析手段。

作者貢獻聲明：

禹海濤：項目負責人、論文構思、指導模型構建及數(shù)據分析、論文修改。

胡曉錕：模型構建、數(shù)據分析及論文撰寫。

李天斌：論文修改。