基于微平面模型和正則化的脆性巖石損傷破裂模擬方法

袁 陽， 徐 濤， 周廣磊， 勒治華

(1. 東北大學(xué) 資源與土木工程學(xué)院， 遼寧 沈陽 110819； 2. 山東科技大學(xué) 能源與礦業(yè)工程學(xué)院， 山東 青島 266590)

在外力和環(huán)境的作用下，由于細觀結(jié)構(gòu)缺陷(如微裂紋、微孔洞等)引起的材料或結(jié)構(gòu)的劣化過程稱為損傷.根據(jù)所采用的損傷力學(xué)理論，巖石損傷本構(gòu)模型可分為宏觀損傷模型和細觀損傷模型[1].針對巖石、混凝土等準脆性材料非彈性力學(xué)行為與微裂紋張開、閉合、摩擦滑移等細觀力學(xué)損傷機制，Bazant等對塑性滑移理論[2]進行拓展，建立了微平面模型(microplane model)[3-4]，其基本思想是直接以應(yīng)力和應(yīng)變的矢量而不是張量建立本構(gòu)關(guān)系.經(jīng)過不斷發(fā)展，微平面模型已經(jīng)成功用于混凝土、土和巖石的本構(gòu)模型中[5].微平面是指在宏觀連續(xù)介質(zhì)體內(nèi)任何一點的無窮小區(qū)域內(nèi)垂直于任意方向的平面.與傳統(tǒng)方法采用宏觀連續(xù)介質(zhì)應(yīng)力、應(yīng)變張量及其不變性來表達材料的本構(gòu)模型不同，微平面模型通過在給定的材料點上作用于所有可能取向的單個微平面上的應(yīng)力和應(yīng)變矢量來描述材料在宏觀尺度上的行為.這種從微觀到宏觀的分析方法為描述巖石的力學(xué)本構(gòu)行為提供了新思路[6].

數(shù)值模擬方法在巖石破壞機理研究中獲得廣泛應(yīng)用，而采用連續(xù)介質(zhì)損傷力學(xué)方法在模擬巖石類材料破壞時會導(dǎo)致局部化也已成為共識[7-8]，這意味著有限元模擬分析中使用應(yīng)變軟化模型時會嚴重依賴于離散化網(wǎng)格的尺寸[9].Mondal等[10]認為這種全局響應(yīng)對空間離散化表現(xiàn)出很強的依賴性的原因是當應(yīng)用網(wǎng)格細化連續(xù)介質(zhì)時，漸進損傷破壞區(qū)域的解不能收斂.崔煥平等[11]將損傷判定的極限拉應(yīng)變值替換成與網(wǎng)格尺寸相關(guān)的函數(shù)來解決網(wǎng)格尺寸效應(yīng).Bazant等[12]將非局部理論用于分析準脆性材料損傷問題，初步發(fā)展了非局部損傷力學(xué)方法，并對積分型非局部模型進行了綜述[13].基于準脆性材料在壓縮和剪切狀態(tài)下拉伸開裂和非線性三軸響應(yīng)的微平面模型，Bazant等[14]提出了一種能夠反映材料損傷和結(jié)構(gòu)尺寸效應(yīng)的非局部微平面模型.文獻[15-16]基于ABAQUS有限元軟件二次開發(fā)平臺將非局部理論[17]嵌入M7微平面模型[18]中，對混凝土應(yīng)變局部化的有限元網(wǎng)格敏感性進行了分析.但目前提出的微平面模型及其損傷均是通過在每個微平面上添加指數(shù)函數(shù)來表征微平面應(yīng)力到達閾值后產(chǎn)生的應(yīng)力軟化從而影響宏觀應(yīng)力[16].

本文提出了一種基于有效應(yīng)力概念的微平面損傷模型，通過引入非均質(zhì)細觀模型的思想[19]，摒棄了微平面中采用的應(yīng)力邊界假設(shè)，大大減少所需的力學(xué)參數(shù).基于有限元軟件COMSOL及MATLAB二次開發(fā)平臺，將微平面非線性本構(gòu)關(guān)系所采用的PDE(partial differential equation)形式嵌入程序中，編寫了定義微平面損傷的相關(guān)程序，并在數(shù)值模型中引入高階梯度正則化方法，通過在方程中定義材料的內(nèi)部特征長度來防止局部化問題，模型中節(jié)點處的應(yīng)力或應(yīng)變狀態(tài)不僅受到該點歷史應(yīng)力應(yīng)變狀態(tài)影響，而且還受到鄰近節(jié)點狀態(tài)的影響，一定程度上解決了有限元計算中的網(wǎng)格依賴性問題.

1 數(shù)值模型與方法

1.1 微平面模型

微平面模型通過應(yīng)力率和應(yīng)變率之間的關(guān)系構(gòu)建微平面上本構(gòu)關(guān)系來描述宏觀應(yīng)力應(yīng)變與時間之間的變化關(guān)系.其基本假設(shè)有靜態(tài)約束和動態(tài)約束，動態(tài)約束指微平面上的應(yīng)變等于宏觀應(yīng)變張量εij的投影分量，其宏觀應(yīng)變張量與微平面上的應(yīng)變分量直接相關(guān)，而靜態(tài)約束則是指微平面上的應(yīng)力等于宏觀應(yīng)力張量σi,j的投影分量，其宏觀應(yīng)力張量與微平面上的應(yīng)力分量直接相關(guān)[20].本文采用的動態(tài)約束適用于描述微平面應(yīng)變沿各方向連續(xù)變化，而微平面應(yīng)力沿各方向變化不連續(xù)的情況，微平面上法向、切向應(yīng)變分量與宏觀應(yīng)變張量之間的關(guān)系可表示為

(1)

式中:ni為微平面空間法向向量；mi,li為切向向量分量，i，j=1,2,3；εN,εM,εL(N為法向，M, L為切向T分解的分量)分別為動態(tài)約束條件下微平面上的法向和切向應(yīng)變，其在微平面上的關(guān)系如圖1所示；Nij,Mij,Lij具有對稱性，其定義如下：

(2)

圖1為2×21個積分點微平面模型，其空間形狀為一正二十面體，積分點(微平面)的位置為其頂點(12個)和棱中點(30個)[21].以應(yīng)變率和應(yīng)力率形式表示的微平面本構(gòu)關(guān)系(式(3))能夠很好地描述模型應(yīng)力應(yīng)變與時間的聯(lián)系[22]，該關(guān)系在COMSOL中通過添加PDE物理場來實現(xiàn).

(3)

其中，EN,ET分別為法向和切向彈性常數(shù)，其與單元的彈性模量關(guān)系為

(4)

式中，υ為單元節(jié)點的泊松比.

圖1 微平面上應(yīng)變率矢量的分量

宏觀應(yīng)力張量為微平面應(yīng)力在各積分點上不同方向進行積分求值之和，為在程序中更好地進行計算，宏觀應(yīng)力張量與微平面應(yīng)力矢量通過不同方向積分后的近似關(guān)系如式(5)所示，其積分點的空間位置及權(quán)重詳細論述可見文獻[6].

(5)

式中：n代表總的積分點數(shù)目;k表示某一微平面；ωk為微平面k的權(quán)重.

1.2 微平面損傷模型

基于應(yīng)變等效假設(shè)[23]的損傷在細觀損傷力學(xué)模型中已得到廣泛的應(yīng)用[24-26]，其基本思想是將無損材料中的名義應(yīng)力以有損材料的有效應(yīng)力代替建立的損傷本構(gòu)關(guān)系.有損材料的彈性模量定義為

E=E0(1-D) .

(6)

式中:E0為材料的初始模量;D為損傷值(0≤D≤1).

本文提出的微平面損傷模型的應(yīng)力狀態(tài)基于屈服準則判定.首先，對于處于拉應(yīng)力狀態(tài)時的節(jié)點，即法向應(yīng)力為正(σN> 0)的情況采用最大拉應(yīng)力準則，認為某個微平面上的最大法向應(yīng)力超過閾值時，該微平面系統(tǒng)對應(yīng)的節(jié)點即發(fā)生損傷，可表示為

(7)

當單元處于壓縮狀態(tài)(σN≤0)時，采用Mohr-Coulomb(MC)屈服準則.在微平面系統(tǒng)中可以理解為當單元的某個微平面上的正應(yīng)力和剪應(yīng)力滿足該條件時發(fā)生剪切破壞，可表示為

(8)

D=D(n) .

(9)

對于單個微平面而言，損傷仍然是一個標量.對于某一單元節(jié)點，發(fā)生拉伸損傷時，其發(fā)生損傷的微平面方向的值可以表示為

(10)

其中：εN為微平面法向應(yīng)變；εt0為開始發(fā)生損傷的拉伸應(yīng)變值；η為殘余強度系數(shù)；εtu為該微平面對應(yīng)的極限拉伸應(yīng)變.發(fā)生剪切損傷時，損傷值可表示為

(11)

其中：εc0為微平面法向應(yīng)變εN小于零時對應(yīng)的最先發(fā)生剪切破壞的微平面所對應(yīng)的壓縮應(yīng)變閾值；σc0為最先發(fā)生剪切破壞的微平面上的極限壓縮應(yīng)力值.為了區(qū)別不同損傷值，對滿足最大拉應(yīng)力破壞準則的拉伸損傷規(guī)定為負值，對滿足MC屈服準則的剪切損傷規(guī)定為正值.

1.3 正則化方法引入

(12)

式中,?n為n階梯度.又因為

(13)

式中:?2為拉普拉斯算子;系數(shù)α,β為高階項系數(shù).由式(13)可得

(14)

將式(14)代入式(12)并略去高階項，可得

(15)

1.4 巖石的非均質(zhì)表征

巖石由各種晶粒、膠結(jié)物和孔隙組成，為表征其非均質(zhì)特性，假設(shè)其內(nèi)部細觀參數(shù)服從Weibull分布：

(16)

式中:u為巖石細觀參數(shù)(彈性模量、強度、內(nèi)聚力等);u0是一個與單元參數(shù)平均值有關(guān)的參數(shù);m為非均質(zhì)系數(shù).從圖2可知，隨著非均質(zhì)系數(shù)的增大，細觀參數(shù)的分布范圍變得愈加狹窄，表示巖石內(nèi)部參數(shù)愈均勻.將Weibull分布函數(shù)引入數(shù)值模型中，其物理意義在于表征巖石的非均質(zhì)性，在數(shù)值模型中可以防止大量單元在同一加載步出現(xiàn)損傷，其次，此方法可以體現(xiàn)出材料的“宏觀塑性”[25].

圖2 非均質(zhì)系數(shù)m對應(yīng)的概率密度分布(彈性模量)

2 模型在COMSOL中的實現(xiàn)及驗證

微平面模型需要針對每一個微平面構(gòu)建本構(gòu)關(guān)系，本文模型基于COMSOL Multiphysics 5.2 和MATLAB平臺進行二次開發(fā).

2.1 模型在COMSOL平臺實現(xiàn)過程

通過在COMSOL中建立微平面損傷模型，在MATLAB平臺中加入循環(huán)函數(shù)以實現(xiàn)巖石參數(shù)(彈性模量、強度、內(nèi)聚力和損傷)的每一步求解后的更新迭代.具體流程(見圖3)為，先對每一個方向上的微平面本構(gòu)關(guān)系在COMSOL界面中通過添加PDE物理場以建立微平面上的應(yīng)力-應(yīng)變隨時間變化的關(guān)系，如需建立2×21個積分點的微平面模型，則需要建立21個物理場，每個物理場中包含3個偏微分方程對應(yīng)式(3)，對于本文中的數(shù)值模型，表征微平面模型共建立63個偏微分方程.計算流程為首先對每個單元節(jié)點上的局部宏觀應(yīng)變增量根據(jù)式(1)進行分解，即微平面k上的微平面法向應(yīng)變矢量為非局部應(yīng)變張量與法向投影張量的乘積(動態(tài)約束).對于每一個應(yīng)變分量，可以建立Helmholtz方程通過式(15)進行局部應(yīng)變和非局部應(yīng)變之間的替換.值得注意的是，本文中采用正則化模型對三點彎曲實驗進行模擬，該實驗預(yù)制裂縫處主要為拉伸破壞.因此為減少計算量，本文只對微平面法向應(yīng)變進行非局部變量的替換，而切向應(yīng)變則不替換.得到微平面非局部應(yīng)變矢量分量后，根據(jù)微平面本構(gòu)關(guān)系式(3)進行微平面應(yīng)力的計算，繼而根據(jù)微平面上的損傷狀態(tài)進行判斷是否需要對單元的彈性模量進行折減，最后通過高斯公式進行數(shù)值積分得到該加載步下的宏觀應(yīng)力張量.

2.2 模型驗證

根據(jù)上述步驟，編寫基于微平面損傷模型的MATLAB程序，為節(jié)省計算資源，選取單位球面積分點數(shù)目為2×21.實驗建立了單軸壓縮模型，試件尺寸為直徑25 mm，高62.5 mm.數(shù)值模型和花崗巖參數(shù)如表1所示.值得注意的是，數(shù)值模擬的彈性模量和強度是單元的平均值[27]，而實驗中的彈性模量和強度是實驗結(jié)果.

圖4繪制了數(shù)值模型和實驗的應(yīng)力-應(yīng)變曲線，實驗開始加載階段為孔隙壓密階段，后線性部分為彈性階段，接近峰值時出現(xiàn)非線性階段，試件達到峰值強度162.3 MPa后破壞，由于實驗過程的應(yīng)變通過應(yīng)變片來采集，故不能顯示峰后階段.

圖3 基于微平面損傷模型的正則化方法計算流程圖

表1 單軸壓縮實驗和模擬的材料參數(shù)

本文提出的微平面損傷模型并未假設(shè)孔隙壓密階段，因此不能表征這一階段.圖4中的點劃線為數(shù)值模擬線性階段平移曲線，可知對于彈性階段和接近峰值的非線性階段以及峰值強度，數(shù)值結(jié)果與實驗結(jié)果較為吻合.

限于篇幅，圖5只給出了單軸壓縮數(shù)值模擬結(jié)果的部分損傷演化過程圖，圖例中不同數(shù)值代表損傷值的種類和大小.根據(jù)1.2節(jié)定義正值為壓縮狀態(tài)下的剪切損傷，負值為拉伸損傷，絕對值越大代表該種損傷越大.在開始加載階段，模型損傷為零.隨著加載步增大，試件隨機出現(xiàn)了剪切和拉伸損傷，在短時間內(nèi)，這些損傷逐漸匯聚貫通，最后與單軸壓縮實驗中的試件相同形成剪切帶.而數(shù)值模擬在第66～68步恰好為峰值強度后出現(xiàn)應(yīng)力跌落.結(jié)果表明，模型能夠較好地反映巖石在應(yīng)力加載條件下的損傷演化過程.

圖4 單軸壓縮數(shù)值模擬和實驗結(jié)果應(yīng)力-應(yīng)變曲線

圖5 單軸壓縮數(shù)值模擬和實驗破壞結(jié)果對比

3 正則化模型分析

采用微平面損傷模型(本文稱之為普通模型)得到的宏觀變形、應(yīng)力及損傷發(fā)展依賴于有限元網(wǎng)格尺寸和數(shù)量.本文采用三點彎曲試件進行正則化方法分析，試件尺寸如圖6所示，底部2個支座固定，試件跨度為70 mm，試件長度為90 mm，試件高度為30 mm，寬度為10 mm，預(yù)制裂縫a為寬1 mm，高4 mm.在數(shù)值模擬中，倘若常采用圓柱作為上下支點，其與試件為線接觸，在進行損傷迭代計算時，易于在支點處產(chǎn)生應(yīng)力集中發(fā)生破壞，與實際情況不符.為獲得良好的計算結(jié)果，本文將圓柱體與試件接觸改為寬2 mm，長10 mm的矩形接觸面，加載速率為1×10-6m/s.

3.1 非均質(zhì)系數(shù)對損傷的影響

首先對不同非均質(zhì)系數(shù)(1，5和10)的微平面損傷普通模型在不同時間加載步下的破壞模式

圖6 三點彎曲實驗試件尺寸(單位：mm)

進行對比分析.除加載速率外數(shù)值模型參數(shù)見表1.為提高模型計算速度并防止網(wǎng)格收斂性問題，模型網(wǎng)格劃分采用COMSOL默認的四面體網(wǎng)格(Mesh1)，最大單元尺寸為3.15×10-3m，最小單元尺寸為1.35×10-4m，最大單元增長率為1.35，單元總數(shù)為70 109個.圖7給出了不同非均質(zhì)系數(shù)m的試件損傷演化.由圖7可見，非均質(zhì)系數(shù)對試件的損傷破壞影響明顯，非均質(zhì)系數(shù)低(m=1)時，巖樣中損傷在預(yù)置裂縫周圍彌散分布；非均質(zhì)系數(shù)高(5，10)時，巖樣中損傷集中于預(yù)置裂縫尖端，形成明顯的局部化損傷區(qū)域.非均質(zhì)系數(shù)為1時，損傷區(qū)域雖然在預(yù)制裂縫處較多，但是在試件其他區(qū)域也有大量損傷區(qū)域存在，與實際情況不符.對于巖石材料而言，非均質(zhì)系數(shù)為5左右時結(jié)果較為合理[25-27].

3.2 特征長度對損傷影響

根據(jù)有限元模型網(wǎng)格尺寸(1.35×10-4～3.15×10-3m)，采用特征長度分別為0.3，0.4，0.5，1，2 mm進行數(shù)值模擬，數(shù)值模型非均質(zhì)系數(shù)均為5，網(wǎng)格劃分方式(Mesh1)是由COMSOL根據(jù)物理場自動生成的四面體形式.圖8為具有不同特征長度的數(shù)值模型的損傷云圖(時步為50)，圖例中的數(shù)值其意義與2.2節(jié)中相同.從圖8可知，采用相同網(wǎng)格劃分方式，正則化模型與普通模型的破壞模式大致相同.因此正則化模型可以看作是對普通微平面損傷模型中某些變量(如應(yīng)變)的修正，并不會導(dǎo)致結(jié)果產(chǎn)生巨大差異.并且，在相同網(wǎng)格劃分下，特征長度越小，越接近最小網(wǎng)格尺寸，損傷效果越好；特征長度為0.4 mm和 0.3 mm 時損傷位置和損傷發(fā)展方向準確，但最終損傷效果與特征長度為0.5 mm的結(jié)果相近，因此模型單元尺寸與特征長度相近時，可認為對損傷結(jié)果影響較小.特征長度(l=2 mm)過大時，由于特征長度大于大部分網(wǎng)格尺寸，導(dǎo)致?lián)p傷面積過大而不符合實際情況.由此可見，在相同網(wǎng)格劃分條件下，存在一個最優(yōu)特征長度值，可以得到最好的裂紋擴展效果.在可預(yù)見的范圍內(nèi)，對相同尺寸的數(shù)值模型而言，隨著網(wǎng)格密度增大，最小單元尺寸減小，最優(yōu)特征長度值也會隨之減小.

圖7 不同非均質(zhì)系數(shù)的試件損傷演化

圖8 不同特征長度的損傷演化

3.3 網(wǎng)格敏感性分析

采用第二種較稀疏的網(wǎng)格(Mesh2)進行數(shù)值模擬，均質(zhì)度為5.圖9為普通和正則化模型三點彎曲數(shù)值模型在兩種網(wǎng)格情況下第50步時局部損傷云圖，圖中Mesh1網(wǎng)格的最大單元為3.15×10-3m，最小單元為1.35×10-4m，總單元數(shù)量為70 109個，Mesh2網(wǎng)格的最大單元尺寸為7.2×10-3m，最小單元尺寸為9×10-4m，總單元數(shù)量為10 772個.圖中正則化模型采用的特征長度為1 mm.可以得知，在改變網(wǎng)格密度后，在相同加載步數(shù)下，普通微平面模型的破壞模式發(fā)生了較大的改變，損傷區(qū)域長度變短，并且與裂縫平行的相鄰單元也發(fā)生了損傷.而對于正則化模型而言，改變網(wǎng)格密度后，其破壞區(qū)域大小和方向基本保持不變，說明引入正則化方法后，模型損傷結(jié)果對網(wǎng)格依賴性較小.

圖9 普通模型和正則化模型兩種網(wǎng)格密度損傷演化

對兩種網(wǎng)格劃分Mesh1和Mesh2的載荷-位移曲線進行普通模型和正則化模型的對比分析.載荷為施加在試件上部區(qū)域的z方向反力，位移為施加載荷的區(qū)域z方向位移.圖10為普通和正則化模型在兩種不同網(wǎng)格密度下對三點彎曲試件進行數(shù)值模擬，其中正則化模型的特征長度為1 mm，非均質(zhì)系數(shù)為5.由圖10可知，隨著網(wǎng)格密度減小(網(wǎng)格尺寸增大)，普通模型模擬的強度發(fā)生了明顯變化，此為有限元計算過程中出現(xiàn)的相同模型由于網(wǎng)格尺寸大小、密度不同產(chǎn)生的尺寸效應(yīng).總體而言，模型在引入正則化方法后，試件的整體強度和應(yīng)變峰值增大，但不同網(wǎng)格密度大小對正則化模型的載荷-位移曲線影響較小，二者趨勢在第40計算步之前基本一致，這是由于正則化微平面損傷模型可以考慮單元與有限范圍內(nèi)單元之間的影響，這樣減少了由于損傷局部化引起的單元網(wǎng)格敏感性，保證了數(shù)值模擬結(jié)果的客觀性.而在峰后階段(第40計算步后)載荷-位移曲線相差較大，其原因可能是本文提出的損傷模型不涉及峰后階段的巖石力學(xué)行為，該部分工作將于今后開展.

圖10 普通和正則化三點彎曲模型的載荷-位移曲線

4 結(jié) 論

1) 對于普通微平面損傷模型，材料越均勻，其損傷區(qū)域越集中，裂紋擴展模式越穩(wěn)定.

2) 普通微平面損傷模型和正則化模型的損傷模式基本一致，但正則化模型通過引入不同特征長度來表征損傷程度，對于給定的數(shù)值模型，存在一個最優(yōu)的特征長度來得到較好的損傷效果.一般而言，特征長度與模型網(wǎng)格長度相近時損傷效果較好.

3) 正則化方法一定程度上消除了局部損傷模型的網(wǎng)格尺寸效應(yīng).該方法通過考慮特征長度范圍內(nèi)的積分點的非局部效應(yīng)，對網(wǎng)格敏感性進行限制，保證了數(shù)值模擬結(jié)果的客觀性.