基于双重网格的混凝土自适应宏细观协同有限元分析方法

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

(河海大学水利水电学院，江苏，南京 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) 与全细观模拟结果的对比分析表明，本文方法可在保证分析精度的前提下，高效分析混凝土损伤开裂的跨尺度演化过程，为开展混凝土材料与结构的精细化破坏分析提供了可行手段。