薛大為，呂璽琳，任中俊，劉泳鋼

（1.同濟(jì)大學(xué)巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海 200092；2.同濟(jì)大學(xué)土木工程學(xué)院，上海 200092；3.湖南科技大學(xué)土木工程學(xué)院，湖南湘潭 411201；4.四川省建筑科學(xué)研究院，四川成都 610084）

自然界中多孔巖體（如砂巖）在受荷變形過程中可能發(fā)生局部化［1］，其通常表現(xiàn)為塑性（或損傷）應(yīng)變和位移場的局部化分布，并伴隨著局部化帶內(nèi)孔隙塌陷，亦可稱為“壓實(shí)局部化失穩(wěn)”［2］。壓實(shí)局部化變形模擬和預(yù)測對(duì)于地球物理和巖土工程問題具有重要意義，如斜坡、隧道的穩(wěn)定性以及天然氣回收、二氧化碳儲(chǔ)存等都可能遇到與之相關(guān)的問題。

近年來有關(guān)加載誘發(fā)的局部化變形研究方面取得較大進(jìn)展，室內(nèi)試驗(yàn)中已成功觀察到多孔巖體在壓實(shí)過程中形成的壓實(shí)帶，其發(fā)生表現(xiàn)為孔隙塌陷和顆粒破碎為特征［3-6］。Liu等［7］的研究表明，多孔巖體中可出現(xiàn)具有非零傾角的剪切增強(qiáng)壓實(shí)帶和完全垂直于最大主應(yīng)力方向的純壓實(shí)帶。當(dāng)前脆性巖土材料的延遲失穩(wěn)也逐漸成為一個(gè)熱門研究領(lǐng)域（如巖石［8-9］及多孔巖體［10］）。多孔巖體在蠕變過程中，由于其隨時(shí)間變化的特性，往往會(huì)發(fā)生以加速變形為特征的自發(fā)性失穩(wěn)［10］。多孔巖體中微結(jié)構(gòu)的存在導(dǎo)致了延遲壓實(shí)局部化的形成，且伴隨顯著的帶內(nèi)應(yīng)變加速演化現(xiàn)象。Heap等［11］研究了蠕變?cè)囼?yàn)中加速變形的時(shí)間演化和純壓實(shí)局部化帶的起始和擴(kuò)展。Brantut等［10］的試驗(yàn)研究結(jié)果表明，具有非零傾斜角度的剪切增強(qiáng)型壓實(shí)帶也能被蠕變激活。

有關(guān)多孔巖體壓實(shí)局部化的理論判別方面，通過彈塑性理論框架下的變形分叉理論或可控性分析，可預(yù)測加載導(dǎo)致的壓實(shí)局部化起始點(diǎn)［12-14］。然而，彈塑性框架不能反映多孔巖體隨時(shí)間變化的力學(xué)特性，因此需要使用彈黏塑性理論。在彈黏塑性框架下，由于無法直接定義增量應(yīng)力和增量應(yīng)變關(guān)系的剛度矩陣，導(dǎo)致分叉理論和可控性理論不再適用。Pisano等［15］通過多個(gè)常微分方程系統(tǒng)提供了一種識(shí)別彈黏塑性固體在穩(wěn)態(tài)外部攝動(dòng)（即蠕變條件）下應(yīng)變加速（即應(yīng)變不穩(wěn)定演化）和減速（即應(yīng)變穩(wěn)定演化）的可能方法。結(jié)合特定的運(yùn)動(dòng)學(xué)約束，該方法已被成功用于分析飽和及準(zhǔn)飽和黏土的延遲分散性型失穩(wěn)［16］以及多孔巖體中純壓實(shí)帶的理論判別中［14］。然而，現(xiàn)有的研究并沒有考慮具有非零傾角的剪切增強(qiáng)壓實(shí)局部化情形。

本文提出一個(gè)黏塑性壓實(shí)局部化失穩(wěn)指數(shù)，用于表征率相關(guān)多孔巖體中擬瞬時(shí)局部化和延遲局部化?；陴に苄钥煽匦岳碚?，建立非線性常微分方程系統(tǒng)，將響應(yīng)變量的加速度與響應(yīng)變量率聯(lián)系起來。推導(dǎo)該系統(tǒng)的失穩(wěn)指數(shù)，用于識(shí)別自發(fā)傳播的壓實(shí)局部化帶內(nèi)不穩(wěn)定加速變形響應(yīng)，并通過本構(gòu)模擬和有限元數(shù)值模擬對(duì)所提出失穩(wěn)指數(shù)的合理性進(jìn)行了驗(yàn)證。

1 蠕變條件下壓實(shí)局部化觸發(fā)條件

在定常外部攝動(dòng)下（即蠕變），多孔巖體可發(fā)生以應(yīng)變加速為特征的延遲壓實(shí)局部化失穩(wěn)。在這一過程中，壓實(shí)帶內(nèi)的本構(gòu)響應(yīng)可表示為

式中：及為加載應(yīng)力率和加載應(yīng)變率；相應(yīng)地，及為響應(yīng)應(yīng)變率和響應(yīng)應(yīng)力率；De為彈性矩陣；及為黏塑性應(yīng)變率。

考慮到廣義壓實(shí)帶以狹窄區(qū)域內(nèi)密集的局部變形為特征，因此可利用簡單剪切模式來描述帶內(nèi)的運(yùn)動(dòng)變形特性，如圖1所示。

圖1 壓實(shí)局部化帶內(nèi)簡單剪切運(yùn)動(dòng)學(xué)示意圖Fig.1 Schematic diagram of simple shear kinemat?ics inside compaction localization bands

從圖1可看出，帶內(nèi)變形受到法向應(yīng)力率及切向應(yīng)力率控制；與此同時(shí)，還需要滿足平面應(yīng)變條件及簡單剪切變形約束。因此，壓實(shí)帶內(nèi)的控制加載變量率和響應(yīng)變量率可表述為

式（1）中，黏塑性應(yīng)變率通過過應(yīng)力方法計(jì)算得到

式中：Φ(F)為黏性核函數(shù)，其大小取決于非彈性應(yīng)力狀態(tài)（即所謂的過應(yīng)力）與屈服面之間的距離；Q(σ)為塑性勢函數(shù)。為考慮具有任意傾角的壓實(shí)帶，需要進(jìn)一步考慮旋轉(zhuǎn)參考系統(tǒng)，旋轉(zhuǎn)角度即壓實(shí)帶的傾角，可通過方向余弦矩陣T對(duì)參考系統(tǒng)施加旋轉(zhuǎn)［16］。相應(yīng)地，旋轉(zhuǎn)應(yīng)力率、應(yīng)變率、彈性矩陣可表示為

式中

式中：θ為壓實(shí)帶傾角；K和G分別為彈性體積模量及剪切模量。

顯然，式（1）-（4）描述了多孔巖體壓實(shí)帶內(nèi)的本構(gòu)特性。利用上述公式推導(dǎo)出響應(yīng)變量率ε?α及σ?β為

式中：Iαα和Oβα分別為單位矩陣及零矩陣。進(jìn)一步地，可推導(dǎo)得到響應(yīng)變量的演化加速度

式中：和為體力項(xiàng)，在蠕變條件下為零。考慮簡單剪切條件（即=0）及定常外部攝動(dòng)條件（即蠕變=0），式（9）可轉(zhuǎn)化為

在蠕變條件下有

式中：H為硬化模量；Hχ與式（2）、式（3）中的控制加載變量及響應(yīng)變量的具體形式相關(guān)，亦可稱為可控性模量。

最終，控制響應(yīng)變量變速演化規(guī)律的常微分方程系統(tǒng)得以建立。式（10）通過依賴時(shí)間變化的矩陣Z（即彈黏塑性本構(gòu)算子）控制響應(yīng)變量的變速演化規(guī)律。若矩陣Z的特征值為負(fù)，則響應(yīng)變量將呈現(xiàn)減速演變，可視為穩(wěn)定響應(yīng)；然而若矩陣Z的特征值為正，將導(dǎo)致響應(yīng)變量速率的加速演變，可視為不穩(wěn)定響應(yīng)。式（10）表明，矩陣Z的特征值由Zαα和Zββ的特征值組成。Zαα為對(duì)角矩陣且特征值的正負(fù)號(hào)由H in控制；而Zββ為一個(gè)對(duì)角矩陣和一個(gè)海森矩陣之和。由于塑性勢函數(shù)為凸函數(shù)，海森矩陣的特征值為恒正，因此Zββ的特征值永遠(yuǎn)小于Zαα的特征值。這就意味著當(dāng)Hin＞0時(shí)，響應(yīng)變量呈現(xiàn)減速演化趨勢，即該系統(tǒng)是恒穩(wěn)定的；相反，當(dāng)Hin＜0時(shí)，響應(yīng)變量速率有加速發(fā)展的潛在可能，意味著系統(tǒng)是潛在不穩(wěn)定的。因此，Hin＞0是彈黏塑性材料穩(wěn)定的充分條件，而Hin＜0則對(duì)應(yīng)潛在壓實(shí)局部化失穩(wěn)。

2 多孔巖體壓實(shí)局部化理論預(yù)測

2.1 多孔巖體彈黏塑性模型

將上述理論與一種應(yīng)變硬化彈黏塑性本構(gòu)模型相結(jié)合，該模型由Nova等［17］提出，包含了硬化和軟化的相互競爭機(jī)制。屈服函數(shù)（當(dāng)h=F時(shí)，代表屈服函數(shù)）及塑性勢函數(shù)（當(dāng)h=Q時(shí)，代表塑性勢函數(shù)）表示為

通過標(biāo)定參數(shù)M h、μh及αh，可靈活控制式（13）中給出的屈服面和塑性勢面的形狀。模型中的硬化變量率和可表達(dá)為

式中：Bp、ρm及ξm為材料參數(shù)。ps可用于捕捉由于孔隙變化導(dǎo)致的硬化或軟化機(jī)制，其中塑性體積應(yīng)變反映孔隙的變化。多孔巖體壓縮變形時(shí)，顆?？紫稌?huì)明顯減小，這一物理過程可由ps反映，而pm則同時(shí)通過體積塑性應(yīng)變和等效塑性剪應(yīng)變引入軟化機(jī)制。

利用式（3），該模型可轉(zhuǎn)變?yōu)閺楌に苄阅Ｐ?，黏塑性?yīng)變可通過式（16）求得：

式中：pc0為pc的初始值；ω為黏性參數(shù)。模型的彈性部分采用線彈性模型。

根據(jù)上述模型，Marinelli和Buscarnera［18］基于三軸試驗(yàn)結(jié)果及變形分叉理論，標(biāo)定了多種多孔巖體材料參數(shù)，標(biāo)定的Berea砂巖材料參數(shù)如表1所示?；谒⒌谋緲?gòu)模型開展有限元數(shù)值模擬，在不同圍壓（p0=75 MPa、150 MPa、200 MPa以及300 MPa）下模擬得到的力學(xué)響應(yīng)與三軸試驗(yàn)結(jié)果對(duì)比如圖2所示，兩者較為符合。

表1 Berea砂巖材料參數(shù)［18］Tab.1 Material parameters of Berea sandstone[18]

圖2 Berea砂巖力學(xué)響應(yīng)Fig.2 Mechanical response of Berea sandstone

2.2 平面應(yīng)變單元壓實(shí)局部化分析

進(jìn)一步地，對(duì)平面應(yīng)變狀態(tài)下延遲壓實(shí)局部化進(jìn)行分析。將平面應(yīng)變單元視為一個(gè)積分點(diǎn)，驗(yàn)證所提出的黏塑性壓實(shí)局部化失穩(wěn)指數(shù)的適用性。分析包括：①各向同性壓縮，使模型具有初始圍壓；②平面應(yīng)變剪切模擬，使模型進(jìn)入非彈性狀態(tài)；③改變約束條件，使模型進(jìn)入蠕變變形狀態(tài)。對(duì)于分析步②，當(dāng)應(yīng)力狀態(tài)達(dá)到初始屈服面時(shí)，對(duì)應(yīng)于某個(gè)角度的失穩(wěn)指數(shù)可能變?yōu)榱隳酥霖?fù)值（即Hin(θ)≤0），此角度即為可能發(fā)生的壓實(shí)局部化帶傾角。當(dāng)存在Hin(θ)≤0時(shí)，Hin(θ)最小值對(duì)應(yīng)的角度為最可能出現(xiàn)的局部化帶傾角，即A(θmin)=minA(θ)。在各種可能的局部化變形中，A(θ=0)＜0時(shí)發(fā)生的變形稱為壓實(shí)局部化變形。這種情況下，傾斜角度為0°的局部化帶成為一種潛在可能。這種壓實(shí)局部化變形又分為2種情況，其一為A(θmin)＜A(θ=0)，具有非0°傾斜角的剪切增強(qiáng)壓實(shí)帶；其二為A(θmin)＜A(θ=0)，具有0°傾斜角的純壓實(shí)帶。在隨后的分析步③中，將約束條件更改為與該傾角對(duì)應(yīng)方向施加簡單剪切蠕變變形。應(yīng)注意到，分析步③所施加的控制條件體現(xiàn)了壓實(shí)變形帶內(nèi)的變形運(yùn)動(dòng)特征。

為研究當(dāng)應(yīng)力路徑達(dá)到不同屈服面區(qū)域后的模型響應(yīng)及失穩(wěn)指數(shù)變化規(guī)律，選擇p1=105 MPa、p2=190 MPa這2種圍壓開展Berea砂巖模擬，計(jì)算結(jié)果如圖3所示。應(yīng)力路徑與屈服面的交點(diǎn)分別為屈服點(diǎn)Ⅰ和屈服點(diǎn)Ⅱ，由于黏塑性理論允許應(yīng)力狀態(tài)超越屈服面，采用圍壓p1和p2計(jì)算所得應(yīng)力狀態(tài)與屈服面變化的不匹配造成了變形特性的差異。當(dāng)應(yīng)力狀態(tài)達(dá)到屈服面時(shí)，應(yīng)力狀態(tài)增長速度大于屈服面擴(kuò)展速度，此階段為不穩(wěn)定蠕變變形階段，可產(chǎn)生由加速變形造成的蠕變壓實(shí)局部化失穩(wěn)；而當(dāng)應(yīng)力狀態(tài)增長速度逐漸小于屈服面擴(kuò)展速度時(shí)，不穩(wěn)定蠕變響應(yīng)轉(zhuǎn)化為穩(wěn)定蠕變響應(yīng)，此時(shí)出現(xiàn)的壓實(shí)局部化為穩(wěn)定變形。針對(duì)屈服點(diǎn)Ⅰ和Ⅱ的所有傾角失穩(wěn)指數(shù)計(jì)算結(jié)果如圖4所示，從圖中可看出，基于所提出失穩(wěn)指數(shù)的壓實(shí)局部化理論預(yù)測與基于傳統(tǒng)分叉理論的預(yù)測結(jié)果相匹配。這2種理論所得結(jié)果表明：屈服點(diǎn)Ⅰ點(diǎn)的應(yīng)力狀態(tài)對(duì)應(yīng)傾角為±21°的剪切增強(qiáng)壓實(shí)帶；屈服點(diǎn)Ⅱ點(diǎn)預(yù)測得到傾角為0°的純壓實(shí)帶。蠕變階段軸向應(yīng)變的演化如圖5a所示，軸向蠕變應(yīng)變最開始呈現(xiàn)為加速不穩(wěn)定發(fā)展，隨后應(yīng)變率逐漸降低至趨于零，軸向應(yīng)變達(dá)到穩(wěn)定值。根據(jù)圖3結(jié)果可知，應(yīng)力狀態(tài)和屈服面演化速率不匹配造成了應(yīng)變加速（或減速）演化，因而應(yīng)力狀態(tài)增長速度等于屈服面擴(kuò)展速度時(shí)對(duì)應(yīng)的時(shí)間點(diǎn)即為從不穩(wěn)定蠕變到穩(wěn)定蠕變的轉(zhuǎn)折點(diǎn)（實(shí)心點(diǎn)標(biāo)出）。圖5b表明，圍壓p1和p2條件下的失穩(wěn)指數(shù)最初為負(fù)值，亦即式（8）中彈黏塑性算子的特征值為正，表明軸向應(yīng)變呈加速變化趨勢，即發(fā)生不穩(wěn)定蠕變。隨著蠕變變形發(fā)展，特征值由負(fù)值轉(zhuǎn)化為正值，即逐漸進(jìn)入減速變形階段。圖5b中失穩(wěn)指數(shù)正負(fù)號(hào)的轉(zhuǎn)折點(diǎn)與圖5a不穩(wěn)定蠕變轉(zhuǎn)折點(diǎn)一致，說明所提出的黏塑性壓實(shí)局部化失穩(wěn)指數(shù)能有效識(shí)別蠕變的加速失穩(wěn)和減速穩(wěn)定階段。

圖3 材料點(diǎn)應(yīng)力路徑分析Fig.3 Stress paths of material point analysis

圖4 分叉理論與失穩(wěn)指數(shù)計(jì)算結(jié)果對(duì)比Fig.4 Comparison of the bifurcation theory and the instability index

圖5 Berea砂巖剪切蠕變分析Fig.5 Analysis of simple shear creep of Berea sandstone

3 多孔巖體壓實(shí)局部化有限元模擬

進(jìn)一步開展Berea砂巖局部化形成過程的平面應(yīng)變有限元數(shù)值模擬。分析模型及網(wǎng)格如圖6a所示，模型長寬比為2：1，由12 800個(gè)C3D8單元組成。底部邊界施加全約束，頂部按應(yīng)變速率10-5s-1施加荷載，以使區(qū)域內(nèi)每個(gè)高斯點(diǎn)的應(yīng)力狀態(tài)達(dá)到屈服面，然后將頂部的牽引力保持恒定，以確保試樣可以在固定的外部擾動(dòng)下隨時(shí)間蠕變。值得注意的是，這里有限元計(jì)算模擬所采用的參數(shù)、軸壓、控制條件以及加載量與第3節(jié)點(diǎn)積分理論分析完全相同。為了觸發(fā)應(yīng)變局部化，在模型中間設(shè)置了一個(gè)薄弱區(qū)域（將pc取為原值的95%）。

Berea砂巖有限元數(shù)值模擬結(jié)果如圖6b。由圖可見，一旦應(yīng)力狀態(tài)進(jìn)入塑性區(qū)（T/Tc=0），在圍壓p1和p2條件下計(jì)算得到的塑性體積應(yīng)變場呈壓實(shí)局部化變形模式，且得到的局部化帶傾角與材料點(diǎn)預(yù)測結(jié)果完全一致。圖6b還體現(xiàn)了恒定軸向應(yīng)力作用下體積塑性應(yīng)變的空間傳播模式。在圍壓p1和p2條件下所得的體積塑性應(yīng)變首先從薄弱區(qū)域出現(xiàn)，然后在模型內(nèi)部形成多條壓實(shí)帶，這種壓實(shí)帶的形成意味著可能發(fā)生局部化失穩(wěn)。從位于壓實(shí)局部化帶內(nèi)的高斯積分點(diǎn)獲得的局部響應(yīng)如圖7所示，圍壓p1和p2計(jì)算結(jié)果中提取的高斯點(diǎn)（圖6中灰色星號(hào)示出）顯示出最初的加速趨勢（即不穩(wěn)定）以及隨后的減速趨勢（即穩(wěn)定）。塑性體應(yīng)變區(qū)域內(nèi)的穩(wěn)定和不穩(wěn)定應(yīng)變響應(yīng)分別與對(duì)應(yīng)傾角的失穩(wěn)指數(shù)的正負(fù)號(hào)有關(guān)（即傾角為21°和0°分別對(duì)應(yīng)Hin(θ=21o)和Hin(θ=0o)）。分析結(jié)果表明，所提出的黏塑性壓實(shí)局部化失穩(wěn)指數(shù)能很好地用于率相關(guān)多孔巖體中壓實(shí)局部化的分析。

圖6 Berea砂巖平面應(yīng)變壓縮試驗(yàn)?zāi)MFig.6 Numerical simulations of plane strain compression test for Berea sandstone

圖7 Berea砂巖有限元高斯積分點(diǎn)分析Fig.7 Finite element Gaussian point analysis of Berea sandstone

4 結(jié)語

基于黏塑性假設(shè)及壓實(shí)帶內(nèi)變形運(yùn)動(dòng)學(xué)特性，建立了描述多孔巖體蠕變條件下應(yīng)變局部化響應(yīng)變量隨時(shí)間演化的常微分方程組，并提出了一個(gè)判別多孔巖體蠕變誘發(fā)延遲局部化的黏塑性壓實(shí)局部化失穩(wěn)指數(shù)。進(jìn)一步地，采用能捕捉不同應(yīng)變局部化模式的彈黏塑性應(yīng)變硬化及服從非關(guān)聯(lián)流動(dòng)法則的本構(gòu)模型，開展了材料點(diǎn)失穩(wěn)分析及有限元模擬，結(jié)果表明：①通過建立的理論，在加載階段不同壓實(shí)應(yīng)變局部化模式的預(yù)測結(jié)果與率無關(guān)彈塑性變形分叉理論預(yù)測結(jié)果一致；②蠕變過程模擬中響應(yīng)變量的加速和減速是由超過彈性域的應(yīng)力狀態(tài)與屈服面之間的相對(duì)運(yùn)動(dòng)引起的，而失穩(wěn)指數(shù)的符號(hào)變化（由負(fù)到正）對(duì)應(yīng)了響應(yīng)變量由加速向減速演化的轉(zhuǎn)折點(diǎn)。這些結(jié)果表明，文中提出的失穩(wěn)指數(shù)能夠合理判別多孔巖體蠕變壓實(shí)局部化現(xiàn)象。

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

薛大為：理論推導(dǎo)、程序編寫及稿件撰寫。

呂璽琳：提出研究思路和方法、稿件審核。

任中?。航Y(jié)果校核、稿件撰寫。

劉泳鋼：方案設(shè)計(jì)、稿件修訂。