2023级本科生“混凝土结构理论”课程研究报告系列
前言 上学期，我完成了《混凝土结构理论》课程的学习，对其中的混凝土单向受力状态比较感兴趣，因而选择了相关的题目：基于ABAQUS的内聚力单元分析混凝土的单轴压缩破坏，并撰写报告如下。
摘要：为探究混凝土单轴压缩破坏的细观机理，本研究基于ABAQUS软件，建立了包含随机分布骨料、砂浆基质及界面过渡区的二维细观模型。通过Python脚本实现骨料随机投放，并利用内聚力单元模拟界面过渡区行为。结果表明，破坏始于界面过渡区，并沿其扩展形成宏观裂缝，应力－应变曲线呈现典型三阶段特征。
关键词： ABAQUS；混凝土细观模型；Cohesive单元；破坏过程
混凝土作为最主要的工程材料之一，其破坏过程涉及复杂的细观损伤累积与裂缝扩展。宏观本构模型难以揭示内部损伤演化细节，因此细观数值模拟成为研究混凝土破坏机理的重要手段。本文基于ABAQUS软件，采用内聚力单元模拟界面过渡区，建立随机骨料细观模型，旨在直观呈现混凝土在单轴压缩下的损伤演化过程。
1 内聚力模型 基于连续损伤力学的内聚力模型采用在实体单元之间嵌入内聚力单元的方法来模拟损伤及断裂行为。内聚力单元是ABAQUS中用于模拟界面分离与裂缝扩展的一种特殊单元类型，该单元采用牵引力分离准则定义单元的损伤演化，并可通过定义断裂能的方式来模拟单元的断裂，而目前已有较多对混凝土I型和Ⅱ型断裂能的研究，并提出了相应的标准测试方法[1]。本研究着重讨论如何模拟界面过渡区（ITZ）的薄弱效应，采用零厚度内聚力单元表示界面过渡区（ITZ）。2 建模过程详解2.1 模型规划 本模型确定为尺寸100mm×100mm的二维平面模型。选择该尺寸，既符合混凝土常见试件的规格，便于与文献结果进行对比，又能将单元数量控制在个人计算机可承受的范围内。2.2 随机骨料生成的代码应用 在ABAQUS中，.rpy文件用于记录用户在CAE界面上的操作命令。这些文件可以帮助用户在进行二次开发时，提取和复用这些命令，从而提高编程效率。从网络开源分享，可找到用于生成随机骨料的Python脚本。其核心思路是：首先定义模型边界和骨料重叠判断函数，然后利用循环和判断语句实现骨料的不重叠随机投放。出于节约计算量的考虑，用脚本生成的骨料为圆形骨料。 混凝土骨料随机位置的确定需满足各骨料不相互接触，且有一定的最小距离，为此，假定每一个骨料颗粒都有一个影响范围，在该范围内其他骨料不能入侵，如图1所示。
图1 圆形骨料间的最小距离
在脚本中该关系用下式表述：
(2-1)
其中k取1.1，保证骨料间存在一定空隙。 模拟混凝土二维模型时需求得二维平面的骨料级配曲线，为保证混凝土骨料在二维级配和三维级配下的等效性，Walraven等建立了混凝土试件空间内骨料级配含量与其内截面内骨料面积之间的关系，即Walraven公式，该公式将三维Fuller级配曲线转化为混凝土试件内二维截面具有的内截面出现的概率[2]：
(2-2)
式中，为骨料占混凝土总体积的百分比，为级配骨料粒径，为不同骨料级配在试件二维截面上出现的概率。 在脚本中通过限制不同尺寸圆形骨料的颗粒数实现对骨料级配的模拟。 在ABAQUS中，通过“File -> Run Script”菜单可载入并执行脚本，可建立如图2所示二维构件。
图2 混凝土二维模型（构件）
2.3 内聚力单元的插入 若要在所有骨料－砂浆边界手动创建内聚力单元，不仅工作量大，且极易出错，故本研究使用相关插件实现在指定位置快速插入内聚力单元。 该插件以自定义工具栏的形式加载。在界面中指定要插入内聚力单元的部件，设置关键参数，即可将零厚度内聚力单元插入骨料与砂浆的接触面（如图3）。
图3 利用插件插入内聚力单元
3 材料参数的计算与输入 在混凝土细观有限元分析中，将宏观材料属性准确“翻译”并输入至分析软件，是决定模型物理真实性与计算可靠性的最核心环节，亦是本研究面临的主要难点。3.1 粗骨料：相对简单的弹性参数 在混凝土细观力学研究中，为平衡计算效率与模型精度，常将粗骨料简化为各向同性的线弹性材料。此简化基于一个基本假设：在单调荷载下，细观尺度的非线性破坏主要集中于界面区域，骨料的塑性变形相对次要。本研究参数主要参考陈涛在细观混凝土建模方面的研究成果[3]。具体取值如下： 粗骨料的弹性模量取为 50 GPa，平面内泊松比取为 0.2。3.2砂浆及其本构模型 熊学玉，肖启晟等在研究中提出可将砂浆基体也视作各向同性的线弹性材料，而将砂浆单元间的接触面设置为零厚度的内聚力单元，并定义其破坏准则。但此做法计算时间长且不易收敛，故采用另一思路：用ABAQUS中的混凝土损伤模型（CDPM）模拟砂浆材料[3]。 填写砂浆对应的CDPM参数应先计算砂浆材料在单轴压缩与单轴拉伸下的应力－应变关系。 本文取，由GB50010-2010规范[4]，单轴压缩下的应力-应变曲线可表示为：
(3-1)
式中：和分别为压应力和压应变； 是峰值应力对应的临界应变；为无损伤时的初始弹性模量；和 n是待确定常数。常数 可按下式估计
(3-2)
式中: 的单位为 MPa。 根据曲线的连续性条件，常数 n 满足如下方程
(3-3)
式中： 。一般地，砂浆的临界应变要比混凝土的大得多。本文中采用 Kohees等[5]的经验公式来估计砂浆和 ITZ 的临界应变，其表达式为：
(3-4)
对于单轴拉伸，应力-位移曲线由断裂能控制，并采用如下表达式[6]:
(3-5)
式中：和分别是拉应力和裂纹张开位移；是峰值拉应力，即抗拉强度；是最终失效时的位移，并有 ；和是两个材料常数，可分别取为3和6.93。 以上便是砂浆的受压与受拉本构。ABAQUS要求输入屈服应力－非弹性应变关系、损伤因子－非弹性应变关系、屈服应力－开裂应变关系及损伤因子－开裂应变关系，已知砂浆单轴压缩与单轴拉伸的应力应变关系，通过以下公式便可计算出所需数据：
(3-6)
式中，为受压非弹性应变；和为压应力和压应变；为初始弹性模量对应的弹性压应变；
(3-7)
式中，为受拉开裂应变；和为拉应力和拉应变；为初始弹性模量对应的弹性拉应变。
关于损伤参数，一般采用Sidoroff基于提出的损伤因子计算方法：
受压：
受拉：
（3-8）
借助Excel计算出所需四组数据，并填入ABAQUS中。最后还需填入如图4混凝土损伤塑性参数，其中粘性参数取0.0005可同时满足精度和收敛性的要求。
图4
3.3 界面过渡区（ITZ） 采用基于牵引－分离本构的零厚度内聚力单元对ITZ进行模拟，其力学行为由弹性阶段、损伤起始与损伤演化三部分定义。
1. 弹性刚度参数的确定 内聚力单元的初始线弹性行为由法向刚度和切向刚度描述（对应ABAQUS中的、与）。这些刚度并非直接测量获得，而是基于其所表征的薄弱层（ITZ）的等效弹性性质计算得出，其计算公式为：
(3-9)
从网络资源及ABAQUS的帮助文档了解到，是一个虚拟的“有效厚度”，用于将三维材料常数转换为面内的刚度密度。通常取1.0 mm以保证数值稳定性。剪切模量由以下公式计算得到，其中泊松比参考砂浆取值
(3-10)
ITZ的弹性模量及后续各项参数参考陈涛研究中的取值[3]。
2. 损伤起始与演化参数的计算与输入 损伤起始判断采用二次名义应力准则（Quadratic Nominal Stress Criterion，在ABAQUS中简称Quads）。当名义应力状态满足公式（3-11）时，损伤开始。其中、、 分别为当前法向与两个切向应力；、、为对应的峰值强度。
(3-11)
损伤演化规律采用基于能量的Benzeggagh-Kenane（BK）混合模式准则。该准则定义了材料从损伤起始到完全失效所消耗的断裂能，其核心参数包括：I型断裂能 ：对应于法向张开破坏。Ⅱ型断裂能 ：对应于第一剪切方向破坏。BK准则材料系数 η：控制混合模式下等效断裂能 在ABAQUS中创建名为“ITZ”的材料并在相应位置输入上述参数，定义其弹性行为、Quads损伤及损伤演化，完成对界面过渡区材料属性的定义。4 分析设置、求解与后处理4.1 分析设置 在ABAQUS中进行隐式静力学计算。设置加载板与试件的接触关系（如图5）与边界条件，荷载以加载板位移的方式施加：上加载板向下位移1.2 mm，下加载板不移动。 此外，对内聚力单元启用单元删除以模拟裂缝形成。
图5 加载板与试件的接触
4.2 求解监控与后处理 在后处理模块中，可查看如下结果：S（应力）、U（位移）、DAMAGEC（砂浆受压损伤）、DAMAGET（砂浆拉伸损伤），结果如图6－图9的云图所示。
图6 S(应力)云图
图7 U(位移)云图
图8 DAMAGEC(砂浆受压损伤)云图
图9 DAMAGET(砂浆拉伸损伤)云图
在“创建XY数据－-ODB历程变量输出”中可分别输出上加载板的反力与位移随时间变化的曲线。 在“创建XY数据－操作XY数据”中用combine（X，X）运算符可换算得到如图10所示的应力应变曲线。
图10 混凝土试件的应力-应变曲线
5 结果观察、分析5.1 细观破坏过程分析 对损伤云图及内聚力单元失效过程的观察表明，混凝土的破坏始于细观层面。初始损伤集中出现在界面过渡区，这与ITZ因其较高孔隙率与脆弱性而被视为混凝土中最薄弱相的理论认知一致。随后，微裂缝主要沿ITZ扩展并相互连接，形成复杂的裂缝网络。由于模型考虑了加载板与试件端部的摩擦约束，裂缝的最终贯通形成了一条明显的宏观倾斜破坏带，试件整体呈现近似“锥形”的破坏形态。 此破坏模式直观地揭示了混凝土破坏是一个由细观损伤萌生、汇聚直至宏观局部化贯穿的渐进过程。5.2 宏观力学响应与模型评估 模拟所得的应力-应变曲线呈现出典型的三个阶段：近似线弹性阶段、非线性上升段以及峰值后的软化下降段。 然而，与典型试验曲线相比，本模型响应存在可辨差异：弹性阶段的斜率（表征弹性模量）偏低，曲线峰值位置偏后；峰值后曲线下降较为平缓，且在应力从峰值约55 MPa下降至30 MPa附近时出现一明显拐点。这些差异可能源于以下多方面因素的共同影响：1）细观参数的不确定性：ITZ界面属性及砂浆本构参数的取值对模拟响应极为敏感，文献参考参数可能不完全适用于本模型，这直接影响初始刚度与软化行为。2）模型简化效应：二维平面应变模型无法反映真实三维约束状态，可能高估了横向约束，从而影响破坏模式与软化段形态。峰值后拐点的出现，可能反映了内部骨料骨架在主要裂缝形成后仍具有一定的残余承载机制。5.3 学习总结与模型意义 本次模拟实践表明，基于随机骨料与Cohesive单元的细观模型能够有效地复现混凝土破坏的核心机理，即损伤始于ITZ并主导破坏过程。尽管模拟结果在定量上与实验存在差距，但此过程深刻揭示了输入参数的代表性与模型简化的合理性是决定数值模拟预测能力的关键。该模型为理解混凝土宏观非线性行为背后的细观物理机制提供了一个有力的可视化分析工具。
参考文献
[1] 熊学玉,肖启晟. 基于内聚力模型的混凝土细观拉压统一数值模拟方法[J].水利学报,2019,50(04):448-462.DOI:10.13243/j.cnki.slxb.20181061.[2] 陈青青,张煜航,张杰,等. 含孔隙混凝土二维细观建模方法研究[J].应用数学和力学,2020,41(02):182-194.[3] 陈涛.基于改进内聚力本构模型的混凝土细观损伤和断裂模拟研究[D].大连理工大学,2023.DOI:10.26991/d.cnki.gdllu.2023.003342.[4] 中华人民共和国住房和城乡建设部. 混凝土结构设计规范 GB50010-2010(2015年版). 北京：中国建筑工业出版社, 2015.[5] Kohees M., Sanjayan J., Rajeev P. Stress-strain relationship of cement mortar under triaxialcompression[J]. Construction and Building Materials, 2019, 220: 456-463.[6] Hordijk D. A. Tensile and tensile fatigue behaviour of concrete; experiments, modelling andanalyses[J]. Heron, 1992, 37(1).
参考网络资料：
https://www.bilibili.com/video/BV1x54y1T7ri/?spm_id_from=333.1387.favlist.content.click&vd_source=8a9a656319b72bfb1e543506789a2409
https://www.bilibili.com/read/cv8532259/?from=search&spm_id_from=333.337.0.0&opus_fallback=1
https://www.bilibili.com/video/BV1qK411n7At/?spm_id_from=333.1387.favlist.content.click&vd_source=8a9a656319b72bfb1e543506789a2409
https://zhuanlan.zhihu.com/p/592627971
