混凝土是一种典型的多相非均质复合材料，其力学性能和破坏特征与其材料组成密切相关。为了解决宏观方法在分析混凝土开裂问题时难以准确捕捉结构裂缝信息的问题，学者们尝试从细观、微观或更低尺度不断探索新的分析方法。目前，混凝土细观分析方法已被广泛应用于特殊环境下混凝土开裂研究^[1]与混凝土非均匀性效应分析^[2]等多类问题。常用的细观模拟方法包括离散元模型^[3]、格构模型^[4]、随机断裂模型^[5]与内聚力模型^[6](Cohesive Zone Model，简称CZM)等。其中，CZM细观模型通过在基质单元间嵌入内聚力单元实现了对混凝土内部微裂缝扩展和闭合过程的模拟，在分析多裂缝问题上具有优越表现^[6]。此外，零厚度内聚力单元还可有效避免界面过渡区(Interfacial Transition Zone，简称ITZ)厚度过小引发的建模问题。

由于拉伸荷载下混凝土破坏机理较为明确，且易于采用牵引力-分离准则描述其损伤演化，CZM细观模型起初主要应用于模拟混凝土的拉伸断裂^{[7− 9]}。REN等^[7]和WANG等^{[8 − 9]}分别在二维CT扫描细观模型和二维随机骨料细观模型中嵌入了零厚度内聚力单元，并分析了拉伸荷载作用下CZM细观模型的破坏过程和宏观力学性能。为了模拟压缩荷载下混凝土复杂的破坏模式，常采用的方法为定义混合断裂模式与颗粒间相互作用。WU等^[10]在CZM模型中引入Benzeggagh-Kenane(B-K)混合断裂判定准则以模拟骨料互锁和表面摩擦等复杂机制。WANG等^[11]在细观模型中引入混凝土损伤塑性模型(CDP)和内聚力单元并建立了CDP-CZM细观模型，分析了CDP的膨胀角参数对细观模型的影响。除定义B-K准则外，XIONG等^[12]在CZM模型中定义了接触行为以模拟颗粒间的碰撞和摩擦，取得了良好的应用效果，并且指出接触行为决定了模型的受压表现。

基于混合断裂模式的CZM细观模型可有效描述混凝土在复杂受力状态下的破坏模式，在混凝土结构细观分析方面具有突出的应用潜力。然而，模型需要定义大量诸如单元刚度、强度与断裂能等内聚力单元参数，且各内聚力单元参数对细观模型整体效应所产生的影响尚未明确，严重限制了CZM细观模型在结构分析中的应用。为此，在验证基于Python脚本生成的CZM细观模型网格敏感性和有效性的基础上，探讨并分析了混合断裂模式下关键内聚力单元参数对混凝土单轴受力性能的影响。

1   内聚力单元细观模型建立方法

1.1   生成随机骨料模型

采用Monte Carlo方法生成随机骨料细观模型，即按预设顺序随机投放各个骨料直至满足级配要求。为了实现这一过程，编写了Python脚本辅助骨料生成与投递，并将生成的骨料信息储存以便导入至Abaqus软件完成模型建立。为了使细观模型具有较优的密实度和宏观强度，假定各级粒径的骨料分布满足WALRAVEN公式^[13]，如式(1)所示。

式中：P_c为当前粒径组的骨料含量占比；P_k为粗骨料含量总占比；t=d_c/D_max表示当前粒径组代表直径与粗骨料最大直径的比值；d_c为当前粒径组的直径；D_max为所有骨料中的最大直径。

二维混凝土细观模型中，常用的骨料形状为圆形、椭圆形与多边形等。其中，圆形骨料多用于轻骨料混凝土细观模型或简化的普通混凝土细观模型^[1]，椭圆形骨料多用于模拟鹅卵石粗骨料^[8]，而多边形骨料则适合模拟如碎石等类型的粗骨料^[12]。综合考虑各粗骨料形状应用场景，将粗骨料形状设置为多边形，并基于基础圆形骨料完成任意多边形骨料生成工作，如图1所示。具体步骤如下所示：

1)将混凝土内粗骨料粒径级配分成[d₀, d₁]，…，[d_n−1, d_n]。根据Walraven公式求解各粒径级配在混凝土中的占比与面积；

2)随机生成满足粒径级配要求的骨料，并检验当前骨料是否存在边界冲突问题，若满足要求将骨料半径与圆心坐标位置[R_k, x_k, y_k]存入骨料信息库；

3)定义多边形顶点数量、角度以及长度参数，根据给定信息随机生成每个多边形骨料的顶点数j、顶点与x轴夹角α_j以及顶点到圆心的距离l_j；

4)在圆心坐标(x_k, y_k)和顶点方位点(x_k + R_k × cosα_j, y_k + R_k × sinα_j)的连线中，按照预设长度l_j计算多边形的各个顶点坐标。在笛卡尔坐标系中，不同顶点坐标可表示为(x_k + l_j × cosα_j, y_k + l_j × sinα_j)；

5)利用式(2)计算随机生成的第i级骨料的总面积A^ith，若满足级配占比要求则进行下一级骨料生成，否则重复上述步骤直至满足当前级配粗骨料占比要求；

6)保存全部的多边形骨料信息，并在Abaqus中利用Python脚本对细观模型进行参数化建模。

图  1  多边形骨料生成与面积计算

Figure  1.  Polygon aggregate generation and area calculation

下载: 全尺寸图片 幻灯片

1.2   生成零厚度内聚力单元

编写Python脚本以实现在随机骨料模型中自动化插入内聚力单元的过程。CZM细观模型的生成主要包含两个阶段，即基础模型信息读取阶段与新模型生成阶段。首先，利用循环命令读取原模型的单元编号与节点编号并保存。随后，根据初始模型节点和单元信息生成新的内聚力单元编号及对应的节点，如图2所示。按照程序在砂浆单元-砂浆单元和砂浆单元-骨料单元(ITZ层)之间插入新的内聚力单元。

图  2  嵌入零厚度内聚力单元

Figure  2.  Embedding zero thickness cohesive elements

下载: 全尺寸图片 幻灯片

1.3   内聚力单元本构关系

CZM细观模型的宏观力学性能及破坏模式由内聚力单元属性控制。混凝土裂缝扩展形态可分为Ⅰ型(张开型裂缝)、Ⅱ型(滑移型裂纹)以及Ⅲ型(撕裂型裂纹)，如图3所示。部分研究人员也将Ⅱ型断裂模式和Ⅲ型断裂模式统称为Ⅱ型，并分别称呼为第一剪切方向和第二剪切方向^[14]。混凝土在受拉和受压状态下的破坏分别以Ⅰ型断裂和复合型断裂为主，且其本身的非均质性和多相性会影响断裂模式。因此，在CZM细观模型中需要定义的内聚力单元属性包括单元刚度、Ⅰ型与Ⅱ型强度和断裂能以及混合损伤模式等参数。

图  3  裂缝扩展的三种形态

Figure  3.  Three modes of crack propagation

下载: 全尺寸图片 幻灯片

内聚力单元本构采用双线性牵引-分离准则描述，假设初始阶段为线弹性而随后是损伤的起始与演化过程，如图4所示^[14]。单元初始损伤采用平方名义应力准则来判定，如式3所示。

式中：t_n、t_s和t_t分别为受拉方向和两个剪切方向上的名义应力；〈 〉为Macaulay符号，当处于受压状态时其数值为0；$t_n^0$、$t_s^0$和$t_t^0$分别为在受拉方向和两个剪切方向上的最大允许名义应力。

图  4  牵引力-分离准则示意图

Figure  4.  Illustration of Traction-Separation Criterion

下载: 全尺寸图片 幻灯片

在双线性牵引-分离准则中引入损伤软化行为表征本构曲线的退化阶段。引入损伤变量D描述内聚力单元的损伤起始状态到最终破坏状态之间的演化性质，其数学表达式为式4^[14]。

式中：$\delta _m^{\max }$为加载过程中所获得的有效位移最大值；$\delta _m^0$为单元应力达到抗拉强度时的位移值；$\delta _m^f$对应的是材料断裂时的等效位移。

在CZM细观模型中引入接触行为与B-K损伤断裂准则以模拟复杂应力状态下材料细部特征，如图5所示。其中，法向接触行为定义为硬接触，切向接触行为定义了与摩擦角相关的摩擦系数。采用B-K损伤准则描述复杂应力状态下混凝土的混合破坏模式，如式(5)所示^[14]。

式中：G_c、$G_n^c$与$G_s^c$分别为混合型断裂能、Ⅰ型断裂能以及Ⅱ型断裂能；G_N、G_S分别为拉伸和剪切变形能；η为材料参数。

图  5  B-K损伤准则

Figure  5.  B-K damage criterion

下载: 全尺寸图片 幻灯片

2   细观模型破坏模式及参数选定

2.1   细观模型参数设置

参照《混凝土物理力学性能试验方法标准》(GB/T 50081−2019)^[15]中混凝土单轴抗压强度标准试件建立了尺寸为300 mm×150 mm的二维CZM混凝土细观模型。利用Python脚本在Abaqus软件中生成了粗骨料占比为49.86%的CZM细观模型，如图6(a)、图6(b)所示。砂浆单元与骨料单元均采用CPS3单元模拟，并假定二者均为理想弹性材料。内聚力单元采用COH2D4单元模拟，假定ITZ层内聚力单元的初始损伤强度和断裂能参数为一般内聚力单元的1/4，材料属性如表1所示。

在进行轴向压(拉)力加载前，对数值模型的一条短边施加由加载速度控制的竖向位移边界条件，并固定另外一条短边的竖向位移与中间节点的横向位移，如图6(c)所示。为了满足静力加载条件限制，在Abaqus/Explicit分析步中严格控制加载速率。根据CEB-FIB MC2010^[17]中的相关规定，压缩加载速率控制为3×10⁻⁵ s⁻¹，而拉伸加载速率则为3×10⁻⁶ s⁻¹。

图  6  CZM细观模型与边界条件

Figure  6.  CZM meso scale model and boundary condition

下载: 全尺寸图片 幻灯片

表  1  内聚力单元材料属性设置

Table  1.  Material property of cohesive elements

参数	砂浆	骨料	内聚力单元
密度	2300 kg/m3	2700 kg/m3	2300 kg/m3
弹性模量	2.5×104 MPa[16]	7.2×104 MPa[6]	−
泊松比	0.2[16]	0.2[16]	−
拉伸刚度	−	−	6×106 MPa[16]
剪切刚度	−	−	6×106 MPa[16]
Ⅰ型强度	−	−	5 MPa
Ⅰ型断裂能	−	−	10 MPa
Ⅱ型强度	−	−	0.1 N/mm
Ⅱ型断裂能	−	−	2.5 N/mm
材料参数	−	−	1.2[6]

下载: 导出CSV 

| 显示表格

2.2   网格尺寸敏感性分析

现有研究表明^{[11 − 12]}，混凝土细观模型的计算精度与其网格尺寸密切相关。尤其是压缩荷载作用下的CZM细观模型中，网格尺寸决定了内聚力单元的数量及分布，更进一步地影响开裂位置与裂缝发展形态。为此，针对同一细观模型分别按照0.5 mm、1.0 mm、1.5 mm、2.0 mm和3.0 mm的尺寸进行网格划分，计算并分析了不同网格密度的CZM模型在轴压、轴拉荷载下的整体响应。五组模型的节点数量、单元数量及计算分析时长如表2所示。

表  2  模型复杂度与计算时长

Table  2.  Model complexity and computation time

网格尺寸/ mm	节点数量	平面单元 数量	内聚力单元 数量	拉伸计算 时长/h	压缩计算 时长/h
0.5	726 709	397 656	320 955	56.98	58.77
1.0	193 177	102 258	87 149	19.62	20.03
1.5	92 886	47 436	42 665	13.56	13.94
2.0	59 754	27 957	29 080	8.53	8.51
3.0	43 334	19 266	20 837	6.04	6.12

下载: 导出CSV 

| 显示表格

网格尺寸对CZM细观模型的整体力学响应如图7所示。在轴拉荷载作用下，不同网格尺寸的细观模型的应力-应变曲线基本相同，仅在加载后期有细小差别。这是由于网格尺寸并不影响结构薄弱层(ITZ层)的分布，内聚力单元数量对结构薄弱层裂缝生成与开展基本无影响，整体表现为模型不依赖于网格尺寸。然而，压缩荷载作用下细观模型表现出了网格尺寸依赖性。0.5 mm、1.0 mm和1.5 mm网格尺寸下的细观模型峰值应力分别为28.38 MPa、28.02 MPa与28.12 MPa，基本维持在同一水平。而2.0 mm和3.0 mm网格尺寸模型承载力计算值分别为27.02 MPa和26.18 MPa，小于上述模型。其中，粗糙网格模型和精细网格模型的峰值应力相差约8.4%。因此，在充分考虑模型准确性与时间成本的基础上，选取CZM模型的网格尺寸为1.5 mm。

图  7  CZM细观模型网格敏感性

Figure  7.  Mesh sensitivity of CZM mesoscale model

下载: 全尺寸图片 幻灯片

2.3   轴向荷载作用下的细观模型破坏形式

图8展示了轴向拉伸和压缩作用下混凝土的应力-应变曲线和破坏过程。由于混凝土力学性能与粒径级配和骨料分布等因素息息相关，选用平均曲线代表该组参数下的细观模型响应。

图  8  单轴应力-应变曲线与破坏特征

Figure  8.  Axial tensile stress-strain curve and failure characteristics

下载: 全尺寸图片 幻灯片

根据CZM细观模型裂缝开展特征，可将轴向加载过程分为三个阶段。第一阶段为弹性加载阶段(OA段)，对于两个不同的加载方向，内聚力单元几乎全部处于弹性阶段，仅有个别薄弱层的内聚力单元达到损伤起始状态，表现为试件的刚度几乎没有降低。随着位移逐渐增大，CZM模型进入第二加载阶段(AB段)。对于拉伸模型，模型出现一条明显的薄弱带，该处的内聚力单元达到极限状态的数量逐渐增多，并且混凝土裂缝逐渐沿横向发展。由于骨料-砂浆界面的内聚力单元属性明显弱于其他内聚力单元，横向开裂区中达到极限状态的单元也多为骨料-砂浆界面内聚力单元，其余位置未见明显损伤。在轴压模型中，由于竖向位移增大，大量内聚力单元达到损伤条件，部分骨料-砂浆界面内聚力单元达到极限状态，混凝土内部产生大量的细微裂缝。在应力退化阶段(BC段)，拉伸试件中部裂缝逐渐变宽、变长，并最终形成一条横向贯穿的裂缝。由于实体单元之间的传力路径被删除，细观模型承载能力降低。压缩荷载作用下，随着内聚力单元的退化与删除，CZM细观模型内部形成一条明显的斜向裂缝与多条微裂缝。随着竖向荷载逐渐增大，裂缝宽度明显增加，试件底部出现明显的横向位移。

3   参数效应分析

3.1   内聚力单元刚度

在牵引-分离准则中，单元刚度和强度决定了内聚力单元的线弹性行为，而强度和断裂能则定义了单元损伤的起始和演化。因此，刚度参数对内聚力单元表现和模型宏观力学性能影响显著。由于细观模型中内聚力单元的物理厚度为零，在某种程度上单元刚度发挥了罚参数的作用^[7]。不合适的单元刚度将造成收敛困难或结果奇异等问题，影响模型计算精度和效率。为此，分别在三组细观模型参数中选取了5×10⁴ MPa/mm、1×10⁵ MPa/mm、5×10⁵ MPa/mm、1×10⁶ MPa/mm、5×10⁶ MPa/mm、1×10⁷ MPa/mm共六个不同刚度数据，分析单元刚度对细观模型力学性能的影响。内聚力单元计算厚度采用软件默认数值1 mm。

在参数组1中，内聚力单元刚度对CZM细观模型影响如图9所示。单元刚度对CZM细观模型的初始刚度、极限强度及其对应的应变值等参数均有不同程度影响。可以明显观察到，随着内聚力单元刚度的增大，模型的初始刚度和抗压强度逐渐增大，而模型的抗拉强度基本不受刚度改变的影响。对比模型强度所对应的特征应变可知，无论是轴拉和轴压荷载，单元刚度对模型特征应变均有影响。由于单元刚度的罚参数作用，过大或过小的单元刚度会导致应力-应变曲线部分奇异。因此，选择合适的内聚力单元刚度十分重要。

图  9  内聚力单元刚度对CZM细观模型影响

Figure  9.  Influence of cohesive element stiffness on CZM Model

下载: 全尺寸图片 幻灯片

参数组二、三条件下的细观模型在轴向荷载作用下的强度和特征应变计算结果如图10所示。在保持其他参数不变的条件下，混凝土单轴抗压强度随着内聚力单元刚度的增加而变大，基本成指数关系。而对于单轴抗拉强度，基本可以认为单元刚度对其没有影响。究其原因，受拉阶段破坏强度主要取决于模型薄弱层内聚力单元强度，改变单元刚度对抗拉强度的影响十分有限。单元刚度对拉伸、压缩特征应变的影响显著，在单元刚度接近1×10⁶ MPa/mm时，曲线形成了拐点。当单元刚度小于该数值时，细观模型特征应变下降迅速。当高于该数值时，特征应变基本趋于平稳。在两个加载方向上，CZM细观模型的初始刚度基本相同，且表现出了与内聚力单元刚度近似线性的关系。

图  10  内聚力单元刚度对CZM细观模型关键参数的影响

Figure  10.  Influence of cohesive element stiffness on critical parameters of CZM microscale model

下载: 全尺寸图片 幻灯片

3.2   Ⅰ型强度

图11展示了在保持其他参数不变、单独改变Ⅰ型强度的条件下，各个数值模型的应力-应变曲线。根据内聚力单元双折线本构关系，在保持单元刚度相同的前提下，单元强度将决定本构曲线上升段。由模拟结果可知，Ⅰ型强度对数值模型的抗拉和抗压强度均有较大影响。这也表明在混合断裂模式下Ⅰ型断裂模式同时影响压缩和拉伸作用下CZM细观模型的力学响应。同时，由于单元强度减小，为了保持本构曲线包围的面积相同，即断裂能保持一致，内聚力单元的极限位移值将更难达到。这也解释了为什么细观模型应力应变曲线的下降段有明显差别。

图  11  Ⅰ型强度对CZM细观模型影响

Figure  11.  Influence of Mode I nominal stress on CZM microscale model

下载: 全尺寸图片 幻灯片

图12展示了各参数组的模拟结果。Ⅰ型强度对CZM细观模型的抗拉强度、抗压强度和特征应变均有影响，且可近似视为线性关系。然而，Ⅰ型强度对数值模型的初始刚度基本没有影响。这是由于细观数值模型的初始刚度由骨料单元刚度、砂浆单元刚度和内聚力单元控制，基本不受其他因素影响。

图  12  Ⅰ型强度对CZM细观模型关键参数的影响

Figure  12.  Influence of Mode I nominal stress on critical parameters of CZM microscale model

下载: 全尺寸图片 幻灯片

3.3   Ⅱ型强度

Ⅱ型强度对CZM混凝土细观模型的影响如图13所示。当改变Ⅱ型强度时，数值模型的单轴抗压强度随参数增大而增大，应力-应变曲线下降段随参数增大而变陡。曲线下降段趋势变化的原因与单独改变Ⅰ型强度的原因基本一致。Ⅱ型强度对细观模型抗拉强度影响较小，归因于拉伸荷载作用下细观模型主要以Ⅰ型断裂为主，Ⅱ型断裂参数对拉伸表现影响有限。整体而言，Ⅱ型强度对抗压强度、抗拉强度和特征应变的影响近似线性关系。

图  13  Ⅱ型强度对CZM细观模型关键参数的影响

Figure  13.  Influence of Mode Ⅱ nominal stress on critical parameters of CZM microscale model

下载: 全尺寸图片 幻灯片

3.4   Ⅰ型断裂能

如图14所示，Ⅰ型断裂能对CZM细观模型抗压强度及特征应变几乎没有影响。对于参数组1，模型抗压强度最大值和最小值的差距仅为3.05%。对于轴拉响应，改变Ⅰ型断裂能参数后，模型的应力-应变曲线在达到极限抗拉强度前基本吻合。降低Ⅰ型断裂能后，模型抗拉强度明显降低。这是由于模型薄弱层处内聚力单元的删除条件会随断裂能降低而更易达到，即潜在断裂层更易发生断裂。Ⅰ型断裂能对模型整体响应参数的影响也基本可视为线性关系。

图  14  Ⅰ型断裂能对CZM细观模型关键参数的影响

Figure  14.  Influence of Mode I fracture energy on critical parameters of CZM microscale model

下载: 全尺寸图片 幻灯片

3.5   Ⅱ型断裂能

与Ⅰ型断裂能对模型响应影响相似，Ⅱ型断裂能对模型应力-应变曲线上升段影响很小，基本只影响极限强度和下降段变化趋势。在图15中，可以明显观察到Ⅱ型断裂对数值模型强度和特征应变的影响呈线性关系。不同数值模型的抗拉强度最大相差约1.40%，抗压强度最大相差约21.4%。表明Ⅱ型参数对模型轴向拉伸荷载响应基本没有影响。同时，Ⅱ型断裂能对曲线下降段影响显著，断裂能参数组在控制曲线下降段方面至关重要。通过对比数值模拟结果与规范推荐本构关系曲线可知，Ⅱ型断裂能推荐为Ⅰ型断裂能的25倍。

图  15  Ⅱ型断裂能对CZM细观模型关键参数的影响

Figure  15.  Influence of Mode Ⅱ fracture energy on critical parameters of CZM microscale model

下载: 全尺寸图片 幻灯片

4   结论

基于Python脚本在Abaqus软件中自动生成了基于零厚度内聚力单元的混凝土细观分析模型，在验证了CZM细观模型有效性和网格尺寸依赖性的基础上，针对重要的内聚力单元参数开展了参数研究工作。研究结论主要如下：

(1)在轴向压缩和拉伸荷载作用下，CZM细观模型的初始刚度主要由砂浆单元刚度、骨料单元刚度以及内聚力单元刚度控制。在保证实体单元属性不变的条件下，初始刚度仅与内聚力单元刚度相关，基本不受强度和断裂能参数的影响。此外，为了满足计算精度和效率要求，需要选择合适的内聚力单元刚度。

(2)细观模型极限强度与特征应变受内聚力单元刚度、强度和断裂能影响。抗压强度与单元刚度呈指数关系，而抗拉强度不受单元刚度影响。极限强度均与单元强度和断裂能之间为线性关系。特征参数与单元刚度之间呈指数关系，而与强度和断裂能之间大致可描述为线性关系。

(3)在混合断裂模式下，轴压作用下的CZM混凝土细观模型响应受Ⅱ型剪切断裂模式和Ⅰ型张开断裂模式共同作用。而轴拉作用则主要受Ⅰ型张开断裂模式参数控制，受Ⅱ型参数影响有限。

(4)内聚力单元断裂能参数对CZM细观模型的单轴应力-应变曲线下降段影响十分明显。为了获得匹配真实混凝土的本构曲线，应当先根据拉伸响应确定Ⅰ型断裂能，后根据轴压曲线确定Ⅱ型断裂能。对比数值模拟结果和规范推荐混凝土本构关系，Ⅱ型断裂能应设置为Ⅰ型断裂能的25倍。