基于雙重網(wǎng)格的混凝土自適應(yīng)宏細(xì)觀協(xié)同有限元分析方法

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

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

混凝土是典型的隨機多尺度準(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)出明顯的隨機、非均勻、局部化與跨尺度特征。因此，準(zhǔn)確分析混凝土從細(xì)觀裂紋萌生、擴展、集聚至宏觀裂縫形成這一復(fù)雜過程需要考慮其細(xì)觀材料結(jié)構(gòu)[7 ? 8]。

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

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

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

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

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

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

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

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

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

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

為建立細(xì)觀有限元模型，需對混凝土細(xì)觀結(jié)構(gòu)進(jìn)行網(wǎng)格剖分，本文利用ABAQUS前處理模塊，通過MATLAB和PYTHON混合編程開發(fā)了混凝土細(xì)觀有限元網(wǎng)格自動剖分程序。具體而言，首先基于模擬區(qū)域和骨料的幾何信息，依據(jù)ABAQUS規(guī)定的編寫規(guī)則[18]，由程序自動編寫可被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)格剖分實例Fig.4 Example of mesoscale finite element meshing

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

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

式中：ω為偏心率，用于描述塑性勢函數(shù)向其漸近線逼近的速度，一般可取為0.1； σt0為單軸抗拉強度；ψ為膨脹角；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]。單軸受拉時，st=0 ，sc=1 ， 故d=dt；單軸受壓時，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)格獨立剖分，在剖分宏觀網(wǎng)格時，混凝土被視為均勻線彈性材料，而在剖分細(xì)觀網(wǎng)格時，則將混凝土視為由(粗)骨料、砂漿和界面過渡區(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ù)分析對象受力狀態(tài)變化動態(tài)更新，具體而言，當(dāng)某宏觀單元內(nèi)任一積分點的應(yīng)力滿足分析尺度自適應(yīng)轉(zhuǎn)換準(zhǔn)則(詳見節(jié)2.2)時，即需將該宏觀單元從協(xié)同模型中消除并激活與之相應(yīng)的細(xì)觀單元集合。

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

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

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

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

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

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

基于上述自適應(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)力計算結(jié)果不符，需要更新宏細(xì)觀協(xié)同有限元模型并重新進(jìn)行該增量步的迭代求解；反之，若不存在要進(jìn)行分析尺寸轉(zhuǎn)換的宏觀單元，則表明當(dāng)前模型的宏細(xì)觀區(qū)域劃分與應(yīng)力計算結(jié)果相符，可進(jìn)行下一個增量步的迭代求解。

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

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

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

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

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

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

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

3.1 數(shù)值求解流程

如圖7所示，由于在基于雙重網(wǎng)格的混凝土自適應(yīng)宏細(xì)觀協(xié)同有限元分析中涉及到細(xì)觀尺度下的材料非線性，故其數(shù)值求解宜采用增量迭代法。但由于在分析中涉及到源于宏細(xì)觀協(xié)同有限元模型動態(tài)更新的變結(jié)構(gòu)非線性，故其數(shù)值求解流程與傳統(tǒng)的增量迭代法又有所區(qū)別，主要體現(xiàn)為對于任一增量步，均需在原平衡迭代的基礎(chǔ)上進(jìn)行“一致性”迭代，以保證該增量步求解完成后的有限元模型宏細(xì)觀分析區(qū)域劃分與應(yīng)力計算結(jié)果保持一致，即各宏觀單元任一積分點的收斂應(yīng)力解均應(yīng)不滿足如式(9)所示的分析尺度轉(zhuǎn)換準(zhǔn)則。此外，由于宏細(xì)觀協(xié)同有限元模型中細(xì)觀模型的位置與范圍是在分析過程中基于宏觀模型應(yīng)力計算結(jié)果自適應(yīng)確定的，故在開始第一個增量步分析時，假定整體有限元模型全部由宏觀模型構(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é)同有限元分析過程中，需動態(tài)更新宏細(xì)觀協(xié)同有限元模型以保證在細(xì)觀尺度下開展混凝土的損傷破壞分析。因此，對于任一需要進(jìn)行分析尺度轉(zhuǎn)換的宏觀單元，均要確定與該宏觀單元關(guān)聯(lián)的細(xì)觀單元集合。在宏觀網(wǎng)格和細(xì)觀網(wǎng)格獨立剖分的前提下，為保證用于替換某宏觀單元的細(xì)觀單元集合完全填充該宏觀單元占據(jù)的空間，細(xì)觀單元集合需包括細(xì)觀模型中全部或部分位于該宏觀單元邊界內(nèi)的所有細(xì)觀單元。具體而言，若某細(xì)觀單元的任一結(jié)點位于該宏觀單元內(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 多點位移約束方程定義

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

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

4 算例分析

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

算例1. 模擬了混凝土L形試件受拉損傷開裂過程。圖10(a)給出了試件尺寸、加載條件及邊界條件，并同時示出了Winkler等[34]通過物理試驗獲取的宏觀裂縫分布范圍。圖10(b)給出了采用AutoGMC軟件生成的試件細(xì)觀結(jié)構(gòu)，骨料粒徑范圍為5 mm～20 mm，體積含量為50%。采用1.2節(jié)所述程序完成了細(xì)觀模型的有限元網(wǎng)格剖分，如圖10(c)所示，該圖中同時示出了宏觀模型的有限元網(wǎng)格，宏細(xì)觀模型的網(wǎng)格剖分均采用三結(jié)點三角形單元，宏觀模型單元數(shù)量為168個，細(xì)觀模型單元數(shù)量為37 423個。表1列出了細(xì)觀模型各相的材料參數(shù)。由于ITZ力學(xué)參數(shù)難以通過試驗手段測得，通常認(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ù)值試驗(骨料粒徑范圍與體積含量與細(xì)觀模型相同，細(xì)觀計算模型如圖11(a)所示)，并基于數(shù)值試驗所獲均勻化應(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ù)值試驗確定了單軸抗壓強度fc(15.2 MPa)，進(jìn)而結(jié)合單軸拉伸數(shù)值試驗確定的單軸抗拉強度ft(1.45 MPa)，并取fb=1.16fc[10]，即可基于式(10)確定C1、C2、C3和K；應(yīng)力放大系數(shù)s取為1.25。

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

在上述基礎(chǔ)上，開展了混凝土L形試件受拉開裂的自適應(yīng)宏細(xì)觀協(xié)同有限元分析，位移荷載分為32個增量步逐級施加。此外，為對比驗證分析成果的合理性，亦開展了相同條件下的全細(xì)觀模型數(shù)值模擬。圖12(a)～圖12(d)給出了自適應(yīng)宏細(xì)觀協(xié)同有限元分析所得的試件損傷開裂過程(為便于觀察，圖中隱去了損傷變量大于0.95的單元)，圖12(e)～圖12(h)給出了相應(yīng)的全細(xì)觀模型數(shù)值模擬結(jié)果。圖13對比了自適應(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)部擴展并逐漸形成宏觀裂縫，裂縫在試件內(nèi)的分布位于Winkler等[34]通過物理試驗獲取的宏觀裂縫分布范圍內(nèi)(見圖10(a))；隨著加載位移的逐漸增大，宏觀分析區(qū)域逐漸減小，細(xì)觀分析區(qū)域逐漸增大，損傷開裂始終發(fā)生在細(xì)觀分析區(qū)域內(nèi)；在不同加載階段，自適應(yīng)宏細(xì)觀協(xié)同有限元分析所得的宏觀裂縫分布特征均非常接近全細(xì)觀模擬結(jié)果，但在宏觀裂縫端部，兩種方法所得的開裂區(qū)分布在細(xì)觀尺度上存在一定差異，原因主要在于上述兩種方法對在自適應(yīng)宏細(xì)觀協(xié)同有限元分析中分析尺度未轉(zhuǎn)化為細(xì)觀尺度的區(qū)域采用了不同尺度的分析模型(自適應(yīng)宏細(xì)觀協(xié)同分析為宏觀線彈性模型，而全細(xì)觀模擬為細(xì)觀模型)，故難以獲得完全一致的分析結(jié)果。此外，與宏觀裂縫分布特征非常接近相應(yīng)的是，自適應(yīng)宏細(xì)觀協(xié)同有限元分析與全細(xì)觀模擬所得的位移加載邊界上的反力(加載邊界上各結(jié)點豎向結(jié)點反力之和)與加載位移關(guān)系曲線亦基本重合(見圖13)，表明自適應(yīng)宏細(xì)觀協(xié)同有限元分析可以達(dá)到與全細(xì)觀模擬相當(dāng)?shù)木取?/p>

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

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

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

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

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

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

圖17 加載點反力-位移曲線Fig.17 Reaction force-displacement curve

圖18 荷載逐級增加過程中自由度數(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平臺上研發(fā)了相應(yīng)的計算程序ACMSC。該方法的主要特點及優(yōu)點是：

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

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

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

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

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