禹海濤，胡曉錕，李天斌

（1. 成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610059；2. 同濟(jì)大學(xué)土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092；3. 同濟(jì)大學(xué)土木工程學(xué)院，上海 200092）

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

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

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

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

1 非常規(guī)態(tài)型近場(chǎng)動(dòng)力學(xué)理論框架

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

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

式中：J為雅可比矩陣行列式的值。根據(jù)第一類P-K應(yīng)力與應(yīng)變能和變形梯度的關(guān)系，得到近場(chǎng)動(dòng)力學(xué)的力態(tài)表達(dá)式

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

2 基于Hoek-Brown 準(zhǔn)則的近場(chǎng)動(dòng)力學(xué)彈塑性模型

2.1 屈服方程

Hoek 和Brown［15］在1980 年通過(guò)對(duì)大量試驗(yàn)進(jìn)行分析得到了H-B準(zhǔn)則，1992年Hoek等［16］對(duì)H-B準(zhǔn)則進(jìn)行改進(jìn)，提出了適用更廣泛的廣義H-B準(zhǔn)則，為

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

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

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

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

2.2 主應(yīng)力空間的塑性流動(dòng)法則

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

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

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

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

2.3 應(yīng)力更新

當(dāng)巖石進(jìn)入屈服產(chǎn)生塑性應(yīng)變時(shí)，采用主應(yīng)力空間的返回映射方法計(jì)算塑性應(yīng)變?cè)隽?，使得更新后的?yīng)力落在屈服面上，可得主應(yīng)力與塑性應(yīng)變?cè)隽勘磉_(dá)的非線性方程組為

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

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

2.4 斷裂準(zhǔn)則

為了考慮巖石塑性損傷的影響，采用基于等效塑性應(yīng)變的斷裂準(zhǔn)則。等效塑性應(yīng)變的表達(dá)式為

2.5 數(shù)值求解方法

為了求解NOSBPD運(yùn)動(dòng)方程，將模型離散為物質(zhì)點(diǎn)的形式，位移邊界條件通過(guò)施加在邊界處的虛擬物質(zhì)點(diǎn)進(jìn)行實(shí)現(xiàn)，即通過(guò)邊界的虛擬物質(zhì)點(diǎn)帶動(dòng)內(nèi)部物質(zhì)點(diǎn)進(jìn)行傳遞，力邊界條件可直接將外力轉(zhuǎn)化為體力施加在邊界處的物質(zhì)點(diǎn)上。

基于物質(zhì)點(diǎn)的控制方程，采用顯式的動(dòng)態(tài)松弛方法，即通過(guò)增加阻尼項(xiàng)來(lái)求解準(zhǔn)靜態(tài)問(wèn)題［17］，如式（23）：

式中：C為人工阻尼。

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

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

3 算例驗(yàn)證

3.1 巖石試樣的壓縮模擬

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

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

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

圖4 觀測(cè)點(diǎn)在垂直方向的應(yīng)力-應(yīng)變關(guān)系Fig.4 Stress-strain relationship of observation points in vertical direction

對(duì)于含圓孔巖石的數(shù)值模型，為了精確刻畫(huà)孔洞形狀，模型離散為40 090個(gè)物質(zhì)點(diǎn)，物質(zhì)點(diǎn)間距為0.005 m，對(duì)模型四周直接施加80 MPa壓力進(jìn)行計(jì)算。圖5給出均勻受壓荷載作用下含中心圓孔巖石模型的主應(yīng)力云圖和等效塑性應(yīng)變?cè)茍D。由圖可見(jiàn)，平面模型在四周均布?jí)毫ψ饔孟?，圓孔周邊主應(yīng)力方向?yàn)檠貜较蚝铜h(huán)向的圓環(huán)形分布，本文NOSBPD模型計(jì)算得到的等效塑性應(yīng)變分布與有限元結(jié)果一致，進(jìn)一步驗(yàn)證了所提出模型的準(zhǔn)確性。

圖5 主應(yīng)力及等效塑性應(yīng)變?cè)茍DFig.5 Contours of principal stress and equivalent plastic strain

3.2 含預(yù)制裂紋巖石的壓縮破壞模擬

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

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

圖6 PD預(yù)測(cè)的裂紋擴(kuò)展模式與試驗(yàn)結(jié)果的對(duì)比Fig.6 Comparison of PD-predicted crack growth modes with experimental results

3.3 拱形洞室的模型試驗(yàn)對(duì)比分析

Vu 等［19］開(kāi)展了不同埋深條件下軟巖拱形洞室模型的變形破壞試驗(yàn)研究，如圖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，單軸抗壓強(qiáng)度為0.95 MPa，GSI為30，擾動(dòng)參數(shù)D為零，mi為12，mb為0.985，s為4.2×10-4，a為0.522，臨界等效塑性應(yīng)變?nèi)?×10-3［20］。限制模型前后表面的移動(dòng)以模擬平面應(yīng)變假設(shè)，模型底部固定，頂部和兩側(cè)壓力由單獨(dú)的液壓千斤頂系統(tǒng)控制，通過(guò)調(diào)整施加的豎向和水平應(yīng)力可以模擬不同埋深和側(cè)向壓力系數(shù)。限于篇幅，本文選取側(cè)壓力系數(shù)為1 的試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比分析，即施加于試驗(yàn)?zāi)Ｐ晚敳亢蛢蓚?cè)壓力均相等，且分別考慮3 種圍壓條件，分別為280 kPa、525 kPa 和700 kPa。

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

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

圖8 塑性區(qū)域?qū)Ρ菷ig.8 Comparison of plastic regions

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

4 結(jié)論

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

作者貢獻(xiàn)聲明：

禹海濤：項(xiàng)目負(fù)責(zé)人、論文構(gòu)思、指導(dǎo)模型構(gòu)建及數(shù)據(jù)分析、論文修改。

胡曉錕：模型構(gòu)建、數(shù)據(jù)分析及論文撰寫(xiě)。

李天斌：論文修改。