基于雙重網格的混凝土自適應宏細觀協同有限元分析方法

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

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

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

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

圖1 結構分區與宏細觀協同計算模型示意圖Fig.1 Sketches of structure decomposition and macro-meso-scale concurrent computational model

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

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

1 混凝土細觀有限元模型

為了應用有限單元法開展細觀尺度下的混凝土損傷開裂分析，首先需要建立混凝土細觀有限元模型，主要涉及細觀結構模擬、有限元網格剖分和細觀力學模型等3個方面，分述如下。

1.1 細觀結構模擬

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

圖2 細觀結構生成實例Fig.2 Examples of meso-structure generation

1.2 有限元網格剖分

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

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

圖4 細觀有限元網格剖分實例Fig.4 Example of mesoscale finite element meshing

1.3 細觀力學模型

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

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

CDP模型采用如下形式的屈服函數：

引入拉伸、壓縮損傷因子dt、dc分別表征拉伸、壓縮損傷導致的剛度退化，其量值分別隨拉伸、壓縮等效塑性應變的變化而變化。進一步考慮應力反向后的剛度恢復效應，即可給出復雜應力狀態下d與dt、dc之間的關系式：

式中：st、sc的取值與應力狀態相關[24]。單軸受拉時，st=0 ，sc=1 ， 故d=dt；單軸受壓時，sc=0，st=1 ， 故d=dc。

2 宏細觀協同有限元分析方法

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

2.1 基于雙重網格的宏細觀協同有限元模型

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

如圖5(b)所示，在宏細觀協同分析中，僅有部分宏觀模型被作為整體模型的一部分，其余部分則被替換為與之相應的細觀模型，從而形成宏細觀協同分析整體有限元模型；上述協同模型需依據分析對象受力狀態變化動態更新，具體而言，當某宏觀單元內任一積分點的應力滿足分析尺度自適應轉換準則(詳見節2.2)時，即需將該宏觀單元從協同模型中消除并激活與之相應的細觀單元集合。

圖5 宏細觀協同有限元模型Fig.5 Macro-meso-scale concurrent finite element model

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

2.2 自適應尺度轉換

混凝土結構的不均勻應力分布與混凝土材料的應變軟化特性決定了在混凝土結構宏細觀協同有限元分析中，僅需對部分區域(損傷區)開展細觀尺度分析[30]。但由于實際混凝土結構受力狀態的復雜性，通常無法在分析前準確確定損傷區的位置與范圍[31 ? 32]，故需在分析過程中依據結構當前受力狀態確定需要將分析尺度從宏觀轉換為細觀的區域并動態更新宏細觀協同有限元模型，這一過程即為分析尺度的自適應轉換。

為在分析過程中實現分析尺度的自適應轉換，本文基于Ottosen多軸強度準則[16]，提出了以積分點應力為指標的混凝土自適應宏細觀尺度轉換準則，如下式所示：

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

基于上述自適應宏細觀尺度轉換準則，即可在某一增量步迭代收斂后，依據各宏觀單元的當前應力狀態判斷是否存在需要進行分析尺寸轉換的宏觀單元，若存在，則表明當前宏細觀協同有限元模型的宏細觀區域劃分與應力計算結果不符，需要更新宏細觀協同有限元模型并重新進行該增量步的迭代求解；反之，若不存在要進行分析尺寸轉換的宏觀單元，則表明當前模型的宏細觀區域劃分與應力計算結果相符，可進行下一個增量步的迭代求解。

2.3 非協調重疊網格連接

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

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

圖6 非協調重疊網格連接的多點位移約束法Fig.6 Multi-point constraint method for overlapping and nonconforming mesh

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

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

3 數值實現方法

3.1 數值求解流程

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

圖7 自適應宏細觀協同有限元求解流程Fig.7 Flowchart of adaptive macro-meso-scale concurrent finite element solution

3.2 宏細觀協同有限元模型更新

如前所述，在自適應宏細觀協同有限元分析過程中，需動態更新宏細觀協同有限元模型以保證在細觀尺度下開展混凝土的損傷破壞分析。因此，對于任一需要進行分析尺度轉換的宏觀單元，均要確定與該宏觀單元關聯的細觀單元集合。在宏觀網格和細觀網格獨立剖分的前提下，為保證用于替換某宏觀單元的細觀單元集合完全填充該宏觀單元占據的空間，細觀單元集合需包括細觀模型中全部或部分位于該宏觀單元邊界內的所有細觀單元。具體而言，若某細觀單元的任一結點位于該宏觀單元內，則該細觀單元即屬于用于替換該宏觀單元的細觀單元集合，如圖8所示。

圖8 與宏觀單元關聯的細觀單元集合Fig.8 Assemble of mesoscale elements associated with a macroscale element

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

3.3 多點位移約束方程定義

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

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

4 算例分析

在上述基礎上，以ABAQUS為有限元求解工具，在MATLAB平臺上研發了基于雙重網格的混凝土自適應宏細觀協同有限元分析程序ACMSC。為驗證本文方法的可行性和程序編制的正確性，進行如下算例分析。

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

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

表1 細觀材料參數Table 1 Mesoscale material parameters

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

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

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

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

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

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

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

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

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

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

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

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

5 結論

準確分析混凝土的損傷開裂跨尺度演化過程需要考慮其細觀材料結構。本文在有限元法框架內，提出了一種基于雙重網格的混凝土損傷開裂自適應宏細觀協同分析方法，并在MATLAB平臺上研發了相應的計算程序ACMSC。該方法的主要特點及優點是：

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

(2) 通過提出從宏觀尺度至細觀尺度的分析尺度自適應轉換準則，實現了依據宏觀應力計算結果的宏觀和細觀分析區域自適應劃分。

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

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

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