徐 磊，崔姍姍，姜 磊，任青文

(河海大學(xué)水利水電學(xué)院，江蘇，南京 210098)

混凝土是典型的隨機(jī)多尺度準(zhǔn)脆性材料[1]，在細(xì)觀尺度上，通常被視為由(粗)骨料、砂漿及兩者之間的界面過渡區(qū)(Interfacial Transition Zone,ITZ)構(gòu)成的三相非均勻復(fù)合材料[2 ? 4]。雖然在彈性階段，可將混凝土作為均勻材料并采用宏觀本構(gòu)模型描述其受力變形行為[5]，但在損傷開裂階段，基于“均勻性”假定的宏觀本構(gòu)模型難以準(zhǔn)確描述混凝土復(fù)雜的非線性力學(xué)行為[6]，主要原因是混凝土的損傷開裂演化與其細(xì)觀材料結(jié)構(gòu)直接相關(guān)，表現(xiàn)出明顯的隨機(jī)、非均勻、局部化與跨尺度特征。因此，準(zhǔn)確分析混凝土從細(xì)觀裂紋萌生、擴(kuò)展、集聚至宏觀裂縫形成這一復(fù)雜過程需要考慮其細(xì)觀材料結(jié)構(gòu)[7 ? 8]。

理論上，雖然直接建立混凝土結(jié)構(gòu)的細(xì)觀計(jì)算模型可以充分體現(xiàn)細(xì)觀材料結(jié)構(gòu)對(duì)宏觀結(jié)構(gòu)行為的影響，但卻由于該方法對(duì)計(jì)算資源的極高要求而難以實(shí)施[9]。而在實(shí)際混凝土結(jié)構(gòu)中，一般僅有范圍較小的局部區(qū)域(損傷區(qū))會(huì)進(jìn)入非線性階段，其他大部分區(qū)域(彈性區(qū))處于彈性階段[10]，如圖1所示。因此，為兼顧分析精度與效率，一種可行方法是在細(xì)觀尺度下建立損傷區(qū)計(jì)算模型，在宏觀尺度下建立彈性區(qū)計(jì)算模型，并通過尺度連接將上述不同尺度下的局部計(jì)算模型連接起來(lái)以形成整體宏細(xì)觀協(xié)同計(jì)算模型(見圖1)，從而用相對(duì)較少的計(jì)算資源實(shí)現(xiàn)混凝土結(jié)構(gòu)跨尺度損傷開裂演化過程的準(zhǔn)確分析。

圖1 結(jié)構(gòu)分區(qū)與宏細(xì)觀協(xié)同計(jì)算模型示意圖Fig.1 Sketches of structure decomposition and macro-meso-scale concurrent computational model

Eckardt和K?nke[11]采用約束方程法實(shí)現(xiàn)宏細(xì)觀尺度連接，在有限單元法框架內(nèi)提出了混凝土損傷分析的非均勻多尺度方法；Unger和Eckardt[12]對(duì)比分析了約束方程法、Mortar法及Arlequin法等尺度連接方法的優(yōu)缺點(diǎn)，建立了混凝土自適應(yīng)宏細(xì)觀協(xié)同多尺度計(jì)算模型，但所采用的以最大拉應(yīng)力為指標(biāo)的分析尺度轉(zhuǎn)換準(zhǔn)則不適用于復(fù)雜應(yīng)力狀態(tài)；Lloberas-Valls和Rixen等[13]在區(qū)域分解法框架內(nèi)，對(duì)比分析了區(qū)域間非重疊網(wǎng)格的強(qiáng)、弱尺度連接方法，并提出了一種改進(jìn)的弱尺度連接方法；Sun和Li[14]在通過采用均勻宏觀網(wǎng)格簡(jiǎn)化宏-細(xì)觀界面動(dòng)態(tài)調(diào)整的基礎(chǔ)上模擬了混凝土柱在動(dòng)力荷載作用下的自適應(yīng)跨尺度破壞過程。為簡(jiǎn)化細(xì)觀建模、便于形成非重疊網(wǎng)格與實(shí)施宏-細(xì)觀尺度連接，以上方法均采用了均勻規(guī)則的宏觀網(wǎng)格。Rodrigues和Manzoli等[15]實(shí)現(xiàn)了基于非協(xié)調(diào)重疊網(wǎng)格的宏-細(xì)觀尺度連接，但需在協(xié)同計(jì)算模型中引入專門用于施加位移約束的耦合單元，增加了數(shù)值實(shí)施的難度。

本文提出了一種基于雙重網(wǎng)格的混凝土自適應(yīng)宏細(xì)觀協(xié)同有限元分析方法，基本思路是通過布置獨(dú)立剖分的兩套有限元網(wǎng)格，在分析域內(nèi)分別形成將混凝土視為均勻材料的宏觀(尺度)模型和非均勻材料的細(xì)觀(尺度)模型；通過提出基于Ottosen多軸強(qiáng)度準(zhǔn)則[16]的分析尺度自適應(yīng)轉(zhuǎn)換準(zhǔn)則，在分析過程中動(dòng)態(tài)更新宏細(xì)觀協(xié)同有限元模型；通過提出基于形函數(shù)插值的多點(diǎn)位移約束方法，實(shí)現(xiàn)宏細(xì)觀非協(xié)調(diào)重疊網(wǎng)格連接；在此基礎(chǔ)上，給出了基于雙重網(wǎng)格的混凝土自適應(yīng)宏細(xì)觀協(xié)同有限元求解流程，并在MATLAB平臺(tái)上完成了程序開發(fā)。算例分析表明，采用本文所提出的方法可在兼顧效率與精度的前提下，實(shí)現(xiàn)考慮細(xì)觀材料結(jié)構(gòu)的混凝土損傷開裂跨尺度演化過程自適應(yīng)分析。

1 混凝土細(xì)觀有限元模型

為了應(yīng)用有限單元法開展細(xì)觀尺度下的混凝土損傷開裂分析，首先需要建立混凝土細(xì)觀有限元模型，主要涉及細(xì)觀結(jié)構(gòu)模擬、有限元網(wǎng)格剖分和細(xì)觀力學(xué)模型等3個(gè)方面，分述如下。

1.1 細(xì)觀結(jié)構(gòu)模擬

混凝土在細(xì)觀尺度上的材料結(jié)構(gòu)主要取決于(粗)骨料的形狀及其含量、粒徑、級(jí)配等控制參數(shù)?；谙壬呻S機(jī)骨料后進(jìn)行骨料投放的隨機(jī)取放法[17]，研發(fā)了混凝土細(xì)觀結(jié)構(gòu)隨機(jī)生成軟件AutoGMC。對(duì)于任意形狀的模擬區(qū)域，該軟件可依據(jù)給定的控制參數(shù)在模擬區(qū)域內(nèi)完成圓形或多邊形骨料細(xì)觀結(jié)構(gòu)的隨機(jī)生成，其中，圓形骨料可用于近似模擬天然(卵石)骨料，多邊形骨料可用于模擬人工(碎石)骨料。圖2給出了不同試件形狀和結(jié)構(gòu)型式下的細(xì)觀結(jié)構(gòu)生成實(shí)例。

圖2 細(xì)觀結(jié)構(gòu)生成實(shí)例Fig.2 Examples of meso-structure generation

1.2 有限元網(wǎng)格剖分

為建立細(xì)觀有限元模型，需對(duì)混凝土細(xì)觀結(jié)構(gòu)進(jìn)行網(wǎng)格剖分，本文利用ABAQUS前處理模塊，通過MATLAB和PYTHON混合編程開發(fā)了混凝土細(xì)觀有限元網(wǎng)格自動(dòng)剖分程序。具體而言，首先基于模擬區(qū)域和骨料的幾何信息，依據(jù)ABAQUS規(guī)定的編寫規(guī)則[18]，由程序自動(dòng)編寫可被ABAQUS前處理模塊執(zhí)行的PYTHON腳本；進(jìn)一步地，通過MATLAB調(diào)用ABAQUS前處理模塊生成僅包含骨料與砂漿單元的兩相網(wǎng)格；在此基礎(chǔ)上，為模擬骨料與砂漿之間的ITZ，收縮兩相網(wǎng)格中的骨料邊界，并在骨料單元與砂漿單元之間嵌入具有一定厚度(取為100 μm)[19]的ITZ單元，從而形成最終的三相網(wǎng)格，如圖3所示。由于在砂漿單元與骨料單元間插入的ITZ單元所占據(jù)的空間是原兩相網(wǎng)格有限元模型中骨料單元的一部分，故為保證三相網(wǎng)格有限元模型中的骨料粒徑與要求的一致，應(yīng)在細(xì)觀結(jié)構(gòu)生成過程中，將界面過渡區(qū)作為骨料的一部分，即通過增大骨料粒徑的方式使得骨料在細(xì)觀結(jié)構(gòu)中占據(jù)的空間既包括骨料自身，又包括骨料周圍的界面過渡區(qū)。采用上述程序完成了圖2(c)所示細(xì)觀結(jié)構(gòu)的有限元網(wǎng)格剖分，見圖4。

圖3 界面過渡區(qū)單元生成Fig.3 Generation of ITZ element

圖4 細(xì)觀有限元網(wǎng)格剖分實(shí)例Fig.4 Example of mesoscale finite element meshing

1.3 細(xì)觀力學(xué)模型

對(duì)于混凝土細(xì)觀各相材料，還需確定其適用的本構(gòu)模型。由于普通混凝土損傷開裂通常是在ITZ中萌生并向砂漿中擴(kuò)展，而(硬)骨料一般不會(huì)發(fā)生破壞[20]。因此，可將骨料視為線彈性材料，但需考慮砂漿與界面過渡區(qū)的非線性力學(xué)行為[21 ? 22]。本文采用塑性損傷模型(CDP模型)[23]作為砂漿與ITZ的本構(gòu)模型。CDP模型應(yīng)力-應(yīng)變關(guān)系表達(dá)式如下：

式中：ω為偏心率，用于描述塑性勢(shì)函數(shù)向其漸近線逼近的速度，一般可取為0.1； σt0為單軸抗拉強(qiáng)度；ψ為膨脹角；J2為有效應(yīng)力張量偏量的第二不變量；I1為有效應(yīng)力張量的第一不變量。

CDP模型采用如下形式的屈服函數(shù)：

引入拉伸、壓縮損傷因子dt、dc分別表征拉伸、壓縮損傷導(dǎo)致的剛度退化，其量值分別隨拉伸、壓縮等效塑性應(yīng)變的變化而變化。進(jìn)一步考慮應(yīng)力反向后的剛度恢復(fù)效應(yīng)，即可給出復(fù)雜應(yīng)力狀態(tài)下d與dt、dc之間的關(guān)系式：

式中：st、sc的取值與應(yīng)力狀態(tài)相關(guān)[24]。單軸受拉時(shí)，st=0 ，sc=1 ， 故d=dt；單軸受壓時(shí)，sc=0，st=1 ， 故d=dc。

2 宏細(xì)觀協(xié)同有限元分析方法

為了在兼顧效率與精度的前提下準(zhǔn)確分析混凝土損傷開裂的跨尺度演化過程，提出一種基于雙重網(wǎng)格的混凝土自適應(yīng)宏細(xì)觀協(xié)同有限元分析方法，詳述如下。

2.1 基于雙重網(wǎng)格的宏細(xì)觀協(xié)同有限元模型

如圖5(a)所示，為建立宏細(xì)觀協(xié)同有限元模型，在分析域內(nèi)布置兩套有限元網(wǎng)格，分別為宏觀(尺度)網(wǎng)格和細(xì)觀(尺度)網(wǎng)格，故稱雙重網(wǎng)格。宏觀網(wǎng)格和細(xì)觀網(wǎng)格獨(dú)立剖分，在剖分宏觀網(wǎng)格時(shí)，混凝土被視為均勻線彈性材料，而在剖分細(xì)觀網(wǎng)格時(shí)，則將混凝土視為由(粗)骨料、砂漿和界面過渡區(qū)組成的非均勻材料。在此基礎(chǔ)上，將線彈性本構(gòu)模型及參數(shù)賦予宏觀網(wǎng)格中的各單元，即可形成分析域的宏觀有限元模型；類似地，將細(xì)觀各組分的本構(gòu)模型和參數(shù)賦予細(xì)觀網(wǎng)格中的相應(yīng)單元，即可形成分析域的細(xì)觀有限元模型。

如圖5(b)所示，在宏細(xì)觀協(xié)同分析中，僅有部分宏觀模型被作為整體模型的一部分，其余部分則被替換為與之相應(yīng)的細(xì)觀模型，從而形成宏細(xì)觀協(xié)同分析整體有限元模型；上述協(xié)同模型需依據(jù)分析對(duì)象受力狀態(tài)變化動(dòng)態(tài)更新，具體而言，當(dāng)某宏觀單元內(nèi)任一積分點(diǎn)的應(yīng)力滿足分析尺度自適應(yīng)轉(zhuǎn)換準(zhǔn)則(詳見節(jié)2.2)時(shí)，即需將該宏觀單元從協(xié)同模型中消除并激活與之相應(yīng)的細(xì)觀單元集合。

圖5 宏細(xì)觀協(xié)同有限元模型Fig.5 Macro-meso-scale concurrent finite element model

由于宏觀網(wǎng)格和細(xì)觀網(wǎng)格的剖分密度差異通常很大且剖分過程相互獨(dú)立，故在宏細(xì)觀協(xié)同有限元模型中，宏觀模型與細(xì)觀模型連接處的有限元網(wǎng)格不但是非協(xié)調(diào)的，而且會(huì)出現(xiàn)一定程度的重疊現(xiàn)象。因此，為實(shí)現(xiàn)協(xié)同有限元分析，需通過非協(xié)調(diào)重疊網(wǎng)格連接(詳見2.3節(jié))來(lái)保證宏觀模型與細(xì)觀模型連接處的變形協(xié)調(diào)[15]。

2.2 自適應(yīng)尺度轉(zhuǎn)換

混凝土結(jié)構(gòu)的不均勻應(yīng)力分布與混凝土材料的應(yīng)變軟化特性決定了在混凝土結(jié)構(gòu)宏細(xì)觀協(xié)同有限元分析中，僅需對(duì)部分區(qū)域(損傷區(qū))開展細(xì)觀尺度分析[30]。但由于實(shí)際混凝土結(jié)構(gòu)受力狀態(tài)的復(fù)雜性，通常無(wú)法在分析前準(zhǔn)確確定損傷區(qū)的位置與范圍[31 ? 32]，故需在分析過程中依據(jù)結(jié)構(gòu)當(dāng)前受力狀態(tài)確定需要將分析尺度從宏觀轉(zhuǎn)換為細(xì)觀的區(qū)域并動(dòng)態(tài)更新宏細(xì)觀協(xié)同有限元模型，這一過程即為分析尺度的自適應(yīng)轉(zhuǎn)換。

為在分析過程中實(shí)現(xiàn)分析尺度的自適應(yīng)轉(zhuǎn)換，本文基于Ottosen多軸強(qiáng)度準(zhǔn)則[16]，提出了以積分點(diǎn)應(yīng)力為指標(biāo)的混凝土自適應(yīng)宏細(xì)觀尺度轉(zhuǎn)換準(zhǔn)則，如下式所示：

式中：θ為應(yīng)力Lode角；K為形狀因子，可按下式計(jì)算：

基于上述自適應(yīng)宏細(xì)觀尺度轉(zhuǎn)換準(zhǔn)則，即可在某一增量步迭代收斂后，依據(jù)各宏觀單元的當(dāng)前應(yīng)力狀態(tài)判斷是否存在需要進(jìn)行分析尺寸轉(zhuǎn)換的宏觀單元，若存在，則表明當(dāng)前宏細(xì)觀協(xié)同有限元模型的宏細(xì)觀區(qū)域劃分與應(yīng)力計(jì)算結(jié)果不符，需要更新宏細(xì)觀協(xié)同有限元模型并重新進(jìn)行該增量步的迭代求解；反之，若不存在要進(jìn)行分析尺寸轉(zhuǎn)換的宏觀單元，則表明當(dāng)前模型的宏細(xì)觀區(qū)域劃分與應(yīng)力計(jì)算結(jié)果相符，可進(jìn)行下一個(gè)增量步的迭代求解。

2.3 非協(xié)調(diào)重疊網(wǎng)格連接

在基于雙重網(wǎng)格的混凝土宏細(xì)觀協(xié)同有限元模型中，細(xì)觀模型網(wǎng)格的外圍結(jié)點(diǎn)位于宏觀單元內(nèi)部，致使宏細(xì)觀模型的有限元網(wǎng)格間存在重疊現(xiàn)象。為保證宏觀模型與細(xì)觀模型之間的變形協(xié)調(diào)，本文提出基于形函數(shù)插值的多點(diǎn)位移約束法來(lái)實(shí)現(xiàn)宏細(xì)觀非協(xié)調(diào)重疊網(wǎng)格之間的連接。為簡(jiǎn)明計(jì)，假定宏觀單元為三結(jié)點(diǎn)三角形單元，闡明該方法的基本思想。

如圖6所示，細(xì)觀模型某外圍結(jié)點(diǎn)P位于宏觀模型與細(xì)觀模型連接處的某宏觀單元e內(nèi)部，其位置坐標(biāo)為(xp,yp)。宏觀單元e各結(jié)點(diǎn)在平面直角坐標(biāo)系(x,y)中的x、y向位移分別為ui、vi，i=1, 2, 3。

圖6 非協(xié)調(diào)重疊網(wǎng)格連接的多點(diǎn)位移約束法Fig.6 Multi-point constraint method for overlapping and nonconforming mesh

式中：Ni(xp,yp)為宏觀單元e在P點(diǎn)處的形函數(shù)(插值基函數(shù))值。

當(dāng)細(xì)觀結(jié)點(diǎn)P的位移滿足上述約束方程時(shí)，宏觀模型與細(xì)觀模型在該點(diǎn)處變形即是協(xié)調(diào)的。需要說明的是，雖然以上是以三結(jié)點(diǎn)三角形單元為例闡述通過基于形函數(shù)插值的多點(diǎn)位移約束實(shí)現(xiàn)宏細(xì)觀非協(xié)調(diào)重疊網(wǎng)格連接的方法，但該方法對(duì)宏觀單元的類型并無(wú)限制。對(duì)于其他類型的宏觀單元，僅需依據(jù)宏觀單元的位移模式調(diào)整式(13)～式(16)中的形函數(shù)表達(dá)式即可。

3 數(shù)值實(shí)現(xiàn)方法

3.1 數(shù)值求解流程

如圖7所示，由于在基于雙重網(wǎng)格的混凝土自適應(yīng)宏細(xì)觀協(xié)同有限元分析中涉及到細(xì)觀尺度下的材料非線性，故其數(shù)值求解宜采用增量迭代法。但由于在分析中涉及到源于宏細(xì)觀協(xié)同有限元模型動(dòng)態(tài)更新的變結(jié)構(gòu)非線性，故其數(shù)值求解流程與傳統(tǒng)的增量迭代法又有所區(qū)別，主要體現(xiàn)為對(duì)于任一增量步，均需在原平衡迭代的基礎(chǔ)上進(jìn)行“一致性”迭代，以保證該增量步求解完成后的有限元模型宏細(xì)觀分析區(qū)域劃分與應(yīng)力計(jì)算結(jié)果保持一致，即各宏觀單元任一積分點(diǎn)的收斂應(yīng)力解均應(yīng)不滿足如式(9)所示的分析尺度轉(zhuǎn)換準(zhǔn)則。此外，由于宏細(xì)觀協(xié)同有限元模型中細(xì)觀模型的位置與范圍是在分析過程中基于宏觀模型應(yīng)力計(jì)算結(jié)果自適應(yīng)確定的，故在開始第一個(gè)增量步分析時(shí)，假定整體有限元模型全部由宏觀模型構(gòu)成。

圖7 自適應(yīng)宏細(xì)觀協(xié)同有限元求解流程Fig.7 Flowchart of adaptive macro-meso-scale concurrent finite element solution

3.2 宏細(xì)觀協(xié)同有限元模型更新

如前所述，在自適應(yīng)宏細(xì)觀協(xié)同有限元分析過程中，需動(dòng)態(tài)更新宏細(xì)觀協(xié)同有限元模型以保證在細(xì)觀尺度下開展混凝土的損傷破壞分析。因此，對(duì)于任一需要進(jìn)行分析尺度轉(zhuǎn)換的宏觀單元，均要確定與該宏觀單元關(guān)聯(lián)的細(xì)觀單元集合。在宏觀網(wǎng)格和細(xì)觀網(wǎng)格獨(dú)立剖分的前提下，為保證用于替換某宏觀單元的細(xì)觀單元集合完全填充該宏觀單元占據(jù)的空間，細(xì)觀單元集合需包括細(xì)觀模型中全部或部分位于該宏觀單元邊界內(nèi)的所有細(xì)觀單元。具體而言，若某細(xì)觀單元的任一結(jié)點(diǎn)位于該宏觀單元內(nèi)，則該細(xì)觀單元即屬于用于替換該宏觀單元的細(xì)觀單元集合，如圖8所示。

圖8 與宏觀單元關(guān)聯(lián)的細(xì)觀單元集合Fig.8 Assemble of mesoscale elements associated with a macroscale element

遵循上述細(xì)觀單元集合確定原則，即可在某增量步的平衡迭代收斂后，通過在宏細(xì)觀協(xié)同有限元模型中將需要進(jìn)行分析尺度轉(zhuǎn)換的宏觀單元替換為與之相應(yīng)的細(xì)觀單元集合，完成宏細(xì)觀協(xié)同有限元模型更新。

3.3 多點(diǎn)位移約束方程定義

基于2.3節(jié)中提出的基于形函數(shù)插值的多點(diǎn)位移約束法，本文利用ABAQUS提供的多點(diǎn)約束(Multi-Point Constraint，MPC)功能[33]來(lái)實(shí)現(xiàn)宏細(xì)觀模型中不同尺度網(wǎng)格間的連接。具體而言，將位于某一宏觀單元內(nèi)部的細(xì)觀結(jié)點(diǎn)作為“從結(jié)點(diǎn)”，將該宏觀單元的結(jié)點(diǎn)作為“主結(jié)點(diǎn)”，并在獲取“從結(jié)點(diǎn)”與“主結(jié)點(diǎn)”坐標(biāo)的基礎(chǔ)上，計(jì)算出“從結(jié)點(diǎn)”位移約束方程(見2.3節(jié))的各個(gè)系數(shù)，從而確定以“主結(jié)點(diǎn)”位移為變量的“從結(jié)點(diǎn)”位移表達(dá)式并按約定格式在ABAQUS輸入文件中定義該細(xì)觀結(jié)點(diǎn)的多點(diǎn)位移約束方程，實(shí)現(xiàn)宏細(xì)觀協(xié)同有限元模型中非協(xié)調(diào)重疊網(wǎng)格的連接。以2.3節(jié)中的細(xì)觀結(jié)點(diǎn)P為例，給出了多點(diǎn)位移約束方程在ABAQUS中的定義格式，如圖9所示。

圖9 多點(diǎn)位移約束方程定義格式Fig.9 Definition format multi-point displacement constraint

4 算例分析

在上述基礎(chǔ)上，以ABAQUS為有限元求解工具，在MATLAB平臺(tái)上研發(fā)了基于雙重網(wǎng)格的混凝土自適應(yīng)宏細(xì)觀協(xié)同有限元分析程序ACMSC。為驗(yàn)證本文方法的可行性和程序編制的正確性，進(jìn)行如下算例分析。

算例1. 模擬了混凝土L形試件受拉損傷開裂過程。圖10(a)給出了試件尺寸、加載條件及邊界條件，并同時(shí)示出了Winkler等[34]通過物理試驗(yàn)獲取的宏觀裂縫分布范圍。圖10(b)給出了采用AutoGMC軟件生成的試件細(xì)觀結(jié)構(gòu)，骨料粒徑范圍為5 mm～20 mm，體積含量為50%。采用1.2節(jié)所述程序完成了細(xì)觀模型的有限元網(wǎng)格剖分，如圖10(c)所示，該圖中同時(shí)示出了宏觀模型的有限元網(wǎng)格，宏細(xì)觀模型的網(wǎng)格剖分均采用三結(jié)點(diǎn)三角形單元，宏觀模型單元數(shù)量為168個(gè)，細(xì)觀模型單元數(shù)量為37 423個(gè)。表1列出了細(xì)觀模型各相的材料參數(shù)。由于ITZ力學(xué)參數(shù)難以通過試驗(yàn)手段測(cè)得，通常認(rèn)為ITZ的力學(xué)性能與水泥砂漿的類似，參數(shù)取值略小于砂漿[2,3, 12, 19]。

圖10 L形試件算例及宏細(xì)觀有限元網(wǎng)格Fig.10 Numerical example of L-shape specimen and its macro-meso-scale mesh

表1 細(xì)觀材料參數(shù)Table 1 Mesoscale material parameters

為保證宏觀模型與細(xì)觀模型在彈性階段力學(xué)行為的一致性，開展如表1所示細(xì)觀材料參數(shù)下的混凝土單軸拉伸細(xì)觀數(shù)值試驗(yàn)(骨料粒徑范圍與體積含量與細(xì)觀模型相同，細(xì)觀計(jì)算模型如圖11(a)所示)，并基于數(shù)值試驗(yàn)所獲均勻化應(yīng)力-應(yīng)變曲線(見圖11(b))，取應(yīng)力從0～0.4ft的割線彈性模量為宏觀模型的彈性模量[35]，量值為28.6 GPa。此外，為確定自適應(yīng)宏細(xì)觀尺度轉(zhuǎn)換準(zhǔn)則參數(shù)C1、C2、C3和K的取值，亦通過開展單軸壓縮數(shù)值試驗(yàn)確定了單軸抗壓強(qiáng)度f(wàn)c(15.2 MPa)，進(jìn)而結(jié)合單軸拉伸數(shù)值試驗(yàn)確定的單軸抗拉強(qiáng)度f(wàn)t(1.45 MPa)，并取fb=1.16fc[10]，即可基于式(10)確定C1、C2、C3和K；應(yīng)力放大系數(shù)s取為1.25。

圖11 單軸拉伸細(xì)觀計(jì)算模型及應(yīng)力-應(yīng)變曲線Fig.11 Mesoscale model of uniaxial tension and the stress-strain curve

在上述基礎(chǔ)上，開展了混凝土L形試件受拉開裂的自適應(yīng)宏細(xì)觀協(xié)同有限元分析，位移荷載分為32個(gè)增量步逐級(jí)施加。此外，為對(duì)比驗(yàn)證分析成果的合理性，亦開展了相同條件下的全細(xì)觀模型數(shù)值模擬。圖12(a)～圖12(d)給出了自適應(yīng)宏細(xì)觀協(xié)同有限元分析所得的試件損傷開裂過程(為便于觀察，圖中隱去了損傷變量大于0.95的單元)，圖12(e)～圖12(h)給出了相應(yīng)的全細(xì)觀模型數(shù)值模擬結(jié)果。圖13對(duì)比了自適應(yīng)宏細(xì)觀協(xié)同有限元分析與全細(xì)觀模擬所得的加載邊界反力與加載位移關(guān)系曲線，其中，ACMSC和DNS分別表示自適應(yīng)宏細(xì)觀協(xié)同有限元分析和全細(xì)觀模型數(shù)值模擬所得的關(guān)系曲線。

圖12 混凝土L形試件受拉開裂過程Fig.12 Tension cracking process of concrete L-shape specimen

圖13 加載邊界反力-位移曲線Fig.13 Reaction force-displacement curve

從圖12中可以看出，在加載過程中，細(xì)觀損傷肇始于L形試件轉(zhuǎn)角處，繼而沿水平略偏上方向向試件內(nèi)部擴(kuò)展并逐漸形成宏觀裂縫，裂縫在試件內(nèi)的分布位于Winkler等[34]通過物理試驗(yàn)獲取的宏觀裂縫分布范圍內(nèi)(見圖10(a))；隨著加載位移的逐漸增大，宏觀分析區(qū)域逐漸減小，細(xì)觀分析區(qū)域逐漸增大，損傷開裂始終發(fā)生在細(xì)觀分析區(qū)域內(nèi)；在不同加載階段，自適應(yīng)宏細(xì)觀協(xié)同有限元分析所得的宏觀裂縫分布特征均非常接近全細(xì)觀模擬結(jié)果，但在宏觀裂縫端部，兩種方法所得的開裂區(qū)分布在細(xì)觀尺度上存在一定差異，原因主要在于上述兩種方法對(duì)在自適應(yīng)宏細(xì)觀協(xié)同有限元分析中分析尺度未轉(zhuǎn)化為細(xì)觀尺度的區(qū)域采用了不同尺度的分析模型(自適應(yīng)宏細(xì)觀協(xié)同分析為宏觀線彈性模型，而全細(xì)觀模擬為細(xì)觀模型)，故難以獲得完全一致的分析結(jié)果。此外，與宏觀裂縫分布特征非常接近相應(yīng)的是，自適應(yīng)宏細(xì)觀協(xié)同有限元分析與全細(xì)觀模擬所得的位移加載邊界上的反力(加載邊界上各結(jié)點(diǎn)豎向結(jié)點(diǎn)反力之和)與加載位移關(guān)系曲線亦基本重合(見圖13)，表明自適應(yīng)宏細(xì)觀協(xié)同有限元分析可以達(dá)到與全細(xì)觀模擬相當(dāng)?shù)木取?/p>

圖14給出了在位移荷載逐級(jí)增加過程中宏細(xì)觀協(xié)同有限元模型自由度數(shù)量的變化過程，為對(duì)比分析，亦示出了全細(xì)觀模型的自由度數(shù)量，可以看出，在加載初期，與全細(xì)觀模型相比，宏細(xì)觀協(xié)同有限元模型的計(jì)算自由度數(shù)量基本可忽略不計(jì)；隨著加載位移的逐漸增大，宏細(xì)觀協(xié)同有限元模型的計(jì)算自由度亦逐漸增加，完成加載時(shí)，宏細(xì)觀協(xié)同有限元模型的計(jì)算自由度約為全細(xì)觀模型的33.28%?？紤]到宏細(xì)觀協(xié)同有限元模型的計(jì)算自由度在加載過程中是逐步增加的，而全細(xì)觀模型的計(jì)算自由度在加載過程中保持不變，故在保證分析精度的前提下，本文方法的分析效率明顯高于全細(xì)觀模擬。

圖14 荷載逐級(jí)增加過程中自由度數(shù)量的變化Fig.14 Variation of the degree of freedom during the loading process

算例2. 模擬了混凝土簡(jiǎn)支梁三點(diǎn)彎曲試驗(yàn)，試件尺寸及加載條件如圖15所示，骨料粒徑范圍與體積含量、宏觀與細(xì)觀模型材料參數(shù)及計(jì)算參數(shù)取值與算例1保持一致，位移荷載為0.2 mm，分為50步逐級(jí)施加。此外，亦開展了相應(yīng)的全細(xì)觀模型模擬。圖16(a)～圖16(d)給出了自適應(yīng)宏細(xì)觀協(xié)同有限元分析所得的簡(jiǎn)支梁彎拉開裂過程，全細(xì)觀模型數(shù)值模擬結(jié)果如圖16(e)～圖16(h)所示。圖17和圖18分別對(duì)比了上述兩種方法所得的加載點(diǎn)反力與位移關(guān)系和逐級(jí)加載過程中自由度數(shù)量變化過程。

圖15 三點(diǎn)彎曲試件尺寸及加載條件 /mmFig.15 The size and loading conditions of three-point bending specimen

圖16 混凝土簡(jiǎn)支梁彎拉開裂過程Fig.16 Flexural-tensile cracking process of concrete simply supported beam

從圖16中可以看出，在加載過程中，細(xì)觀損傷首先出現(xiàn)于梁底跨中部位，繼而沿豎向向試件內(nèi)部擴(kuò)展并逐漸形成宏觀裂縫，隨著加載位移的增大，通過自適應(yīng)分析尺度轉(zhuǎn)換進(jìn)入細(xì)觀分析尺度的區(qū)域逐漸增大；在不同加載階段，自適應(yīng)宏細(xì)觀協(xié)同有限元分析所得的宏觀裂縫分布特征均非常接近全細(xì)觀模擬結(jié)果，且加載點(diǎn)反力-位移曲線亦基本重合(見圖17)。與全細(xì)觀模型相比，在加載初期，宏細(xì)觀協(xié)同有限元模型的計(jì)算自由度數(shù)量基本可忽略不計(jì)；雖然隨著加載位移的逐漸增大，宏細(xì)觀協(xié)同有限元模型的計(jì)算自由度數(shù)量會(huì)逐漸增加(見圖18)，但直至完成加載，宏細(xì)觀協(xié)同有限元模型的計(jì)算自由度數(shù)量?jī)H為全細(xì)觀模型的11.21%。上述分析表明，采用本文方法可在兼顧效率與精度的前提下，實(shí)現(xiàn)考慮細(xì)觀材料結(jié)構(gòu)的混凝土損傷開裂跨尺度演化過程自適應(yīng)分析。

圖17 加載點(diǎn)反力-位移曲線Fig.17 Reaction force-displacement curve

圖18 荷載逐級(jí)增加過程中自由度數(shù)量的變化Fig.18 Variation of the degree of freedom during the loading process

5 結(jié)論

準(zhǔn)確分析混凝土的損傷開裂跨尺度演化過程需要考慮其細(xì)觀材料結(jié)構(gòu)。本文在有限元法框架內(nèi)，提出了一種基于雙重網(wǎng)格的混凝土損傷開裂自適應(yīng)宏細(xì)觀協(xié)同分析方法，并在MATLAB平臺(tái)上研發(fā)了相應(yīng)的計(jì)算程序ACMSC。該方法的主要特點(diǎn)及優(yōu)點(diǎn)是：

(1) 通過在分析域內(nèi)布置獨(dú)立剖分的宏觀與細(xì)觀網(wǎng)格和建立相應(yīng)的宏觀與細(xì)觀有限元模型，避免了在分析過程中剖分細(xì)觀網(wǎng)格和建立細(xì)觀模型的困難。

(2) 通過提出從宏觀尺度至細(xì)觀尺度的分析尺度自適應(yīng)轉(zhuǎn)換準(zhǔn)則，實(shí)現(xiàn)了依據(jù)宏觀應(yīng)力計(jì)算結(jié)果的宏觀和細(xì)觀分析區(qū)域自適應(yīng)劃分。

(3) 通過提出基于形函數(shù)插值的多點(diǎn)位移約束方法，解決了宏細(xì)觀非協(xié)同重疊網(wǎng)格的連接問題。

(4) 通過形成包括彈性區(qū)宏觀模型和損傷區(qū)細(xì)觀模型的宏細(xì)觀協(xié)同有限元模型，實(shí)現(xiàn)了考慮細(xì)觀材料結(jié)構(gòu)的混凝土損傷開裂跨尺度演化過程分析。

(5) 與全細(xì)觀模擬結(jié)果的對(duì)比分析表明，本文方法可在保證分析精度的前提下，高效分析混凝土損傷開裂的跨尺度演化過程，為開展混凝土材料與結(jié)構(gòu)的精細(xì)化破壞分析提供了可行手段。