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

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

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

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

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

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

1.1 微平面模型

微平面模型通過應(yīng)力率和應(yīng)變率之間的關(guān)系構(gòu)建微平面上本構(gòu)關(guān)系來描述宏觀應(yīng)力應(yīng)變與時(shí)間之間的變化關(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個(gè)積分點(diǎn)微平面模型，其空間形狀為一正二十面體，積分點(diǎn)(微平面)的位置為其頂點(diǎn)(12個(gè))和棱中點(diǎn)(30個(gè))[21].以應(yīng)變率和應(yīng)力率形式表示的微平面本構(gòu)關(guān)系(式(3))能夠很好地描述模型應(yīng)力應(yīng)變與時(shí)間的聯(lián)系[22]，該關(guān)系在COMSOL中通過添加PDE物理場來實(shí)現(xiàn).

(3)

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

(4)

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

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

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

(5)

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

1.2 微平面損傷模型

基于應(yīng)變等效假設(shè)[23]的損傷在細(xì)觀損傷力學(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)基于屈服準(zhǔn)則判定.首先，對于處于拉應(yīng)力狀態(tài)時(shí)的節(jié)點(diǎn)，即法向應(yīng)力為正(σN> 0)的情況采用最大拉應(yīng)力準(zhǔn)則，認(rèn)為某個(gè)微平面上的最大法向應(yīng)力超過閾值時(shí)，該微平面系統(tǒng)對應(yīng)的節(jié)點(diǎn)即發(fā)生損傷，可表示為

(7)

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

(8)

D=D(n) .

(9)

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

(10)

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

(11)

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

1.3 正則化方法引入

(12)

式中,?n為n階梯度.又因?yàn)?/p>

(13)

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

(14)

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

(15)

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

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

(16)

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

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

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

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

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

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

2.2 模型驗(yàn)證

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

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

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

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

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

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

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

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

3 正則化模型分析

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

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

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

圖6 三點(diǎn)彎曲實(shí)驗(yàn)試件尺寸(單位：mm)

進(jìn)行對比分析.除加載速率外數(shù)值模型參數(shù)見表1.為提高模型計(jì)算速度并防止網(wǎng)格收斂性問題，模型網(wǎng)格劃分采用COMSOL默認(rèn)的四面體網(wǎng)格(Mesh1)，最大單元尺寸為3.15×10-3m，最小單元尺寸為1.35×10-4m，最大單元增長率為1.35，單元總數(shù)為70 109個(gè).圖7給出了不同非均質(zhì)系數(shù)m的試件損傷演化.由圖7可見，非均質(zhì)系數(shù)對試件的損傷破壞影響明顯，非均質(zhì)系數(shù)低(m=1)時(shí)，巖樣中損傷在預(yù)置裂縫周圍彌散分布；非均質(zhì)系數(shù)高(5，10)時(shí)，巖樣中損傷集中于預(yù)置裂縫尖端，形成明顯的局部化損傷區(qū)域.非均質(zhì)系數(shù)為1時(shí)，損傷區(qū)域雖然在預(yù)制裂縫處較多，但是在試件其他區(qū)域也有大量損傷區(qū)域存在，與實(shí)際情況不符.對于巖石材料而言，非均質(zhì)系數(shù)為5左右時(shí)結(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進(jìn)行數(shù)值模擬，數(shù)值模型非均質(zhì)系數(shù)均為5，網(wǎng)格劃分方式(Mesh1)是由COMSOL根據(jù)物理場自動生成的四面體形式.圖8為具有不同特征長度的數(shù)值模型的損傷云圖(時(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 時(shí)損傷位置和損傷發(fā)展方向準(zhǔn)確，但最終損傷效果與特征長度為0.5 mm的結(jié)果相近，因此模型單元尺寸與特征長度相近時(shí)，可認(rèn)為對損傷結(jié)果影響較小.特征長度(l=2 mm)過大時(shí)，由于特征長度大于大部分網(wǎng)格尺寸，導(dǎo)致?lián)p傷面積過大而不符合實(shí)際情況.由此可見，在相同網(wǎng)格劃分條件下，存在一個(gè)最優(yōu)特征長度值，可以得到最好的裂紋擴(kuò)展效果.在可預(yù)見的范圍內(nèi)，對相同尺寸的數(shù)值模型而言，隨著網(wǎng)格密度增大，最小單元尺寸減小，最優(yōu)特征長度值也會隨之減小.

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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