王慶，姚俊，譚文祿，潘惠惠

(中電建生態(tài)環(huán)境集團(tuán)有限公司，廣東深圳518102)

巖石等復(fù)合材料是在小到一定的尺度上都具有非均質(zhì)性[1]，這些微細(xì)觀尺度上的非均質(zhì)性對(duì)材料的宏觀物理力學(xué)性質(zhì)具有很大的影響。由于微細(xì)觀結(jié)構(gòu)的復(fù)雜性和多樣性，采用實(shí)驗(yàn)的方法測(cè)量不同微細(xì)觀結(jié)構(gòu)的宏觀等效性質(zhì)需要花費(fèi)大量的時(shí)間和代價(jià)，因此許多學(xué)者傾向于尋求多尺度計(jì)算方法[2,3]來(lái)建立復(fù)合材料的微細(xì)觀結(jié)構(gòu)和宏觀等效性質(zhì)之間的關(guān)系。

多尺度計(jì)算方法有許多不同的分類方法[4-6]，根據(jù)宏細(xì)觀模型間信息的傳遞方式不同，主要分為3類：遞階法[7](Hierarchical multiscale method)，準(zhǔn)協(xié)同法[8](Semi-concurrent multiscale method)和協(xié)同法[9](Concurrent multiscale method)。遞階法在細(xì)觀模型數(shù)值試驗(yàn)之前，首先假定一種典型的宏觀本構(gòu)模型，如彈性、塑性、黏塑性等，這些宏觀本構(gòu)模型中的未知參數(shù)，可以根據(jù)在細(xì)觀模型上進(jìn)行的數(shù)值試驗(yàn)結(jié)果數(shù)據(jù)擬合得到，遞階法只能處理線性問(wèn)題和宏觀本構(gòu)關(guān)系簡(jiǎn)單的非線性問(wèn)題。準(zhǔn)協(xié)同方法的宏觀模型和細(xì)觀模型是分開(kāi)計(jì)算的，在計(jì)算過(guò)程中，宏細(xì)觀模型之間信息是雙向傳遞的，不需要事先知道宏觀的本構(gòu)關(guān)系，因此可以處理復(fù)雜的非線性問(wèn)題。協(xié)同方法的特點(diǎn)是能將細(xì)觀模型的網(wǎng)格直接嵌套在宏觀模型的網(wǎng)格中，細(xì)觀模型的求解直接在宏觀模型中進(jìn)行，所以不需要進(jìn)行尺度分離，但需要大量的存儲(chǔ)空間和計(jì)算代價(jià)。因此選擇遞階法既可以處理?yè)p傷、裂紋擴(kuò)展等復(fù)雜非線性問(wèn)題，也可以避免大量的存儲(chǔ)空間和計(jì)算代價(jià)。

最常見(jiàn)遞階法是計(jì)算均勻化方法，包括線性和非線性，其中線性計(jì)算均勻化方法只能處理線性問(wèn)題，但是巖石等材料在變形過(guò)程中存在細(xì)觀結(jié)構(gòu)的非均質(zhì)性以及不同組分之間的相互作用，這種不同組分之間的相互作用在斷裂多尺度計(jì)算中必須要考慮，因此本文選擇用非線性計(jì)算均勻化方法來(lái)研究巖石等材料的裂紋擴(kuò)展規(guī)律。

1 非線性多尺度控制方程

1.1 增量控制方程

非線性計(jì)算均勻化方法是基于漸進(jìn)展開(kāi)[10]的思想得到的，其與線性計(jì)算均勻化方法控制方程的推導(dǎo)過(guò)程[11]最大的區(qū)別是細(xì)觀結(jié)構(gòu)在變形過(guò)程中材料參數(shù)會(huì)發(fā)生變化，所以控制方程必須是增量形式?？紤]一種非均質(zhì)的材料的宏觀空間區(qū)域是Ωε，邊界為Γε。該材料具有某種周期性的細(xì)觀結(jié)構(gòu)，該細(xì)觀結(jié)構(gòu)具有表征體積單元RVE，其細(xì)觀空間區(qū)域表示為Θ。Ωε中該材料的基本控制方程可以表示為：

(1)

(2)

(3)

位移和力邊界條件為，

(4)

(5)

1.2 增量位移的漸進(jìn)展開(kāi)

按照漸進(jìn)展開(kāi)理論，宏細(xì)觀模型坐標(biāo)是相關(guān)的[12]，假設(shè)宏觀模型坐標(biāo)為x，細(xì)觀模型坐標(biāo)為y，增量形式的位移、應(yīng)力和應(yīng)變可以用尺度因子ε形式表示為：

(6)

(7)

(8)

(9)

(10)

(11)

將增量應(yīng)變式(7)代入到本構(gòu)關(guān)系式(2)中可以得到增量應(yīng)力：

(12)

將增量應(yīng)力式(8)代入到平衡方程式(1)中，根據(jù)尺度因子ε進(jìn)行組合，可以得到不同階的平衡方程：

(13)

(14)

(15)

1.3 宏細(xì)觀控制方程

求解不同階的平衡方程式(13)—(15)，得到宏細(xì)觀控制方程和信息傳遞方程如下。

a) 宏觀問(wèn)題

(16)

b) 細(xì)觀問(wèn)題

(17)

c) 宏細(xì)觀信息傳遞

(18)

(19)

式中N + 1〈σij〉y、N〈σij〉y——第N、N+1步的宏觀應(yīng)力。

2 非線性計(jì)算均勻化方法的實(shí)現(xiàn)

2.1 一般過(guò)程

非線性計(jì)算均勻化方法的一般過(guò)程可以分為4步，見(jiàn)圖1。利用細(xì)觀模型計(jì)算出來(lái)的等效切線模量和等效應(yīng)力求解宏觀方程(16)，得到宏觀應(yīng)變；將第1步中得到的宏觀應(yīng)變作為邊界條件施加到細(xì)觀模型上，求解細(xì)觀方程(17)，得到細(xì)觀模型的應(yīng)力等狀態(tài)量；根據(jù)第2步的計(jì)算結(jié)果按式(18)、(19)求解宏觀模型的等效切線模量和等效應(yīng)力；回到第1步進(jìn)行下一個(gè)增量步的循環(huán)，直到加載完成。

2.2 周期性RVE的生成

本文的數(shù)值案例中使用的材料是水泥砂漿，具體成分是砂子、水泥和水，其水灰比為0.4，砂子重量占40%，水泥砂漿和砂子的密度分別為2 000、2 650 kg/m3，在不考慮氣孔的情況下，細(xì)觀結(jié)構(gòu)中基質(zhì)和砂子的體積比分別為67%和33%。砂子采用的粒徑小于0.5 mm的河砂，假設(shè)砂子的粒徑滿足在區(qū)間[0.3, 0.5] mm之間的均勻分布。

從非線性多尺度控制的邊界條件中可以看出，RVE的細(xì)觀結(jié)構(gòu)必須具有周期性。選擇用第順序投放方法[13]來(lái)生成水泥砂漿的細(xì)觀模型，但是為了使產(chǎn)生的細(xì)觀模型具有周期性，所以要在順序投放方法基礎(chǔ)上進(jìn)行改進(jìn)。周期性RVE中的球體主要分為4類，見(jiàn)圖2。

a) 內(nèi)部球體。遠(yuǎn)離RVE邊界的球體是內(nèi)部球體，可以用傳統(tǒng)的順序投放方法產(chǎn)生。

b) 面上的球體。當(dāng)1個(gè)球體與RVE的面接觸時(shí)，會(huì)在對(duì)應(yīng)的面上復(fù)制1個(gè)相同的球體，原來(lái)的球體稱為主體，復(fù)制的球體稱為從體。

c) 棱上的球體。當(dāng)1個(gè)球體與RVE的棱接觸時(shí)，會(huì)在其它3條棱上復(fù)制3個(gè)從體。

d) 頂點(diǎn)上的球體。當(dāng)1個(gè)球體與RVE的頂點(diǎn)接觸時(shí)，會(huì)在其它7個(gè)頂點(diǎn)上復(fù)制7個(gè)從體。

以大小為(4×4×4) mm3的RVE舉例說(shuō)明周期性RVE的生成，水泥砂漿的細(xì)觀結(jié)構(gòu)見(jiàn)圖3。

2.3 一般過(guò)程的實(shí)現(xiàn)

在非線性計(jì)算均勻化方法中，宏觀模型的本構(gòu)關(guān)系是不需要知道的，可以在加載過(guò)程中通過(guò)細(xì)觀模型確定的。宏觀應(yīng)力和切線模量是隨著RVE的變形而變化的，在不知道顯式本構(gòu)關(guān)系的情況下求解宏觀問(wèn)題，必須依靠ABAQUS中的用戶子程序UMAT。UMAT有2個(gè)變量需要用戶返回給ABAQUS主程序，一個(gè)是Jacobian矩陣?Δσ/?Δe，這與宏細(xì)觀傳遞信息中的切線模量相對(duì)應(yīng)；另一個(gè)是增量應(yīng)力，這與宏細(xì)觀傳遞信息中的增量應(yīng)力相對(duì)應(yīng)。UMAT的調(diào)用都是在單元的高斯積分點(diǎn)上進(jìn)行的，在第k+1步的宏觀模型中的每個(gè)單元的高斯積分點(diǎn)中，ABAQUS調(diào)用UMAT的流程見(jiàn)圖4。

3 非充填預(yù)裂紋的單軸壓縮試驗(yàn)

3.1 試樣及加載系統(tǒng)

預(yù)裂紋壓縮試驗(yàn)中用到的試件材料是水泥砂漿，其配比在前面有所介紹，其它試驗(yàn)條件和結(jié)果可見(jiàn)Zhuang的預(yù)裂紋試驗(yàn)[14]。試樣均采用70 mm×70 mm×140 mm尺寸的標(biāo)準(zhǔn)長(zhǎng)方體，預(yù)裂紋為貫通型，試樣及預(yù)裂紋的幾何形態(tài)見(jiàn)圖5。裂紋厚度為1 mm，長(zhǎng)度為15 mm，傾角α在有15、30、45、60、75°幾種，主要研究不同傾角下的裂紋擴(kuò)展規(guī)律。試驗(yàn)采用無(wú)側(cè)限單軸加載，加載速率為0.25 kN/s，預(yù)加荷載為1 kN，加載5 kN為一步，每步完成后停止加載5 s，進(jìn)行裂隙周邊拍照記錄。

3.2 非充填預(yù)裂紋的擴(kuò)展特征

按裂紋擴(kuò)展的幾何形態(tài)和破壞機(jī)制的不同，將原始裂隙周邊新生裂紋分為3類[15]：翼裂紋、反翼裂紋與次生裂紋，見(jiàn)圖6。新生裂紋的擴(kuò)展特征主要包括起裂點(diǎn)，起裂順序和起裂角，其中翼裂紋的起裂角定義見(jiàn)圖7。

4 數(shù)值算例

4.1 多尺度模型

非充填預(yù)裂紋的單軸壓縮試驗(yàn)的多尺度模型見(jiàn)圖8。由于非線性多尺度計(jì)算是在每個(gè)單元的高斯積分點(diǎn)上進(jìn)行，需要的計(jì)算代價(jià)比較大，因此只考慮預(yù)裂紋附近的材料的非均勻性，用紅色單元表示；遠(yuǎn)離預(yù)裂紋(距離超過(guò)3倍預(yù)裂紋厚度)的材料則是均勻材料，用綠色單元表示。多尺度材料的細(xì)觀結(jié)構(gòu)見(jiàn)圖3。宏觀模型的下邊界約束y方向位移，左右邊界自由，上邊界施加位移加載方式。

4.2 本構(gòu)關(guān)系和參數(shù)

細(xì)觀模型中的基質(zhì)和顆粒采用混凝土塑性損傷本構(gòu)模型[16]，其單軸拉伸和壓縮應(yīng)力應(yīng)變曲線見(jiàn)圖9。

出現(xiàn)損傷之后，材料的彈性模量會(huì)出現(xiàn)退化，彈性模量的退化程度取決于等效塑性應(yīng)變，拉壓應(yīng)力應(yīng)變關(guān)系如下：

(20)

為了得到損傷演化規(guī)律，將式(20)變形為：

(21)

拉壓應(yīng)變?chǔ)與和εt可以分解為彈性和非彈性應(yīng)變的和：

(22)

(23)

d=1-(1-dt)(1-dc)

(24)

混凝土塑性損傷本構(gòu)中參數(shù)主要有，彈性模量E0，泊松比v，抗壓強(qiáng)度σcu，抗拉強(qiáng)度σtu，常數(shù)bt和bc。這些參數(shù)可以通過(guò)水泥砂漿室內(nèi)單軸壓縮試驗(yàn)與非計(jì)算均勻化方法的對(duì)比進(jìn)行調(diào)試，細(xì)觀組分的材料參數(shù)見(jiàn)表1，bt和bc分別取0.1和0.7。

表1 水泥砂漿細(xì)觀組分的材料參數(shù)

室內(nèi)試驗(yàn)和數(shù)值模擬的單軸壓縮的應(yīng)力-應(yīng)變曲線見(jiàn)圖10。從圖中可以看出，數(shù)值模擬結(jié)果與室內(nèi)試驗(yàn)結(jié)果整體比較吻合，但是由于實(shí)際中水泥砂漿中包含了初始裂紋和損傷，在彈性初始階段會(huì)出現(xiàn)裂紋壓密階段，這個(gè)在塑性損傷本構(gòu)關(guān)系中無(wú)法再現(xiàn)。另外由于是剛性試驗(yàn)機(jī)，材料達(dá)到抗壓強(qiáng)度之后會(huì)突然破壞，導(dǎo)致下降段很陡。

圖8中的非均勻材料用多尺度方法進(jìn)行計(jì)算，其細(xì)觀模型的材料參數(shù)見(jiàn)表1；均勻材料采用單一尺度模型進(jìn)行計(jì)算，其本構(gòu)關(guān)系也采用混凝土塑性損傷模型，材料參數(shù)取圖10中的試驗(yàn)數(shù)據(jù)(虛線)。

4.3 結(jié)果分析

由于非線性計(jì)算均勻化方法不能處理非連續(xù)問(wèn)題，將用漸進(jìn)損傷過(guò)程來(lái)代替裂紋擴(kuò)展路徑的研究。由于計(jì)算均勻化方法的宏觀本構(gòu)模型是未知的，需要一個(gè)代表?yè)p傷的變量來(lái)研究漸進(jìn)損傷過(guò)程。幾何損傷理論認(rèn)為[17]，損傷模型應(yīng)力的折減等于未損傷模型受力面積的折減。根據(jù)均勻化準(zhǔn)則，宏觀損傷因子可以表示為：

(25)

預(yù)裂紋傾角為45°的多尺度模型的損傷因子見(jiàn)圖11。傾角為45°的試件的最終破壞模式見(jiàn)圖12，可以看出數(shù)值模擬和室內(nèi)試驗(yàn)的裂紋擴(kuò)展路徑相似。在數(shù)值模型中有2條反翼裂紋，而在室內(nèi)試驗(yàn)中只有1條，主要原因是室內(nèi)試驗(yàn)中的試件不是均勻的。不同地方的強(qiáng)度可能不一樣，對(duì)應(yīng)的細(xì)觀結(jié)構(gòu)也不一樣，本文的非線性計(jì)算均勻化方法可以考慮不同的細(xì)觀結(jié)構(gòu)，但由于缺少細(xì)觀結(jié)構(gòu)的信息，所以只能暫時(shí)假定所有RVE結(jié)構(gòu)都相同。雖然如此，數(shù)值計(jì)算結(jié)果與室內(nèi)試驗(yàn)基本一致，說(shuō)明非線性計(jì)算均勻化方法處理預(yù)裂紋擴(kuò)展的適用性。模型中非均勻材料處的有限元單元用本文的多尺度模型進(jìn)行計(jì)算時(shí)，每個(gè)增量步中每個(gè)單元高斯積分點(diǎn)的計(jì)算時(shí)間需要30 min左右，如果單元和增量步過(guò)多，模型加載完成需要數(shù)周；另外，如果某個(gè)單元對(duì)應(yīng)的細(xì)觀模型在某個(gè)增量步不收斂，將會(huì)導(dǎo)致整個(gè)多尺度模型的計(jì)算終止，因此本文的計(jì)算方法計(jì)算代價(jià)較大，需要在多尺度單元數(shù)量和計(jì)算代價(jià)間進(jìn)行平衡，而且受細(xì)觀模型的收斂性影響，需要注意選擇更加容易收斂的細(xì)觀本構(gòu)模型。

不同原始裂隙傾角下的翼裂紋起裂角見(jiàn)表2。與室內(nèi)試驗(yàn)結(jié)果相同，隨著預(yù)裂紋傾角α的增大，翼裂紋起裂角逐漸減小，見(jiàn)圖13，即翼裂紋向加載方向的偏轉(zhuǎn)程度逐漸增大。

表2 不同原始裂隙傾角下的翼裂紋起裂角對(duì)比 單位：(°)

5 結(jié)論

本文提出了一種非線性計(jì)算均勻化方法，將非線性計(jì)算均勻化的一般步驟在ABAQUS中實(shí)現(xiàn)，并將其運(yùn)用到非填充預(yù)裂紋擴(kuò)展的多尺度模擬，得出以下結(jié)論。

a) 數(shù)值模擬和室內(nèi)試驗(yàn)的單軸壓縮應(yīng)力—應(yīng)變曲線擬合的很好，說(shuō)明混凝土塑性損傷模型可以作為水泥砂漿等脆性材料的細(xì)觀組分的本構(gòu)模型。

b) 在預(yù)裂紋擴(kuò)展數(shù)值算例中，不同傾角的裂紋的擴(kuò)展路徑與試驗(yàn)結(jié)果比較吻合，表明非線性計(jì)算均勻化方法在處理多尺度非線性問(wèn)題是一個(gè)有效的方法，且用損傷來(lái)模擬裂紋擴(kuò)展是可行的。

c) 隨著原始裂隙傾角的增大，翼裂紋起裂角逐漸減小，即翼裂紋向加載方向的偏轉(zhuǎn)程度逐漸增大。