基于内聚力模型的混凝土细观拉压统一数值模拟方法

熊学玉，肖启晟

（1.同济大学 土木工程学院，上海 200092；2.先进土木工程材料教育部重点实验室（同济大学），上海 200092）

1 研究背景

混凝土的破坏主要是因为初始缺陷发展而导致的准脆性断裂［1］。为描述和分析混凝土的断裂行为，早在1960年代，Kaplan［2］就将经典断裂力学理论引入到了混凝土结构的分析中，并已在实际工程中证明了线性断裂力学在大体积混凝土结构中的适用性。但由于混凝土的微观裂缝和亚临界裂缝尺寸相较于一般的金属材料大几个量级，且裂缝尖端周围的应力场分布也呈现出一定的不均匀性，故弹塑性断裂力学的COD和J积分理论亦不适用于描述混凝土在细观尺度上的断裂［3］。因此普通大小的混凝土构件的断裂问题不太适合用经典断裂力学理论进行分析。

鉴于混凝土的宏观断裂行为是在较低尺度上损伤和断裂累计的表现［4-5］，学者们尝试在细观尺度上通过数值模拟来分析和研究混凝土的宏观断裂。近年来，细观尺度的数值模拟逐渐成为热点［6］，并已在静载［7-9］和动载［10-12］作用下的混凝土损伤机理方面取得了一些研究成果。

目前围绕混凝土的断裂行为，已提出了多种断裂仿真模型，其主要分为连续和非连续裂缝模型两类［12-13］。连续裂缝模型通过单元的刚度损伤来描述裂缝的发展，裂缝位置处模型保持几何连续，即不能在几何外观上直接表现裂缝的展开。经典连续裂缝模型包括弥散裂纹模型［14］、钝裂缝带模型［15］、旋转裂缝模型［16］等。非连续裂缝模型则允许定义节点的分离准则来模拟裂缝的发展，即能够模拟裂缝位置处的几何不连续。经典的非连续裂缝模型包括晶格模型［17］、刚体弹簧模型［18］、界面模型［18］等。虽然这些裂缝模型能给出直观的裂缝模拟效果，但计算结果的网格敏感性较大，收敛性较差。内聚力模型［20］和扩展有限元法［21］因能克服上述其他不连续裂缝模型算法上的缺陷，而更常见应用于混凝土裂缝的数值模拟中。其中，内聚力模型在模拟多裂缝展开问题上更具优势。

基于连续损伤力学的内聚力模型是采用在实体单元之间嵌入内聚力单元的方法来模拟损伤以及断裂行为［20］。内聚力单元采用牵引力分离准则定义单元的损伤演化，并可通过定义断裂能的方式来模拟单元的断裂。国内外已有较多对混凝土Ⅰ型和Ⅱ型断裂能的研究［22-28］，并提出了相应的标准测试方法［29-32］。相较于基于应力场和应力强度因子的断裂判定方法，内聚力模型在细观尺度上从能量的角度来模拟断裂，既节省计算成本，同时又能够契合前述混凝土断裂行为的特殊性。

由于混凝土的单轴受拉断裂破坏机理较明确，且易于采用牵引力分离准则来描述损伤演化，故目前内聚力模型已在混凝土受拉断裂模拟中得到了较多的应用［33-38］。而混凝土在压力作用下的断裂破坏机理相较更为复杂［39］，目前采用该模型来模拟混凝土的受压损伤断裂行为的报道较少。

本文在混凝土内聚力模型的基础上，进一步定义内聚力单元混合（Ⅰ+Ⅱ）型断裂准则以及全局的接触行为，将模型的适用范围从原来的混凝土细观受拉断裂模拟拓展至细观受压断裂模拟，并实现细观拉压数值模拟方法的统一。通过对混凝土在单轴拉压作用下的受力全过程模拟分析以及参数研究，探讨混凝土细观断裂损伤与宏观受力响应的关联，讨论断裂和接触行为对混凝土在拉压作用下断裂损伤行为的影响，以期对混凝土细观断裂模拟以及破坏机理研究提供参考。

2 混凝土拉压统一模型的建立

2.1 模型组成与网格划分模型分为4个组成部分：骨料、砂浆颗粒、界面过渡区（ITZ）、砂浆内界面（MII），如图1所示。其中，骨料以及砂浆颗粒均采用实体单元模拟，界面过渡区（ITZ）以及砂浆内界面（MII）采用零厚度的内聚力单元模拟。为生成更多的内聚力单元以减少模型的网格敏感性，模型在平面内划分为三角形网格，实体单元均为楔形单元。为建模方便以及节省计算成本，模型内所建骨料的最小粒径限制为2.5mm，更小粒径的骨料和水泥基体合并考虑为均质的砂浆单元（由砂浆颗粒和MII单元组成）。现暂假定骨料在受力过程中不发生损伤和断裂，保持线弹性受力状态。对于考虑骨料损伤和断裂的情况，将于6.4节中讨论。

图1 数值模型组成部分

模型中的骨料抽象化地用任意多边形来模拟。采用定义在极坐标系上的随机函数来生成任意多边形顶点［40］：

2.2 接触行为本文采用有限元软件ABAQUS 进行数值模拟。鉴于内聚力单元非线性的损伤演化，断裂引起的模型几何不连续以及全局单元之间的碰撞接触容易造成计算的不收敛，为提高模型的收敛性，采用显式求解模块进行计算。由于混凝土在受压开裂后会发生颗粒间的碰撞与错动，故采用全局接触来定义接触行为。由于全局接触算法只能应用于三维面接触［41］，因此虽然进行的是混凝土二维的单轴拉压模拟，但采用三维模型来建模分析。模型的法向接触行为定义为“硬接触”，即假定接触面之间法向仅能传递压力，且不限制压力大小。采用库伦摩擦模型定义切向接触行为，即假定摩擦系数为一定值，且摩擦力大小不受限。摩擦系数μ根据混凝土的内摩擦角φ来计算取定［42］：

2.3 材料模型

2.3.1 实体单元 由于采用三维模型进行二维混凝土单轴拉压模拟，因此需要消除模型平面内与平面外的变形耦联。此外，在高精度的网格划分下，实体单元的损伤由其周边的内聚力单元来描述，故不定义实体单元的损伤行为。赋予实体单元各向异性的线弹性材料模型，单元类型选择为六节点三角形棱柱体（C3D6）单元。材料出平面方向的泊松比设为0以消除出平面方向的泊松效应。

2.3.2 内聚力单元 砂浆与骨料之间的界面过渡区（ITZ）以及砂浆内界面（MII）均采用内聚力单元模拟，单元类型为八节点三维内聚力（COH3D8）单元。鉴于双线型模型适用于描述脆性材料的断裂［42］，且拥有较高的计算效率［44］，故采用双线型模型描述内聚力单元的牵引力分离关系。单元的初始损伤采用平方名义应力准则来判定：

式中：tn、ts和tt分别为在受拉方向和两个剪切方向上的名义应力；分别为在受拉方向和两个剪切方向上的最大允许名义应力。

由于混凝土在拉力和压力作用下分别会产生Ⅰ型和Ⅱ型断裂为主的破坏［45］，且混凝土的非均质性和多相性会导致混合型断裂，因此需定义内聚力单元的Ⅰ型、Ⅱ型以及混合型断裂的损伤演化。目前，采用不同类型的混凝土Ⅱ型断裂能测量方法所得的测量结果［28-29，31］，以及通过混凝土直剪试验结果推算混凝土的Ⅱ型断裂能［46］，均显示混凝土的Ⅱ型断裂能约为Ⅰ型断裂能的25倍。而目前对界面过渡区以及砂浆断裂的研究较少，故暂参考这一倍数关系来定义内聚力单元的断裂能。对于介于Ⅰ型和Ⅱ型断裂之间的混合型断裂，采用Benzeggagh-Kenane（BK）混合断裂准则［47］进行描述，混合型断裂能GC的计算表达式如下：

在采用平方名义应力准则判定初始损伤和采用BK混合断裂准则描述混合型断裂的定义下，内聚力单元的损伤演化本构关系如图2所示。由图2可见，通过定义单元的断裂能和混合断裂准则可确定单元断裂时的极限裂缝宽度［48］和剪切变形。

3 细观混凝土建模与校调

图2 内聚力单元损伤演化本构关系

本文选取Ren 等［48］所做的混凝土单轴拉压试验作为数值模拟的对象。混凝土的平均单轴受拉和受压弹性模量为42.0和37.5 GPa，平均极限抗拉强度和抗压强度为3.24和51.20 MPa。试件的平面几何尺寸为150 mm×150 mm。根据混凝土材料配比计算得骨料的体积占比为45.1%［37］。

由于在细观层次上对混凝土建模需定义内部材料组分以及材料边界之间的本构关系，而目前在混凝土内聚力模型基础上展开的材料试验研究较少。再者，数值模拟中的断裂路径与网格划分直接相关，过粗或过精细的网格划分可能会导致网格敏感性［50］或计算成本过大。故进行模型参数校调以保证的数值模拟方法的可行性与准确性。

为在保证计算结果精度的前提下提升计算效率，对同一模型划分了3种不同尺寸的网格进行对比计算。单元近似尺寸s分别为1.5、1.0和0.5 mm，网格精度的对比情况如图3所示，单元统计信息汇总于表1。3种不同网格精度模型的计算成本之比为1∶2.5∶11.6，大致与单元数量呈正比。其中，调用3.8 GHz主频CPU的16个线程计算1.5 mm网格精度的受拉模型耗时约6.8 h，计算受压模型耗时约8.9 h。

图3 不同单元近似尺寸下模型网格精度

表1 单元数统计

不同网格精度模型的拉压应力、应变关系如图4所示。从图4可以看出，相比于受拉数值模拟结果，混凝土受压数值模拟结果的网格敏感性较为明显。由于模型的损伤以及塑性变形均通过内聚力单元来模拟，因而随着内聚力单元在单位面积上数量的增加，局部区域拥有更均匀的塑性变形能力。在损伤单元的周边区域内产生内力重分布，试件的承载能力拥有一定的上升空间。由于混凝土受拉强度一般取决于最弱截面，且试件内单一截面的面积几乎不受内聚力单元数量的影响，故网格敏感性较不明显。而受压行为则取决于试件的整体响应［45］，内聚力单元数量过少会导致局部损伤集中，进而限制模型的内力重分布，故网格精度对混凝土的受压行为的影响更为显著。

图4 不同网格精度下的拉压应力-应变曲线

表2 模型材料参数

鉴于单元近似尺寸s为1 mm时，即单元尺寸与模型边长尺寸之比为1/150时，模型拥有较高的计算效率和精度，故下文的研究中采用这一网格精度来划分网格以及内嵌入内聚力单元。

在参数校调过程中发现，即使不同骨料分布的模型试件符合相同的骨料级配分布，但会因断裂路径的不同而导致模型在非线性受力阶段的力学响应以及破坏形态上产生明显差异。为取得更具普适性的研究结果，本文建立了一套包含10 种不同骨料分布的模型来进行参数校调以及后文研究分析。各模型的骨料的粒径分布均服从富勒级配曲线，其在二维平面的分布计算式如下［51］：

基于第1节中所述的材料试验结果以及模型单元特征，分别对模型中材料、界面以及接触关系进行了参数校调。经模型校调验证，摩擦角为55°时能取得较好的力学响应和断裂模拟效果，其他材料参数见表2。标准模型的ITZ单元和MII单元的纯Ⅰ型断裂的极限裂缝宽度为0.019和0.050 mm。

4 受拉数值模拟结果

受拉数值模拟的破坏形态如图5所示，应力、应变曲线如图6所示。为获得更具代表性的计算结果，图6中绘制了平均曲线来代表试件组的整体力学响应趋势。曲线上每一点代表10个不同骨料分布的模型试件在该全局平均应变下的应力平均值。图6中也绘制了标准差与离散系数曲线来表示模型组的绝对和相对误差程度。

由图6可见，标准差曲线的上升段主要集中在试件达到最大抗拉强度状态和其后的一小段变形范围内，这说明了因骨料分布而引起的受拉性能的差异主要体现在软化段初期，其后骨料分布对混凝土软化段剩余强度的退化速率影响较小。

图5 受拉破坏形态对比

选择模型组中一个试件进行变形全过程分析，试件的应力、应变曲线如图7所示，结合试件的受力变化过程，在曲线上选择5个特征点代表试件在变形过程中的5个特征状态。如图8所示，为方便观察模型在拉力作用下的损伤演化过程，对模型单元节点在加载方向上的位移按比例Dv进行放大，令试件的宽高比由原来的1∶1放大为1∶2。

图6 受拉数值模拟统计结果

图7 分析试件的受拉应力-应变曲线

图8 模型试件在各特征状态下放大的变形（红色为MII单元，浅蓝色为ITZ单元）

由图8可见，试件在弹性受力阶段中的变形均匀地分布在骨料和砂浆中。自状态A开始，试件逐渐开始出现塑性变形。尽管在该状态下混凝土的平均拉应力尚未达到过渡区界面（ITZ）单元的初始损伤应力，但由于砂浆和骨料颗粒之间的弹性模量差异所导致的应力集中，部分ITZ单元已开始出现损伤。在状态B下，混凝土的应变已逐渐集中至ITZ单元以及其临近的砂浆内界面（MII）单元，在骨料与砂浆之间形成纺锤状的变形集中区域。当试件达到最大的抗拉强度时，部分变形集中区域明显扩张，形成了两条潜在的断裂带（FZ1与FZ2）。在状态D下，潜在的断裂带FZ2逐渐发展成熟。在断裂带位置处，变形集中区域出现联通。其中一个ITZ单元发生完全损伤，即出现了宏观的断裂。MII单元的变形迅速变大，变形集中区域的纺锤状特征几乎消失。值得注意的是，在状态D之后，主断裂带的右下角出现了一块新的变形集中区域（FZ4），并且该区域于状态E时出现了宏观断裂。而该区域的上方是原来在状态D时即形成的变形集中区域（FZ3），该区域的变形随新区域的扩展逐渐萎缩。随着宏观裂缝的扩展，试件最后的宏观断裂路径大致会穿过受载过程中所形成的主变形集中区，但会因为内力重分布而导致其绕过原变形集中区域，形成并穿过新的断裂带。

5 受压数值模拟结果

试件的受压应力应变曲线如图9所示。与受拉数值模拟结果相同的是，在达到平均极限强度状态之前各试件的力学响应保持较低的离散性。而在强度软化阶段，虽然各模型的力学响应出现明显差异，但各模型的应力、应变曲线在剩余强度大致为极限强度的1/2至1/3时出现了“捏缩”的现象。由离散系数曲线上可见，在“捏缩”的位置处曲线出现了一个局部极小值的“谷”。这个“谷”表明混凝土在强度软化段存在一个相对稳定的受力状态。

图9 受压数值模拟统计结果

图10 分析试件的受压应力-应变曲线

图11 分析试件在各状态下的裂缝扩展

图10为模型组中一个试件的应力-应变曲线，在曲线上选择了A 至F 共6 个标志性状态进行分析，对应各状态下模型的宏观裂缝如图11所示。图12和图13给出了不同损伤程度下的ITZ和MII单元的损伤长度。

图12 界面过渡区（ITZ）单元损伤长度-应变曲线

图13 砂浆内界面（MII）单元损伤长度-应变曲线

损伤程度在ABAQUS内采用SDEG值描述，SDEG值表示初始刚度损伤率（损伤因子），当一个单元SDEG值等于1时表示该单元断裂。损伤长度为在对应损伤程度范围内的界面单元的长度之和，单元长度通过提取界面单元的节点坐标信息计算得到，其表达式如下：

式中：L( D )为损伤长度；D 为损伤因子，用SDEG值表示；n为界面单元的总个数；li为第i个界面单元的长度；Di为第i个界面单元的损伤因子。

在初始损伤发生前，试件保持线弹性。由图12 和图13 可见，初始损伤最早出现在ITZ 单元中，其后再出现在MII单元中。在损伤发生后，出现损伤的单元数量的增长速率逐渐增大并达到一个相对稳定值。与受拉数值模拟所不同的是，在试件达到极限承载能力之前，于状态A便有ITZ出现完全损伤，形成了第一条平行于力加载方向的宏观裂缝。在状态B与状态C之间的受力阶段，裂缝逐渐出现在试件底边中部与右上角连线的ITZ单元位置处，并形成了第一个剪切带的雏形。

在状态C与状态E之间，试件的强度迅速下降，砂浆内部出现裂缝，并导致界面过渡区的裂缝产生分支与联通。多数的宏观裂缝在状态E之前便已形成，而状态E位于应力应变曲线的“捏缩”位置。因此推断前述的相对稳定状态是试件的损伤演化至相对成熟后的一个状态。宏观裂缝形成多条剪切带，各剪切带首尾相连并将试件分割成多个楔形块，整个试件的破坏形态已初步成型。此后试件的变形逐渐转化成各混凝土楔形块之间的滑移。

6 参数分析

6.1 断裂能以往研究中一般认为混凝土在单轴拉力作用下的破坏主要是Ⅰ型断裂为主，因此在受拉数值模拟中仅采用一种断裂能描述内聚力单元的断裂［33-38］。虽然这种数值模拟方法在受拉数值模拟中能够得到较满意的模拟结果，但由于混凝土属于多相非均质材料，即使单轴受拉也无可避免地会发生混合型断裂。故仅采用一种断裂能描述内聚力单元的断裂，并以此来模拟混凝土的断裂行为是有待商榷的。此外，混凝土在压力作用下所产生的破坏很大一部分是由剪切作用所造成的，破坏以Ⅱ型断裂为主［45］，故采用一种断裂能无法同时满足混凝土受拉和受压的数值模拟。由于Ⅰ、Ⅱ型断裂能决定混合型断裂能的大小、单元断裂时极限裂缝宽度和剪切变形以及损伤演化行为，因此Ⅰ、Ⅱ型断裂能对混凝土的受压和受拉行为会无可避免地产生影响，故有必要围绕断裂能问题展开参数分析研究。

本文在以表2 为材料参数的标准模型的基础上，分别建立了Ⅰ型和Ⅱ型断裂能为标准模型的60%、80%、120%和140%的对比模型组，其中Ⅱ型断裂能的参数分析中还考虑了Ⅱ型断裂能等于Ⅰ型断裂能的情况计算结果如图14和图15所示。

图14 Ⅰ型断裂能对应力-应变的影响

图15 Ⅱ型断裂能对应力-应变的影响

由图14可见，Ⅰ型断裂能主要影响混凝土在拉、压作用下强度软化段的受力性能，略微影响混凝土的极限强度。由于混凝土在单轴拉伸下主要发生的是Ⅰ型断裂，因此混凝土在拉力作用下的剩余抗拉强度对于I型断裂能的变化更为敏感。

由图15（a）可见，Ⅱ型断裂能与Ⅰ型断裂能在对混凝土的受拉性能的影响上相似。当Ⅱ型断裂能等于Ⅰ型断裂能时，相比于标准模型，整个试件的能量消耗下降了48.2%。这说明混凝土在单轴拉力作用下虽然发生的是偏向于Ⅰ型断裂的混合型断裂，但其整体受拉破坏的能量消耗远大于细观尺度的Ⅰ型材料断裂能。因此在细观尺度下，Ⅱ型断裂能对混凝土的宏观受拉断裂有不可忽略的影响。

由图15（b）可见，混凝土在非线性受力阶段的受力性能会明显地受Ⅱ型断裂能的影响。正如第4节中对受压全过程所分析的，在试件达到极限抗压强度之前即会产生宏观的剪切裂缝，即抗压强度受限于试件的变形能力。故较小的Ⅱ型断裂能会降低试件的变形能力，限制试件在出现局部损伤后的内力重分布，进而降低极限抗压强度。从各试件的应力应变曲线上升段的对比可见，随着Ⅱ型断裂能的增加，混凝土极限抗压强度的上升逐渐放缓，并且达到极限强度时所需的塑性变形增多。这说明了，混凝土试件的极限抗压强度逐渐从变形控制转化为强度控制。

6.2 界面过渡区与砂浆内界面强度从5.1节的分析中可知，混凝土的极限拉压强度不仅受断裂能控制，还受单元的强度控制。此外，由于混凝土的非均质和多向性所导致的混合型断裂同时受剪切与拉伸作用控制。因此本节在标准模型的基础上，于20%～200%的参数取值范围内，分别研究了界面单元抗拉和抗剪强度对混凝土极限抗拉和抗压强度的影响（见图16）。

由图16（a）可见，在界面单元抗拉强度较小时，混凝土的极限抗拉强度与内聚力单元的抗拉强度近乎呈线性关系。由于模型的断裂能不变，因此在较小强度下，内聚力单元拥有较强的变形能力，混凝土的抗拉强度主要受单元强度控制。随着单元强度的增大，单元的变形能力下降，混凝土在达到极限强度前的塑性变形空间减小，因而强度增长空间也会逐渐受限。一般情况下，单元的剪切强度几乎不影响混凝土的抗拉强度。当剪切强度为标准模型的20%时，混凝土的抗拉强度下降仅为7.4%。

由图16（b）可见，单元的抗拉和抗剪强度的上升均能提升混凝土的抗压强度，但抗拉和抗剪强度与混凝土抗压强度之间呈不同形式的非线性关系。在较小强度范围内，单元的抗拉强度的提升能够更明显提升混凝土的抗压强度，而在较大的强度范围内，混凝土的抗压强度更敏感于单元的剪切强度的变化。

图17提取了20%与200%单元抗拉和抗剪强度下的破坏形状，并通过与标准模型的破坏形态（图11（f））对比，可见剪切强度较小的情况下混凝土内部几乎没有形成贯通全截面的断裂带。而随着剪切强度的上升，混凝土内部会形成更多贯通全截面的剪切带，且局部裂缝发展也会越加充分。相比之下，单元抗拉强度的提升只会略微增加裂缝数量，几乎不影响内部断裂带的走向。这也说明了混凝土的受压破坏在细观层次上主要表现为剪切破坏。

图16 界面单元强度对混凝土抗拉和抗压强度的影响

图17 不同界面强度模型的破坏形态

6.3 骨料形状由于圆形骨料建模相对便捷，因此较常应用于细观混凝土的建模中［22，52］。相比于不规则形状骨料，圆形骨料会减少混凝土内部的应力集中程度，进而导致混凝土强度相对偏大［53］。本文建立了圆形、正六边形、正五边形和正方形骨料模型试件组，进行骨料形状的参数研究，每个不同形状骨料模型组包括4种不同骨料分布的模型。

各模型组的平均应力应变曲线如图18所示。从图18可见，随着正多边形骨料边数的增加，混凝土的极限抗拉和抗压强度均有不同程度的上升。其中，正五边形骨料模型的极限抗拉和抗压强度与任意多边形骨料模型相当。

在拉力作用下，骨料边数越多，模型达到极限抗拉强度时所对应的平均应变也越大。相比于任意多边形骨料模型，圆形骨料和正方形骨料模型的抗拉强度分别提高和降低了3.1%和2.5%。产生这一规律可能有以下两点原因：（1）角缘效应的影响［54］，多边形骨料周边ITZ单元的面积大于同体积圆形骨料周边ITZ单元的面积（即多边形骨料周边的薄弱面更多），试件在相同的应力水平下，会有更多的多边形骨料周边ITZ单元会出现损伤和断裂；（2）应力集中的影响［53］，多边形骨料周边的应力分布相对更加集中，多变形骨料周边的ITZ单元更易出现损伤和断裂。

在压力作用下，模型的极限抗压强度也随着骨料边数的增加而上升，其同样可能受到角缘效应和应力集中的影响。相比于任意多边形骨料模型，圆形骨料和正方形骨料模型的抗压强度分别提高和降低了9.3%和4.4%。此外，虽然骨料边数越多混凝土抗压强度越高，但同时其软化段的剩余强度的退化也越快。

图19为各骨料形状模型的破坏形态。从图19可见，随着骨料边数增加，剪切带的倾斜角度减小，剪切带分布的无序性降低，被混凝土分割的楔形块大小更均匀。这说明骨料越光滑则模型整体的受力和出现的损伤会更趋向于平均分布，但出现损伤后的骨料咬合效应越弱，进而导致强度退化更快。

图18 骨料形状对应力-应变的影响

图19 不同骨料形状模型的破坏形态

表3 骨料内界面单元（AII）材料参数

6.4 骨料的材料类型上文的数值模拟假定混凝土受力过程中骨料保持线弹性，虽然一般认为混凝土的受拉性能由ITZ和砂浆控制，较少出现骨料拉断的现象［55］，但混凝土在压力作用下会出现骨料的损伤和断裂［56］。为研究骨料的损伤和断裂对混凝土受压性能的影响，本节在原标准模型组的基础上,于骨料的网格间也嵌入了内聚力单元来模拟骨料的内部界面AII，并分别赋予了石灰岩［57］、大理岩［58］、花岗岩［59］和玄武岩［60］共4种不同类型骨料的材料模型，对应AII单元的材料参数如表3所示。骨料的网格间嵌入内聚力单元后，相比于原模型，内聚力单元数量增加了约70%，计算耗时增加了约115%。

各模型组的受压应力-应变曲线如图20所示。相比于采用线弹性骨料简化的标准模型，采用石灰岩、大理石、花岗岩和玄武岩骨料的混凝土模型最大抗压强度分别下降了1.1%、1.0%、0.23%、0.18%。这一点说明模型组极限抗压强度的主要控制因素仍是ITZ的强度和变形能力，而骨料材料类型的影响则主要表现在强度软化段的受力性能上。其中，采用大理岩骨料的混凝土模型在软化段的剩余强度更高，强度退化更慢。相比与其他骨料，大理岩骨料的AII 单元的强度较低，略大于砂浆MII单元强度，因此更容易发生损伤和断裂。图21为提取的一组骨料内界面单元的损伤长度-应变关系曲线。从图21可见，大理岩骨料的损伤最为明显，其在受力过程中的塑性变形有利于局部区域的内力重分布，因此软化段的强度退化较缓。

图20 骨料的材料类型对混凝土受压应力-应变曲线的影响

图21 骨料内界面单元损伤长度-应变曲线（SDEG＞0.95）

图22 不同骨料材料类型模型的破坏形态

图22为一组相同骨料分布模型的破坏形态。从图22可见，玄武岩骨料模型（图22（d））的整体破坏形态与标准模型（图11（f））较相似，而其他模型左半侧断裂路径的发展与标准模型差异明显。正如第5节所分析的，标准模型在达到极限抗压强度时于右半侧形成了第一个剪切带的雏形（图11（b）），因而模型组右半侧的断裂带（FZ5—FZ8）成形更早，且断裂路径的发展更加稳定。由于左半侧的断裂带发展迟滞于右半侧，且石灰岩、大理岩和花岗岩骨料模型中的骨料均出现了不同程度的损伤和宏观断裂，因此各模型演化至最后的破坏形态在左半侧部分出现了明显的差异。

7 结论

本文在内聚力模型的基础上建立了混凝土拉压统一的细观数值模型，讨论了模型的网格敏感性问题，对混凝土进行了单轴拉压的全过程模拟分析，对断裂能、界面强度、骨料形状和骨料类型进行了参数研究，得到以下结论：（1）由于混凝土的受压性能取决于试件整体响应，过粗的网格划分会降低局部区域的变形能力，限制内力重分布，进而令极限抗压强度下降。由于混凝土的受拉性能取决于最弱截面，网格划分形式对单一截面面积的影响较小，因此网格敏感性较低。单元尺寸与模型边长尺寸之比为1/150时可在较经济的计算成本下取得较理想的模拟结果。（2）相比于受拉性能，混凝土的受压性能因骨料分布不同而表现出更明显的离散性。受拉分析中，骨料分布的影响主要集中在强度软化段的初期。在受压分析中，骨料分布会影响整个非线性受力阶段的受力性能。在混凝土的剩余抗压强度大致等于1/2 ～1/3倍的最大抗压强度时，主剪切带成型，混凝土会进入一个相对稳定的受力状态。（3）混凝土在单轴拉力作用下发生的是偏向于Ⅰ型断裂的混合型断裂。其中，Ⅰ型和Ⅱ型断裂能之差异所致的混合型断裂能大约占受拉总断裂能的一半。混凝土在单轴压力作用下发生的是偏向于Ⅱ型断裂的混合型断裂。提升Ⅱ型断裂能能够增强混凝土在压力作用下的变形和内力重分布能力，进而提升极限抗压强度和剩余抗压强度。（4）界面过渡区和砂浆内界面单元的抗拉和抗切强度均会影响混凝土的抗压强度。虽然整体上混凝土抗压强度取决于界面的剪切强度，但在较小强度范围内，单元的抗拉强度相比于剪切强度对混凝土的极限抗压强度的影响更为明显。（5）随着骨料边数的增加，混凝土的极限抗拉和抗压强度均会有一定的上升，其可能受到角缘效应和应力集中的影响。正五边形骨料模型的极限抗拉和抗压强度与任意多边形骨料模型相当。骨料边数越多可能会令混凝土受压开裂后的骨料咬合效应减弱，进而导致软化段的强度退化率增大。（6）在砂浆和界面过渡区强度小于骨料内界面强度的情况下，骨料的材料类型对混凝土受压性能的影响主要集中在强度软化段。在断裂演化方面，骨料的损伤和断裂会引起断裂路径的分化，进而导致破坏形态出现不同程度的差异。首条受压剪切带临近区域的断裂演化较为相似，而后成型的远端区域的破坏形态差异较为明显。