自2000年以来，中国基础设施建设飞速发展，大跨度桥梁与海港等混凝土基础设施建筑不断增多。根据《2022年交通运输行业发展统计公报》，截至2022年底，全国公路桥梁103.32万座、总长度为8 576.49万m，桥梁数量连续多年位居世界第一，其中混凝土结构桥梁达到90%以上[1]。但随着混凝土的大量应用，其自身的裂缝问题逐步凸显。同时，混凝土裂缝带来的结构破坏日益严重，混凝土维修成本居高不下[2]。随着大规模混凝土结构的加固和维修，研究者对混凝土开裂的研究更加深入。从细观尺度研究混凝土结构裂缝萌生和发展机理，优化混凝土抗裂性能具有重要的工程意义。随着计算机技术的发展，对混凝土的相关研究已经从传统的宏观分析发展到细观、微观分析[3-4]。在宏观尺度上，混凝土通常被认为是各向均质材料，其应力与应变关系通常呈近似线性关系，进而忽略其细观结构对两者关系的非线性影响[5]。对传统的混凝土力学性能分析仅仅停留在对裂纹的定性分析。但随着混凝土应用场景的多样化、复杂化，实际工程对于混凝土的性能要求进一步提高，对混凝土力学性能的相关研究不能仅仅停留在宏观层次[6-7]。在细观尺度下，混凝土通常认为是由砂浆、骨料与这两者过渡区构成[8-9]，这三部分力学性质的差异导致了混凝土力学特性的非均质性与裂纹的随机性[10]。
对混凝土进行细观断裂数值模拟时，首先需要建立细观尺度的混凝土随机骨料模型。目前，国内外学者提出了多种细观混凝土数值模型的构建思路。高郑国等[11-12]提出了二维混凝土随机骨料模型和三维凸性混凝土骨料随机投放算法。秦川等[13]提出了基于背景网格的预处理方法，提高了骨料含量与投放效率。ZHANG等[14]提出了基于缩放骨料的处理手段，可以同时生成包含凸骨料与凹骨料的随机模型。随机凸多边形骨料的生成主要有圆形基质法[15]和极坐标法[16]。圆形基质法主要通过在圆上随机分布顶点，随后连线形成多边形骨料，该方法能够直接控制骨料粒径。极坐标法通过随机控制点的极径和极角来确定端点的位置，无法直接控制骨料粒径，但生成效率较高。同时，影响骨料生成效率的另外一个重要因素就是干涉判断。常见的干涉判断方法有中心距法[15]、面积标度法[11]和点阵法[13]。中心距法通过圆形基质的圆心距离来判断骨料是否发生重叠，该方法限制了骨料的可投放区间，导致较大体积分数骨料投放时间过长且难以成功。面积标度法则难以对特殊重叠情况的骨料进行判断，且运算量较大。点阵法在一定程度上增加了骨料的投放面积，但往往需要大量的基点数据进行判断，运行效率低。因此，这几种干涉判断方法仍需进一步改进。
目前，国内有较多学者利用内聚力单元来研究细观层次混凝土的断裂行为。林力等[17]基于内聚力模型研究了粉煤灰掺量和孔隙率对混凝土开裂的影响。熊学玉等[18]将混凝土内聚力模型的适用范围扩大到受压断裂领域，较好地模拟了受压情况下混凝土的断裂行为。李坤罡[19]研究了骨料参数以及界面力学参数对混凝土力学性能的影响。张如毅[20]分析了内聚力模型断裂参数的敏感性。谢浩等[21-22]在混凝土细观模型中引入孔隙结构和畸形骨料，并分析了其对混凝土物理力学性能的影响。张龙飞等[23]通过激光扫描的方法生成三维骨料并对骨料进行投放。王杨[24]研究了速率相关作用和动荷载作用下混凝土的强度、变形及破坏模式的变化。吴昊等[25]研究了低温条件下再生沥青混凝土的断裂性能。虽然学者对混凝土细观断裂影响因素进行了探讨，但并未提出混凝土细观结构的变异性模型。为此，本文首先对二维多边形随机骨料的投放算法进行改进，以提高骨料的投放体积和投放效率。同时，在内聚力模型的基础上编制全局内聚力单元插入程序，以便在任意实体单元之间插入内聚力单元并对内聚力单元进行分类，克服传统的内聚力单元仅能在预定路径与单元之间进行插入的局限。最后，基于全局插入内聚力模型，研究三相混凝土细观断裂行为，分析骨料体积分数和加载速率等对混凝土受拉本构模型的影响，并基于大量模拟结果建立混凝土受拉概率本构模型。
1 细观模型的建立
1.1 二维随机多边形骨料模型
随机骨料模型由水泥砂浆、骨料以及ITZ过渡区构成。混凝土骨料通常分为卵石骨料与碎石骨料。对卵石骨料通常采用圆形来模拟，碎石骨料多采用多边形模拟。其中，对于任意一个多边形骨料，其大小和形状都是随机确定的，其生成过程主要分为3部分：确定骨料的边数n、随机获得每个边对应的极角α和每个极角对应的极半径r。本文将多边形划分为多个三角形来计算骨料面积，通过在MATLAB中引入POLYAREA函数，将多边形骨料的顶点位置按顺序写入向量X和向量Y，可直接求得准确的骨料面积。
为提高混凝土多边形骨料模型的生成效率，引入射线法对骨料之间是否存在干涉进行判断，将2个面域是否干涉的问题简化为点是否处于面域内部的问题。射线法是一种判断点与面位置关系的高效方法，其基本原理为：对于给定的某点和多边形区域，以该点为端点向多边形区域做射线，保证产生1个或1个以上交点。若射线与多边形有奇数个交点，则表明该点在多边形内部，若有偶数个交点，则表明该点在多边形外部，如图1(a)所示。显然，若2个多边形区域相交，则必定有1个多边形的端点或其边线上的某些点在另一多边形内部。建立细观尺度混凝土随机骨料模型时，判断新生成骨料与已存在骨料是否干涉的具体步骤为：
1) 选取新生成多边形骨料的端点以及每边的10个等分点做控制点，判定每个控制点是否处于多边形骨料的外部。
2) 另选取已生成多边形骨料的端点做控制点，判定每个控制点是否处于新生成多边形骨料的外部。该步骤的目的是合理排除新生成多边形中每边相邻2个等分点恰好处在多边形骨料外部的极端情况，如图1(b)所示。当上述条件均满足时，可判断新骨料与已生成的骨料并未发生重叠，新骨料投放成功。
其他常见骨料干涉判断方法还有点阵法(图2(a))、面积标度法(图2(b)和2(c))和中心距法(图2(e))。与点阵法相比，射线法无需逐个判断骨料内部各个基点，极大缩短了判断流程。与面积法相比，射线法可有效克服新生成多边形骨料端点即使全部位于外部仍无法判断骨料是否重叠的缺陷，如图2(d)所示。与中心距法相比，射线法可实现骨料体积分数较高时的快速投放。基于射线法原理，建立混凝土骨料随机模型的主要步骤如图3所示。投放实例结果表明，射线法可将骨料投放体积分数从40%提升至70%，且投放时间控制在3 min以内，投放效率明显提升。
本文采用多边形骨料来模拟REN等[28]试验采用的石灰岩骨料，利用基于PYTHON的多边形骨料随机投放程序，生成长×宽为150 mm×150 mm的二维随机骨料模型。根据试验条件，粗骨料最小粒径为5 mm，最大粒径为15 mm，骨料的粒径分布符合二维平面的富勒级配分布曲线，可按照下式对混凝土粒径级配进行精确控制[26]：
式中：fi(D<Di)表示粒径小于Di的骨料体积占总混凝土体积的比例；fa表示骨料体积占混凝土体积的比例；Dmax为最大骨料直径，取15 mm。根据混凝土配合比反算可求得fa为45.1%[27]。按照[5，10) mm和[10，15] mm的骨料粒径范围划分，相应的骨料体积分数分别为28.6%和16.5%，得到的二维随机骨料多边形模型如图4所示。
1.2 有限元数值模型
1.2.1 内聚力单元的嵌入
利用ABAQUS对细观骨料断裂行为进行数值模拟。首先，对实体单元进行分割，并使用内聚力单元将其连接；然后，将内聚力单元具体分为砂浆内聚力单元和过渡区内聚力单元；最后，进行有限元断裂模拟分析。内聚力单元是在相邻实体单元之间插入新的零厚度单元，通过其自身的损伤和失效来模拟裂纹的发展情况。为了模拟混凝土断裂时裂纹发展的随机性，需要在有限元模型中所有实体单元之间插入内聚力单元并且赋予材料特性。ABAQUS仅能够实现特定路径(边界或层面)的内聚力单元插入，无法在任意实体单元之间插入内聚力单元，需要编制基于PYTHON的全局内聚力单元插入程序在所有实体单元之间插入内聚力单元。全局内聚力单元插入过程如图5所示，具体步骤如下。
1) 整理模型信息。获取骨料和砂浆实体单元信息，将单元及其所属的节点信息进行分类并建立二者映射关系。
2) 单元分离。根据节点与单元信息，计算每一个节点所连接的单元个数，并将节点分离为其对应的单元个数，将各个单元的公共节点替换为不重复的对应位置节点，实现单元分离。
3) 内聚力单元插入。根据节点与单元的映射关系，分析单元的相邻关系及公共边，根据得到的公共边节点信息，插入内聚力单元。
4) 模型文件生成。将内聚力单元和实体单元根据其位置进行分类，分为骨料单元、砂浆单元、骨料砂浆过渡区内聚力单元和砂浆内聚力单元，生成对应模型文件。
对二维随机骨料模型进行处理，可生成插入的二维模型内聚力单元，如图6所示。选用全局边长为1 mm的三角形网格对单元进行划分，得到28 788个砂浆单元、23 112个骨料单元、46 817个内聚力单元。对试块进行拟静态加载，上边界节点设置0.09 mm的竖向拉伸位移，加载时间为60 s，应变加载速率为1×10-5 s-1，对下边界全部节点施加竖向约束，仅最左边节点施加横向约束，如图7所示。
1.2.2 内聚力单元本构模型
利用有限元模型进行细观混凝土数值模拟，需要选用恰当的内聚力单元牵引力分离本构准则，目前主要有双线型准则、指数形准则和梯形准则。水泥砂浆的断裂通常采用双线形牵引力分离定律进行数值模拟，如图8所示。随着加载进行，界面单元的应力不断增大，待达到最大弹性应力后，界面单元的刚度按照裂纹扩展准则不断减弱，相应的承载能力也逐渐降低，直到界面单元的刚度降为零，此时，界面单元失去承载力并从模型中删除该单元，形成裂纹。
内聚力单元损伤起始准则主要有最大名义应力准则(MAXS)、最大名义应变准则(MAXE)、二次名义应力准则(QUADS)以及二次名义应变准 则(QUADE)。本文采用二次名义应力准则作为内聚力单元从损伤起始阶段达到损伤演化阶段的判断依据，即当内聚力单元名义应力满足式(2)时，内聚力单元出现损伤。
式中：
当满足二次名义应力准则时，内聚力单元就进入损伤演化阶段，该阶段的刚度kn按下式计算：
式中：
式中：
式中：un和us分别为拉剪混合模式下单元位移沿法向和切向的分量。
损伤演化分析主要有基于能量分析和基于位移分布2种。基于位移分析损伤演化时需要指定最大破坏位移
式中：
2 试验验证
为便于验证本文模拟计算的可靠性，拟采用与REN等[28]确定的相同试验参数。采用ASTM I型水泥作为砂浆，细骨料为天然硅质砂，粗骨料为破碎的石灰岩，最大粒径为15 mm，试件长×宽×高为150 mm×150 mm×50 mm，多组试验测得的平均混凝土抗拉极限强度为3.24 MPa。建模时，考虑到数值模拟的效率和模型的最小骨料粒径，初步设定全局网格尺寸为1 mm。在对内聚力单元进行划分之后，细观骨料数值模型的组成主要有砂浆、骨料、界面过渡区(ITZ)以及砂浆内界面(MII)，对内聚力参数进行合理调校，确定各细观组分的材料特性见表1。
| 参数 | 骨料 | 砂浆 | 边界内聚力单元 | 砂浆内聚力单元 |
|---|---|---|---|---|
| 弹性模量/GPa | 72 | 28 | 24 | 26 |
| 平面内泊松比 | 0.16 | 0.2 | ||
最大容许拉应力 | 2.3 | 4 | ||
最大容许切应力 | 10 | 30 | ||
Ⅰ型断裂能 | 0.01 | 0.1 | ||
Ⅱ型断裂能 | 0.25 | 2.5 | ||
| BK准则材料参数η | 1.2 | 1.2 |
基于蒙特卡罗理论，骨料的随机性影响混凝土结构力学特性，为了保证数值模拟结果准确，需对随机模型进行足够多的数值模拟。为此，先生成20组随机骨料模型，分别进行单向拉伸数值仿真，得到混凝土平均抗拉强度与模拟次数之间 的关系曲线如图9所示。从图9可以看出：1) 随着模拟次数增加，数值仿真得到的混凝土平均抗拉强度逐步收敛于试验测得的平均抗拉强度3.24 MPa，初步证明了数值模拟方法的可靠性；2) 当模拟次数达到10次时，数值仿真结果波动降低，并逐渐收敛。
3 混凝土断裂的随机性
3.1 细观断裂特征
为了分析混凝土受拉破坏全过程机理，研究细观断裂特征，随机选取3组混凝土样本，针对其中3个加载阶段，竖向最大变形分别为9×10-4 mm(第1阶段)、2.3×10-2 mm(第2阶段)和9×10-2 mm(第3阶段)，绘制混凝土试件内部不同位置的变形云图，如图10所示。图10中，颜色表示变形大小，试件最下端为蓝色，表示竖向变形为0，为竖向约束位置；试件最上端为红色，表示竖向变形达到最大值。由于图10中还绘制了粗骨料的位置，故试件上端未完全显示为红色。混凝土试件的变形云图可表征混凝土受拉过程的内部结构破坏情况。从图10可以看出：1) 在拉伸开始时，混凝土应变较小，内部结构没有发生损伤，混凝土内部变形分布较均匀；2) 随着应变增加，混凝土内部(主要是ITZ区域)出现变形断层，微小裂纹开始出现，此时，应力集中在微小裂纹的尖端，进一步促使裂纹扩大，导致应力的增长速度变小；3) 随着试件变形进一步增加，裂纹进一步扩展，混凝土内部变形出现明显断层，形成垂直于拉伸方向的贯通大裂纹，结构出现明显破坏，同时，应力随着应变增加逐步减小直到为0。由于骨料砂浆接触区抗拉强度较低，该区域的内聚力单元首先发生破坏，随后，裂纹两端出现明显的应力集中，导致裂纹进一步扩展直到贯通。该现象同时表明骨料的随机特性在细观层次上会影响混凝土裂纹的萌生和发展，进而导致混凝土断裂过程出现变异性。
3.2 细观受拉本构的变异性
由试验验证结果可知，为保证数值模拟的计算效率和计算精度，可仅选取其中10组混凝土随机骨料结构模型进行受拉全过程数值仿真计算，得到对应的混凝土应力-应变关系曲线中的变化带，并将单轴拉伸试验的数值模拟均值曲线与试验均值曲线相对比，如图11所示。从图11可以看出：1) 试验应力-应变均值曲线与数值模拟均值曲线较吻合，二者的峰值应力几乎相同；2) 各数值模拟应力-应变曲线在弹性阶段非常接近，而在达到极限强度后开始离散，随着应变进一步增大，离散性越显著。这是由于在初始变形阶段，混凝土结构作为一个统一的整体呈现弹性变化特征；随着变形增加，骨料形状与空间分布的随机性往往导致混凝土内部变形不均匀，混凝土内部出现明显的应力集中现象，导致部分内聚力单元被拉坏，此时，混凝土应力与应变之间的关系不再是确定的线性关系，数值模拟的抗拉应力峰值呈现离散状态；随着破坏的内聚力单元逐渐增加，裂纹发展也受到骨料内部不均匀应力的影响，数值模拟的裂纹发展的离散性更显著，应力-应变曲线的下降段出现明显的变异性。因此，仅通过1条确定性混凝土拉伸应力-应变曲线来描述混凝土的拉伸本构行为具有一定的局限性。
4 抗拉强度的变异性和模型化
4.1 加载速率对抗拉本构模型的影响
除了混凝土结构参数对其本构模型产生影响外，混凝土的力学响应也与加载速率密切相关。为模拟不同加载速率下混凝土拉伸应力-应变曲线、抗拉强度和断裂行为，深入研究加载速率对混凝土拉伸断裂的影响，选取10组骨料体积分数为45%的多边形骨料混凝土模型进行加载模拟，应变加载速率(简称应变率)分别取10-5、10-4和10-3 s-1，ABAQUS求解设置为隐式拟动态，得到不同加载速率下混凝土拉伸应力-应变曲线(见图12)和应变为2.4×10-4时混凝土的裂纹发展情况(见图13)。为清晰显示裂纹的发展，将试件变形放大系数统一取为100。从图12和图13可见：1) 据应力-应变关系曲线，随着应变加载速率增加，混凝土抗拉强度显著增高，当应变率由10-5 s-1提高到10-3 s-1时，其抗拉强度由3.22 MPa提高到7.52 MPa。 2) 从裂纹发展看，在相同应变下，随着应变加载速率的增加，混凝土裂缝长度与宽度愈小，对应的混凝土应力越大。这是由于高应变率下的混凝土裂缝仍集中在混凝土骨料周围，砂浆区域仍未贯通，能保持较高应力，混凝土结构难以发生内部损伤，混凝土抗拉强度较高。在低应变率下，混凝土裂缝已经在砂浆区域相互连接，基本形成1条贯通的大裂纹，混凝土应力开始明显下降，混凝土抗拉强度也相应降低。
4.2 骨料体积分数对混凝土本构模型的影响
混凝土拉伸数值模拟表明混凝土试件受拉性能主要受模型中骨料体积分数的影响。为研究骨料体积分数对混凝土抗拉性能的影响，首先，通过控制多边形骨料顶点与边数来控制骨料的畸变程度并生成随机骨料模型，取骨料体积分数为15%、25%、35%和45%的混凝土试件进行单向拉伸数值模拟。为减少骨料粒径对计算结果的影响，仅取骨料粒径范围为5~10 mm。生成4组不同骨料体积分数的混凝土模型并对其本构关系进行统计，得到不同骨料体积分数下混凝土的应力-应变曲线和抗拉强度与骨料体积分数的关系曲线，如图14所示。由图14(a)可知：随着骨料体积分数增加，混凝土结构初始弹性模量逐渐变大，但抗拉强度反而逐渐变小。其主要原因可能是随着骨料体积分数提高，骨料与砂浆的接触界面即ITZ区域也逐步增多，该区域抗拉强度较低，更加容易发生损伤并相互贯通，进而破坏混凝土结构。同时，由于砂浆相对于骨料的弹性模量较小，随着骨料体积分数增加，在相同应变下，混凝土弹性阶段的应力更大，最终导致混凝土结构弹性模量增加。图14(b)所示为混凝土抗拉强度均值ft,r与骨料体积分数fa之间的拟合曲线，两者呈现较好线性关系，其表达式为
4.3 细观概率本构模型
由于加载速率、骨料体积分数以及骨料的形状和空间分布对混凝土抗拉强度和本构关系影响显著，故混凝土强度的数值模拟结果呈现离散性。由图11可知：骨料形状和空间分布的随机性对于混凝土本构关系的影响主要体现在达到峰值拉应力之后，由于骨料形状与空间分布的随机性导致裂纹发展的随机性，最终导致本构关系出现高度变异性，因此，仅通过宏观混凝土本构模型难以解释和描述本构变异性现象。为此，选取70个体积分数相同(fa=45%)但空间分布随机的细观骨料模型，对其应力-应变曲线的峰值应力进行统计分析，以建立混凝土的细观概率本构模型。根据 GB 50010—2010《混凝土结构设计规范》[29]，宏观混凝土的拉伸本构关系可表示为
式中：
在混凝土拉伸本构关系中，αt与εt, r均可通过ft, r确定，且αt控制了拉伸本构曲线的下降段，因此，ft, r是该本构模型的关键参数。对70个随机混凝土模型的抗拉强度进行统计分析，得到ft, r的统计分布直方图及其拟合曲线，如图15所示。正态分布和对数正态分布均能较好地描述峰值拉应力的分布规律，正态分布略好，其相应的概率密度函数为
通过ft,r的统计分布可得到αt与εt, r的统计分布，从而将确定性的宏观混凝土抗拉本构曲线的参数概率化，得到概率型的混凝土抗拉本构模型。上述统计是根据骨料体积分数为45%的细观模型得到的，当骨料体积分数发生变化时，抗拉强度均值和标准差将发生变化。从图15可知：抗拉强度的标准差较小(σ<0.1 MPa)，可认为骨料体积分数变化时标准差不变，而抗拉强度的均值可由式(9)确定。根据表1，水泥砂浆(即不考虑粗骨料)的抗拉强度为4 MPa，将图14(b)中的抗拉强度量纲一化，得到量纲一化的抗拉强度
式中：
根据《混凝土结构设计规范》[29]，混凝土抗拉强度的均值
式中：α1为混凝土实体强度与立方体试件强度的差异影响系数，α1=0.88；α2为混凝土脆性影响系数，对C40混凝土，取1.00，对C80混凝土，取0.87；fcu,k为混凝土立方体抗压强度。
综合式(12)和式(13)，可得考虑骨料体积分数影响的混凝土抗拉强度均值为
因此，可得任意体积分数的混凝土的概率型峰值拉应力为
5 结论
1) 基于MATLAB，引入射线法判断骨料是否发生干涉，改进了混凝土随机多边形骨料生成算法，将骨料投放体积分数从40%提升至70%，投放效率明显提升。通过编制PYTHON程序，对ABAQUS软件进行二次开发，突破内聚力单元仅能在特定路径进行插入的局限性，实现了所有实体单元之间的内聚力单元插入，方便实现混凝土非均质特性的细观数值模拟。
2) 混凝土结构的数值模拟结果与试验结果较吻合，验证了该数值方法的可行性。随着模拟次数增加，数值仿真得到的混凝土平均抗拉强度逐步收敛于试验测得的平均抗拉强度3.24 MPa，验证了数值模拟方法的可靠性。当模拟次数达到10次时，数值仿真结果能达到相对收敛状态。
3) 细观混凝土结构的断裂行为和本构曲线呈现明显变异性。在初始变形阶段，混凝土结构作为一个统一的整体呈现弹性变化特征。随着变形增加，骨料随机性往往导致混凝土内部变形不均匀，混凝土内部出现明显的应力集中现象，导致部分内聚力单元被拉坏，此时，混凝土应力与应变之间的关系不再是线性关系，数值模拟的抗拉应力峰值呈现离散状态。随着破坏的内聚力单元逐渐增加，裂纹发展也受到骨料内部不均匀应力的影响，数值模拟的裂纹发展的离散性越显著，应力-应变曲线的下降段出现明显的变异性。
4) 基于70个骨料体积分数相同(均为45%)但空间分布随机的细观骨料模型进行统计分析，建立了混凝土的细观概率力学本构模型，较好地反映了骨料随机性对混凝土力学特性的影响。当应变率从10-5 s-1提高到10-3 s-1时，混凝土抗拉强度从3.22 MPa提高到7.52 MPa；随着骨料体积分数增加，混凝土结构初始弹性模量逐渐变大，但抗拉强度反而逐渐变小。
郭文华, 冶桐杰, 陈定市. 基于内聚力模型的混凝土受拉断裂数值模拟与概率本构模型[J]. 中南大学学报(自然科学版), 2025, 56(2): 586-597.
GUO Wenhua, YE Tongjie, CHEN Dingshi. Cohesive model-based numerical simulation and probabilistic constitutive modeling of concrete tensile fracture[J]. Journal of Central South University(Science and Technology), 2025, 56(2): 586-597.