基于非線性均勻化的類巖石材料多尺度計算方法研究

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

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

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

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

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

1 非線性多尺度控制方程

1.1 增量控制方程

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

(1)

(2)

(3)

位移和力邊界條件為，

(4)

(5)

1.2 增量位移的漸進展開

按照漸進展開理論，宏細觀模型坐標是相關的[12]，假設宏觀模型坐標為x，細觀模型坐標為y，增量形式的位移、應力和應變可以用尺度因子ε形式表示為：

(6)

(7)

(8)

(9)

(10)

(11)

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

(12)

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

(13)

(14)

(15)

1.3 宏細觀控制方程

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

a) 宏觀問題

(16)

b) 細觀問題

(17)

c) 宏細觀信息傳遞

(18)

(19)

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

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

2.1 一般過程

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

2.2 周期性RVE的生成

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

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

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

b) 面上的球體。當1個球體與RVE的面接觸時，會在對應的面上復制1個相同的球體，原來的球體稱為主體，復制的球體稱為從體。

c) 棱上的球體。當1個球體與RVE的棱接觸時，會在其它3條棱上復制3個從體。

d) 頂點上的球體。當1個球體與RVE的頂點接觸時，會在其它7個頂點上復制7個從體。

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

2.3 一般過程的實現(xiàn)

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

3 非充填預裂紋的單軸壓縮試驗

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

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

3.2 非充填預裂紋的擴展特征

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

4 數(shù)值算例

4.1 多尺度模型

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

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

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

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

(20)

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

(21)

拉壓應變εc和εt可以分解為彈性和非彈性應變的和：

(22)

(23)

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

(24)

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

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

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

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

4.3 結(jié)果分析

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

(25)

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

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

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

5 結(jié)論

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

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

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

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