基于相场正则化黏聚区模型的混凝土混合型细观断裂行为

Mixed⁃mode mesoscale fracture behavior of concrete based on a phase field regularized cohesive zone model

摘要

为研究混凝土混合型断裂及损伤演化机理，分析了细观结构特性的影响机制，借助相场正则化黏聚区模型开展了数值仿真研究。结果表明：相场正则化黏聚区模型可准确模拟混凝土细观断裂行为，且不受相场尺度参数与网格尺寸影响。在混合型断裂过程中，预制切口尖端处首先萌生损伤并起裂，致使混凝土软化，裂缝呈一定角度扩展，最终结构失效。降低集料体积分数可提升抗裂性能，减少裂缝偏转次数；优化界面性能亦可起到改善作用，但其效果相对偏弱。

Abstract

To study the mechanism of mix-mode fracture and damage evolution in concrete， and analyze the influence of the mesoscale structural properties， a phase field regularized cohesive zone model （PF-CZM） was used to perform numerical simulations. The results indicated that the PF-CZM can accurately model the mesoscale fracture behavior of concrete， which is independent from the phase field scale parameter and meshing size. In the process of mix-mode fracture， damage and cracking occur at the notch tip， leading to concrete softening， cracks to expand at an angle， and eventual failure. Reducing the aggregate volume content can improve the cracking resistance. Optimizing the interfacial properties also contributes to these improvements， but the effect is weaker compared to a decrease in aggregate volume fraction.

Graphical abstract

关键词

道路工程 / 混凝土 / 相场正则化黏聚区模型 / 相场理论 / 混合型断裂

Key words

road engineering / concrete / phase-field regularized cohesive zone model / phase field theory / mixed-mode fracture

0 引 言

混凝土在荷载作用下易诱发开裂，进而导致服役结构性能下降^［1-4］。故此，揭示混凝土中裂缝起裂与演化的机理，对提高混凝土结构的耐久性和安全性具有重要意义。目前，已有学者通过室内试验的方式研究了混凝土的断裂行为^［5-7］，但难以观测材料中损伤和裂缝的演化过程，并且未能建立宏观断裂与细观各相组成特性之间的关系。此外，前人研究关注Ⅰ型和Ⅱ型断裂居多^［8-11］。而混凝土结构在服役过程中多处于拉-剪同时作用的复杂应力状态，属于混合型断裂，但目前该类型受关注较少。故此，有必要从细观尺度出发，借助数值仿真，对更符合真实受力状态的混凝土混合型断裂开展深入研究。

传统混凝土断裂仿真理论主要分为离散裂缝模型^［12］和弥散裂缝模型^［13］。前者通过在单元边界处设置界面单元表征开裂引起的强不连续性^［14］。而后者则基于连续力学理论，将开裂引起的强不连续性弥散至邻近单元内部予以近似表征^［15］。上述模型各有优点，但均具有较强的网格依赖性^{［16， 17］}。为解决该问题，后续有学者相继提出了裂缝带模型^［18］、梯度增强损伤模型^［19］等理论，并引入了扩展有限元^［20］、近场动力学^［21］等数值方法。但这些模型与方法在数值实现上存在一定难度或需定义烦琐的裂缝追踪算法和开裂准则。

相场模型是近年来发展较快的断裂理论模型。该模型将尖锐裂缝近似为一定宽度的正则化裂缝，并引入相场变量表征材料的损伤状态，且自带基于能量原理的裂缝扩展及方向判据，可有效弥补上述传统仿真模型的不足^{［22， 23］}。遗憾的是，相场模型目前仍主要应用于脆性断裂问题，且具有很高的相场尺度参数敏感性^［24］。直至2017年，Wu^［25］提出了相场正则化黏聚区模型（Phase field regularized cohesive zone model，PF-CZM），其统一了脆性和准脆性断裂问题，并消除了网格及相场尺度参数敏感性，已成功应用于玻璃、橡胶等材料的断裂仿真研究^［26-31］。但是，目前应用于混凝土混合型细观断裂的研究仍较少。

故此，本文拟基于有限元方法对PF-CZM进行数值实现，探讨其在细观尺度下模拟混凝土断裂的可行性，并揭示混凝土的混合型断裂行为及损伤演化机理，建立集料体积分数及界面力学性能与混凝土细观断裂的量化关系，以期为混凝土材料与工程结构的断裂破坏分析提供科学依据。

1 PF-CZM的理论背景

如图1所示，含尖锐裂缝的固体记为\Omega ={\mathbb{R}}^{n_{\mathrm{d}\mathrm{i}\mathrm{m}}}(n_{\mathrm{d}\mathrm{i}\mathrm{m}}=\mathrm{1,2},\mathrm{3})，其外边界为\partial \Omega \subset {\mathbb{R}}^{n_{\mathrm{d}\mathrm{i}\mathrm{m}}-\mathrm{1}}，外法向矢量为 n。为保证边值问题的适定性，外边界\partial \Omega被分为两个独立的部分\partial {\Omega }_{\mathrm{u}}和\partial {\Omega }_{\mathrm{t}}，并分别施加特定的位移边界 u^*（ x ）与力边界 t^*（ x ）。同时，在仅考虑小应变的情况下，位移场 u （ x ）和应变场 ε （ x ）可用于描述固体的变形状态。

1.1　尖锐裂缝的相场正则化

固体中的尖锐裂缝记为\mathrm{𝒮}\subset {\mathbb{R}}^{n_{\mathrm{d}\mathrm{i}\mathrm{m}}-\mathrm{1}}，其外法向矢量记为{\mathit{n}}_{\mathrm{𝒮}}，如图1（a）所示。在相场模型中，尖锐裂缝\mathrm{𝒮}被弥散为一定宽度的正则化裂缝\mathrm{ℬ}\subset \Omega，其宽度由相场尺度参数b确定，并引入相场变量d(\mathit{x})\in [\mathrm{0},\mathrm{1}]描述材料的损伤状态，如图1（b）所示。正则化裂缝\mathrm{ℬ}的外边界为\partial \mathrm{ℬ}，其外法向矢量为{\mathit{n}}_{\mathrm{ℬ}}。故此，根据图像分割理论，尖锐裂缝A_{\mathrm{𝒮}}可以近似地被正则化裂缝A_{\mathrm{d}}(d)表示^［23］：

\underset{\mathrm{尖锐}\mathrm{裂缝}}{\underbrace{A_{\mathrm{𝒮}}={\int }_{\mathrm{𝒮}}\mathrm{d}A}}\approx \underset{\mathrm{正则}\mathrm{化裂}\mathrm{缝}}{\underbrace{{\int }_{\mathrm{ℬ}}\gamma (d,\nabla d)\mathrm{d}V=A_{\mathrm{d}}(d)}}（1）

式中：γ为裂缝面积密度泛函数，其表达式为：

\left\{\begin{array}{l}\gamma (d,\nabla d)=\frac{\mathrm{1}}{c_{\alpha }}\left[\frac{\mathrm{1}}{b}\alpha (d)+b{\left|\nabla d\right|}^{\mathrm{2}}\right]\\ c_{\alpha }=\mathrm{4}{\int }_{\mathrm{0}}^{\mathrm{1}}\sqrt[]{\alpha (\beta )}\mathrm{d}\beta \end{array}\right.（2）

式中：\alpha (d)\in [\mathrm{0},\mathrm{1}]为裂缝几何函数；\nabla d(x)为相场变量d(x)的梯度；c_{\alpha }为与裂缝几何函数有关的缩放系数。可以发现，当相场尺度参数b趋近于0时，式（1）中的正则化裂缝A_{\mathrm{d}}(d)收敛于尖锐裂缝A_{\mathrm{𝒮}}^［25］。

同时，为保证固体断裂过程中局部能量耗散逐渐递增，并明显区分d=0（未受损伤）和d=1（完全损伤）两种损伤状态，裂缝几何函数\alpha (d)应满足如下条件^［25］：

\alpha \text{'}(d)\ge \mathrm{0},\alpha (\mathrm{0})=\mathrm{0},\alpha (\mathrm{1})=\mathrm{1}（3）

另一方面，固体在发生开裂后，其Helmholtz自由能\psi会出现退化，如式（4）所示：

\psi (\epsilon ,d)=\psi (\epsilon ,\omega (d))（4）

式中：\omega (d)\in [\mathrm{0,1}]为能量退化函数。

故此，为确保固体断裂过程中裂缝宽度不会持续增大，Miehe等^［23］提出能量退化函数，满足如下条件：

\left\{\begin{array}{l}\omega (\mathrm{0})=\mathrm{1},\omega (\mathrm{1})=\mathrm{0}\\ \omega \text{'}(d)\le \mathrm{0},\omega \text{'}(\mathrm{1})=\mathrm{0}\end{array}\right.（5）

1.2　控制方程

在PF-CZM中，根据热力学能量等效原理，裂缝相场d（ x ）的强形式控制方程为^［25］：

\left\{\begin{array}{l}\nabla \cdot \mathit{q}+Q\le \mathrm{0},\mathrm{i}\mathrm{n}\mathrm{ℬ}\\ \mathit{q}\cdot {\mathit{n}}_{\mathrm{B}}\ge \mathrm{0},\mathrm{o}\mathrm{n}\partial \mathrm{ℬ}\end{array}\right.（6）

式中：热力学广义 q 和源项Q为：

\left\{\begin{array}{l}\mathit{q}=\frac{\mathrm{2}b}{c_{\alpha }}{G}_{\mathrm{f}}\nabla d\\ Q=-\omega \mathrm{\text{'}}(d)\overline{Y}-\alpha \mathrm{\text{'}}(d)\frac{G_{\mathrm{f}}}{c_{\alpha }b}\end{array}\right.（7）

式中：Y为裂缝驱动力，也称为能量释放率；\overline{Y}为有效裂缝驱动力；G_f为材料的断裂能。

在不考虑材料愈合的情况下，裂缝扩展具有不可逆性。故此，Miehe等^［32］引入历史变量对源项Q进行修正，如式（8）所示：

\left\{\begin{array}{l}Q=-\omega \mathrm{\text{'}}(d)\mathrm{\mathscr{H}}-\alpha \mathrm{\text{'}}(d)\frac{G_{\mathrm{f}}}{c_{\alpha }b}\\ \mathrm{\mathscr{H}}=\mathrm{m}\mathrm{a}\mathrm{x}({\overline{Y}}_{\mathrm{0}},\underset{n\in [\mathrm{0},T]}{\mathrm{m}\mathrm{a}\mathrm{x}}{\overline{Y}}_n)\end{array}\right.（8）

式中：\mathrm{\mathscr{H}}为有效裂缝驱动力\overline{Y}的历史最大变量；{\overline{Y}}_{\mathrm{0}}为有效裂缝驱动力的初始阈值。当\overline{Y}>{\overline{Y}}_{\mathrm{0}}时，材料发生开裂。通过该历史变量的定义，控制方程式（6）也将转化为等式，从而可以通过常规求解器进行计算。

另一方面，根据经典的力学平衡方程，位移场 u （ x ）的强形式控制方程为^［25］：

\left\{\begin{array}{l}\nabla \cdot \mathit{\sigma }+{\mathit{b}}^{\mathrm{*}}=\mathrm{0},\mathrm{i}\mathrm{n}\Omega \\ \mathit{\sigma }\cdot \mathit{n}={\mathit{t}}^{\mathrm{*}},\mathrm{o}\mathrm{n}\partial \Omega \end{array}\right.（9）

式中： σ 为应力张量； b^* 为体积均布力。

1.3　本构关系

在上述控制方程中还需补充应力张量 σ 和有效裂缝驱动力\overline{Y}的本构关系。考虑混凝土在拉、压状态下的不同力学行为，基于Rankine破坏准则，定义其本构关系为^［25］：

\left\{\begin{array}{l}\mathit{\sigma }=\omega (d)\overline{\mathit{\sigma }}=\omega (d)E_{\mathrm{0}}:\mathit{\epsilon }\\ Y=-\omega \text{'}(d)\overline{Y}\\ \overline{Y}=\frac{\mathrm{1}}{\mathrm{2}{E}_{\mathrm{0}}}〈{\overline{\sigma }}_{\mathrm{1}}〉\\ {\overline{Y}}_{\mathrm{0}}=\frac{f_{\mathrm{t}}^{\mathrm{2}}}{\mathrm{2}{E}_{\mathrm{0}}}\end{array}\right.（10）

式中：〈x〉=\mathrm{m}\mathrm{a}\mathrm{x}(x,\mathrm{0})为McAuley括号；{\overline{\sigma }}_{\mathrm{1}}为最大有效主应力；E₀和f_t分别为材料的杨氏模量及抗拉强度。

1.4　特征函数

由前文可知，PF-CZM包括两个重要的特征函数，即裂缝几何函数\alpha (d)与能量退化函数\omega (d)，而这两类函数正是诸多相场模型之间的主要区别。在PF-CZM中，Wu^［25］采用参数化的多项式函数作为特征函数，如式（11）所示。此次改进较好地解决了计算结果对相场尺度参数b的敏感性，并统一了脆性断裂与准脆性断裂问题。

\left\{\begin{array}{l}\alpha (d)=\mathrm{2}d-d^{\mathrm{2}}\\ \omega (d)=\frac{{(\mathrm{1}-d)}^p}{{(\mathrm{1}-d)}^p+a_{\mathrm{1}}d\cdot P(d)}\\ c_{\alpha }=\mathrm{\pi }\\ P(d)=\mathrm{1}+a_{\mathrm{2}}d+a_{\mathrm{3}}{d}^{\mathrm{2}}\end{array}\right.（11）

式中：p、a_{\mathrm{2}}和a_{\mathrm{3}}通过材料的软化曲线予以确定；a_{\mathrm{1}}则通过式（12）计算：

\left\{\begin{array}{l}{a}_{\mathrm{1}}=\frac{\mathrm{4}}{\mathrm{\pi }}\cdot \frac{l_{\mathrm{c}\mathrm{h}}}{b}\\ l_{\mathrm{c}\mathrm{h}}=\frac{E_{\mathrm{0}}{G}_{\mathrm{f}}}{f_{\mathrm{t}}^{\mathrm{2}}}\end{array}\right.（12）

式（13）中给出了常见的4种软化曲线参数取值，其有效性已在文献［25，33-36］中得以验证。对于混凝土材料，一般建议采用Cornelissen软化曲线。

\left\{\begin{array}{l}\mathrm{线性}:p=\mathrm{2.0},a_{\mathrm{2}}=-\frac{\mathrm{1}}{\mathrm{2}},a_{\mathrm{3}}=\mathrm{0}\\ \mathrm{指数}:p=\mathrm{2.5},a_{\mathrm{2}}={\mathrm{2}}^{\mathrm{5}/\mathrm{3}}-\mathrm{3},a_{\mathrm{3}}=\mathrm{0}\\ \mathrm{双曲}\mathrm{线}:p=\mathrm{4.0},a_{\mathrm{2}}={\mathrm{2}}^{\mathrm{7}/\mathrm{3}}-\frac{\mathrm{9}}{\mathrm{2}},a_{\mathrm{3}}=\mathrm{0}\\ \mathrm{C}\mathrm{o}\mathrm{r}\mathrm{n}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{e}\mathrm{n}:p=\mathrm{2.0},a_{\mathrm{2}}=\mathrm{1.386}\mathrm{8},a_{\mathrm{3}}=\mathrm{0.910}\mathrm{6}\end{array}\right.（13）

2 PF-CZM的数值实现

2.1　求解模块

本文基于COMSOL Multiphysics平台，通过有限元方法对PF-CZM进行数值实现。由前文所述，裂缝驱动力是确定裂缝起裂、扩展、演化的重要指标。本文采用状态变量（State variables）模块计算每个增量中初始时的裂缝驱动力^［37］，其计算公式为：

\mathrm{\mathscr{H}}=\mathrm{m}\mathrm{a}\mathrm{x}\left[{\overline{Y}}_{\mathrm{0}},\mathrm{m}\mathrm{a}\mathrm{x}({\overline{\mathit{𝒴}}}_n,\overline{Y})\right]（14）

式中：{\overline{\mathit{𝒴}}}_n为直至最后一个增量所达到的历史最大裂缝驱动力；\overline{Y}和{\overline{Y}}_{\mathrm{0}}分别为当前增量的有效裂缝驱动力及其阈值。在当前增量末尾时，最大裂缝驱动力\overline{\mathit{𝒴}}被式（15）更新，并应用于下一个增量中的计算。

\overline{\mathit{𝒴}}=\mathrm{m}\mathrm{a}\mathrm{x}({\overline{\mathit{𝒴}}}_n,\overline{Y})（15）

此外，本文采用固体力学（Solid mechanics）和泊松方程（Poisson's equation）模块分别对PF-CZM的控制方程式（9）和（6）进行定义。在固体力学模块中，应力-应变本构关系式（10）通过线弹性材料模型予以表征，定义杨氏模量为\omega (d)E_{\mathrm{0}}，泊松比为{\nu }_{\mathrm{0}}。在泊松方程模块中，根据式（6），定义扩散系数c和源项f为式（16）。之后，COMSOL Multiphysics会自行通过加权余量法将上述定义的强形式控制方程及本构关系转换为便于有限元离散求解的等效弱形式控制方程。

\left\{\begin{array}{l}c=\frac{\mathrm{2}b}{c_{\alpha }}{G}_{\mathrm{f}}\\ f=-\omega \text{'}(d)\mathrm{\mathscr{H}}-\frac{G_{\mathrm{f}}}{c_{\alpha }b}\alpha \text{'}(d)\end{array}\right.（16）

2.2　求解算法

在完成PF-CZM的数值实现后，对该类强场耦合问题，一般可选择整体全耦合算法（Monolithic algorithm）或分离交错迭代算法（Staggered algorithm）进行求解。但COMSOL Multiphysics中内嵌的基于Newton-Raphson的整体算法在数值健壮性方面有所欠缺^［38］。故此，本文采用交错迭代求解法，将节点未知变量分为若干组，每次仅求解一组中的场变量，其他场变量被固定。然后，该组变量根据上次计算结果更新，第二组的场变量被求解，其他场变量被固定。以此类推，直至最后的解收敛。该方法虽收敛速度较慢，但具有较好的数值健壮性，对求解纯力学场的断裂问题十分有效。

2.3　模型验证

为验证上述模型与求解算法的可靠性，本文采用单轴拉伸的经典案例进行计算，并将模拟结果与解析解^［39］进行比较，其几何尺寸及边界条件如图2所示。试件尺寸为200 mm × 10 mm × 10 mm，左侧被固定，右侧被施加单调递增的位移荷载。材料属性选取混凝土的典型参数，杨氏模量为30 GPa，泊松比为0.2，抗拉强度为3.0 MPa，断裂能为0.12 N/mm。软化曲线选用Cornelissen曲线，即参数取值为p=2，a₂=1.386 8，a₂=0.656 7。同时，为讨论相场尺度参数b与网格尺寸h对计算结果的影响，本节选取4个相场尺度参数b取值：5、10、20、50 mm，以及4个网格尺寸h取值：b、5b、10b、20b，并选用四边形单元对模型进行离散。此外，试件两端被赋予d=0的Dirichlet边界条件，以避免试件在边缘产生损伤，导致无法完整获取致使试件断裂的断裂能^［35］。

图3为模拟得到的荷载-位移曲线。可以发现，当b/h=1时，计算结果在弹性阶段和峰值荷载处与解析解吻合度较高，但在软化阶段产生了不同的下降趋势。分析认为是网格尺寸h较大造成正则化裂缝所弥散的损伤过宽，导致计算结果失真。而当b/h≥5时，不同网格尺寸h与相场尺度参数b对应的荷载-位移曲线基本相同，且均与解析解具有较高的一致性。这表明，只要相场尺度参数b与网格尺寸h的比值大于等于5，PF-CZM就可精准地模拟混凝土的准脆性断裂，且计算结果不受网格尺寸h和相场尺度参数b的影响。

3 混凝土细观混合型断裂模拟

3.1　细观数值模型建立

为揭示混凝土在细观尺度下的混合型断裂行为，本文引入Gálvez等^［40］所开展的混合型三点弯曲小梁试验进行模拟，其几何尺寸与边界条件如图4所示。试件尺寸为337.5 mm × 75 mm × 50 mm，跨中预制长37.5 mm、宽2 mm的缺口。试件底部设置两个不对称的、跨度为187.5 mm的约束，并在顶部施加一个距中心线84.975 mm的偏心荷载，以实现混合型断裂。试件集料级配取自文献［41］，如表1所示，其体积分数为30%。

本文参照文献［42-44］的方法首先构建含多边形随机骨料的虚拟试件。该方法通过Matlab中内嵌的Rand随机函数及Convhull凸包函数随机生成不规则多边形骨料，并在骨料周围形成一定宽度的缓冲区，然后将生成的骨料投放至矩形区域内，且指定相邻骨料的缓冲区之间没有相交，以确保骨料之间的最小间距。然而，上述方法未考虑界面过渡区（Interfacial transition zones，ITZs）的生成，所以本文在此基础上，将缓冲区宽度定义为ITZs的厚度加上骨料之间的最小间距（见图5），并保留ITZs区域，以建立包含骨料、砂浆、ITZs的三相混凝土细观数值模型。其中，相邻骨料之间的最小间距取0.1 mm。而ITZs的厚度根据已有试验表明通常在10~50 μm^［45］，但为提高计算效率，本文假定ITZs的厚度为0.5 mm。同时，为避免计算结果出现偶然性，本节建立3个虚拟试件进行计算，其中虚拟试件3的数值模型如图6所示。

细观数值模型中，水泥砂浆的材料属性取Gálvez等^［40］的试验数据，但ITZs的材料属性通常难以通过试验测定，故本文参考前人的试验与数值研究^［46-49］，取其杨氏模量、抗拉强度、断裂能为水泥砂浆的75%，如表2所示。集料被假定为线弹性体，其杨氏模量取70 GPa，并认为不会发生开裂。各相的泊松比取混凝土的经验值0.2。水泥砂浆及ITZs采用Cornelissen软化曲线表征，即参数取值为p=2，a₂=1.3868，a₂=0.6567。

此外，前文已证明，当相场尺度参数b与网格尺寸h的比值大于等于5时，PF初CZM可以获得准确解。因此，本文在满足精度要求的同时，为更明显地观察裂缝起裂、扩展及演化的过程，选取网格尺寸h为0.5 mm，相场尺度参数b为2.5 mm（5h），并采用三角形单元对模型进行离散。同时，为提高计算效率，本节对集料部分采用稍作放大的网格尺寸。最终模型的网格分布如图7所示。

对荷载施加，本文采用控制裂缝张口位移（Crack mouth opened displacement，CMOD）的间接位移加载方法。相比直接位移加载，该方法可更为完整地捕捉开裂过程，避免在某一步长内出现较快的裂缝扩展，从而导致荷载-位移或荷载-CMOD曲线的部分信息缺失^［50］。

3.2　混合型断裂行为及损伤演化机理分析

图8为计算得到的荷载-CMOD曲线。结果表明：不同虚拟试件之间的荷载-CMOD曲线几乎完全相同，且与试验数据吻合度较高，仅在峰值荷载处及软化阶段存在细微差异。此外，当CMOD较小时，荷载-CMOD曲线呈线性增长，该阶段可近似地将混凝土视作为线弹性体。随着CMOD继续增大，曲线斜率开始下降，进入非线性阶段。分析认为是混凝土内部有损伤（微裂缝）产生，致使试件性能退化。随后，曲线达到峰值荷载，并随着CMOD的增大，荷载开始减小，进入软化阶段。这是因为混凝土内部损伤恶化，并发生宏观开裂，导致试件失稳，抵御荷载能力下降。最后，荷载随CMOD的增大而持续减小，直至试件完全失效。

图9则为不同虚拟试件的裂缝扩展路径。可见，不同虚拟试件之间存在一定差异，但均与试验数据接近，其主要原因为开裂区域内不同的集料含量和形状导致，该影响因素也将在后文进行详细讨论。此外，在混凝土混合型断裂中，试件会在缺口尖端处由于应力集中作用首先发生开裂，并且裂缝在拉-剪应力同时作用下沿着ITZs区域逐渐向顶部扩展，最终形成与预制缺口呈一定角度的裂缝路径。

该数值算例再次验证了本文所建立的PF-CZM模拟混凝土细观混合型断裂行为的可靠性，并为后文研究细观组成特性对混凝土断裂的影响机制奠定了扎实的理论基础。

为进一步分析混凝土在混合型断裂过程中的损伤及裂缝演化机理，本节以虚拟试件3为例，选取弹性阶段、非线性阶段、软化阶段和失效阶段中的代表性时刻进行详细讨论，分别为图8中点A（CMOD=0.007 mm）、点B（CMOD=0.025 mm）、点C（CMOD=0.1 mm）和点D（CMOD=0.5 mm）。图10为上述代表性时刻下混凝土中的损伤分布图。结果表明：当CMOD较小时，混凝土处于弹性阶段，内部几乎没有损伤，仍保持较好的完整性，如图10（a）所示。当CMOD逐渐增大，混凝土进入非线性阶段，首先在缺口尖端处萌生损伤，即微裂缝，试件性能开始衰减，但仍无宏观裂缝形成，如图10（b）所示。随后，CMOD继续增大，混凝土进入软化阶段，在缺口尖端处发生起裂，裂缝沿着薄弱的ITZs区域扩展，导致结构失稳，如图10（c）所示。最后，裂缝在荷载的持续作用下继续向小梁顶部扩展，直至贯穿整个试件，导致结构完全失效，丧失抵御荷载能力，如图10（d）所示。需注意的是，ITZs由于其较低的抗拉强度及断裂能，是裂缝扩展的主要区域。故此，后文将进一步讨论ITZs界面力学性能对混凝土细观开裂的影响。

4 细观组成特性对混凝土混合型断裂的影响

4.1　集料体积分数的影响

为了建立集料体积分数与混凝土细观断裂行为之间的关系，本节在3.1节的基础上，模拟了不同体积分数下的混凝土荷载-CMOD曲线，并分析了体积分数对峰值荷载与峰值CMOD的影响，如图11和图12所示。

结果表明：集料体积分数对混凝土的细观断裂行为具有较大的影响。具体来说，在荷载-CMOD曲线的弹性阶段，随着集料体积分数增加，曲线斜率会增大，即混凝土弹性模量增大。这是因为，集料的弹性模量远大于水泥砂浆，随着集料含量的增多，混凝土的弹性模量也会随之增大。但是，集料体积分数的增加反而会导致峰值荷载和与之对应的峰值CMOD近似呈线性关系减小。分析认为是随着集料含量增多，ITZs的体积占比也会提高，从而增加了混凝土内的断裂薄弱区，导致抗裂性能下降。由此可见，适当增加混凝土的砂浆含量，减小集料体积分数，可以实现混凝土峰值荷载和峰值CMOD的双提升，进而获得更为优异的抗裂性能。

此外，图13呈现了不同体积分数下混凝土的裂缝扩展路径。由此可见，当体积分数较小时，裂缝扩展路径较为平滑，如图13（a）所示。随着体积分数的增加，裂缝的曲折程度也会随之增加，并且会改变裂缝扩展的方向，如图13（c）~（d）所示。这是因为，集料通常不会发生开裂，但集料周围的ITZs具有较弱的抗裂性能。因此，当裂缝扩展至集料处时，会沿着集料边缘的ITZs继续演化，导致裂缝路径发生偏转，并且偏转方向与集料的形状有关。故此，集料的体积分数越高，裂缝在扩展过程中发生的偏转次数则越多，进而导致最终形成的裂缝路径就越为曲折。

4.2　ITZs界面力学性能的影响

为了探讨ITZs界面力学性能对混凝土细观断裂的影响规律，本节在3.1节的基础上，取ITZs的抗拉强度和断裂能分别为50%、75%、100%的砂浆，并以体积分数为30%的虚拟试件为例，绘制了不同ITZs界面力学性能下的荷载-CMOD曲线，如图14所示。可见，随着ITZs界面力学性能的下降，荷载-CMOD曲线在弹性阶段与软化阶段变化较小，但在峰值荷载处产生了较大的差异。

故此，本节进一步分析了不同集料体积分数下ITZs界面力学性能与峰值荷载和峰值CMOD的关系，以揭示两者之间的耦合作用机制，如图15所示。结果表明：在不同的体积分数下，ITZs界面力学性能对峰值荷载与峰值CMOD产生了不同程度的影响。对峰值荷载，在不同的体积分数下，其与ITZs界面力学性能近似呈线性增长关系，且增长斜率随着体积分数的增加而增加。这是因为，ITZs是混凝土发生开裂的首要区域，当该处的力学性能提升时，混凝土的断裂阈值会随之提高，从而可承受更高的峰值荷载。并且，在集料体积分数较大时，ITZs亦具有较高的体积占比，此时提升界面力学性能带来的改善效果则更加显著。

另一方面，不同体积分数下峰值CMOD与ITZs界面力学性能的线性拟合曲线斜率均在0±0.000 12以内，偏差部分分析认为是由模型随机性及数值误差导致。故此，可以认为界面力学性能对混凝土峰值CMOD几乎没有影响，而主要由体积分数控制。综上所述，改善水泥砂浆与集料之间ITZs的界面力学性能是提升混凝土抗裂强度的关键，但对峰值CMOD影响较小，其主要由水泥砂浆含量及性能决定。

此外，以体积分数为30%的虚拟试件为例，图16绘制了不同ITZs界面力学性能下的裂缝扩展路径。结果表明：ITZs界面力学性能越差，裂缝发生偏转的次数越多，裂缝扩展路径越曲折。以图16中A处为例，当ITZs界面力学性能较低时，该处会发生开裂，如图16（a）所示；而当界面力学性能提升时，该处的断裂阈值会随之提高，在荷载作用下会产生部分损伤，但不会形成宏观裂缝，进而不会导致裂缝偏转，如图16（b）和（c）所示。此外，当界面力学性能较低时，邻近裂缝主路径的ITZs会萌生细微损伤，如图16（a）中B处所示。该处损伤虽不会演化为宏观裂缝，但会致使混凝土整体性能下降，并且在二次加载时存在较大的开裂隐患。

5 结 论

（1）基于有限元法对PF-CZM完成了数值实现，并通过随机算法构建了混凝土的细观数值模型，结合解析解及混合型三点弯曲小梁试验验证了模型的准确性。结果表明：当相场尺度参数b与网格尺寸h的比值大于等于5时，数值模拟结果与解析解及试验数据均基本吻合，PF-CZM可实现对混凝土裂缝与损伤演化全过程的精准模拟以及抗裂性能的定量预测，并且不受相场尺度参数b和网格尺寸h的影响。

（2）混凝土发生混合型断裂的过程中，试件在弹性阶段几乎没有损伤产生，结构保持完整，荷载-CMOD曲线呈线性增长。随着CMOD逐渐增大，预制缺口尖端处开始萌生损伤并发生开裂，裂缝在拉-剪应力同时作用下扩展，致使混凝土软化。最终，裂缝扩展至小梁顶部，结构完全失效，试件丧失抵御荷载能力。

（3）细观尺度下，降低集料体积分数可实现混凝土峰值荷载与峰值CMOD的双提升，并可获得更为平滑的裂缝路径。而提升ITZs的界面力学性能亦可提高整个试件的断裂阈值，减少裂缝偏转次数，但其改善效果要弱于体积分数的下降，且对混凝土峰值CMOD影响较小。

参考文献

[1] Wang Y, Zhang B, Gao S H, et al. Investigation on the effect of freeze-thaw on fracture mode classification in marble subjected to multi-level cyclic loads[J]. Theoretical and Applied Fracture Mechanics, 2021, 111: No.102847.

[2] Zhang K, Yuan Q, Huang T, et al. Predicting the cracking behavior of early-age concrete in CRTS III track[J]. Construction and Building Materials, 2022, 353: No.129105.

[3] 梁宁慧, 缪庆旭, 刘新荣, 等. 聚丙烯纤维增强混凝土断裂韧度及软化本构曲线确定[J]. 吉林大学学报:工学版, 2019, 49(4): 1144-1152.

[4] Liang Ning-hui, Miao Qing-xu, Liu Xin-rong, et al. Determination of fracture toughness and softening traction⁃separation law of polypropylene fiber reinforced concrete[J]. Joural of Jilin University (Engineering and Technology Edition), 2019, 49(4): 1144-1152.

[5] 左新黛, 张劲泉, 赵尚传. 在役混凝土T梁疲劳刚度退化及寿命预测方法[J]. 吉林大学学报:工学版, 2023, 53(9): 2563-2572.

[6] Zuo Xin-dai, Zhang Jin-quan, Zhao Shang-chuan. Fatigue stiffness degradation and life prediction method of in⁃service concrete T⁃beams[J]. Joural of Jilin University(Engineering and Technology Edition), 2023, 53(9): 2563-2572.

[7] Grégoire D, Rojas-Solano L B, Pijaudier-Cabot G. Failure and size effect for notched and unnotched concrete beams[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2013, 37(10): 1434-1452.

[8] Hordijk D A. Local approach to fatigue of concrete [D]. Delft: Faculty of Civil Engineering and Geosciences, Delft University of Technology, 1991.

[9] Zhang H Z, Xu Y D, Gan Y D, et al. Experimentally validated meso-scale fracture modelling of mortar using output from micromechanical models[J]. Cement & Concrete Composites, 2020, 110: No.103567.

[10] Chen Y W, Feng J L, Li H, et al. Effect of coarse aggregate volume fraction on mode Ⅱ fracture toughness of concrete[J]. Engineering Fracture Mechanics, 2021, 242: No.107472.

[11] Yang D, Dong W, Liu X F, et al. Investigation on mode-I crack propagation in concrete using bond-based peridynamics with a new damage model[J]. Engineering Fracture Mechanics, 2018, 199: 567-581.

[12] 童谷生, 姚良发, 徐攀. C40自密实混凝土Ⅱ型断裂韧度的确定[J]. 力学季刊, 2021, 42(4): 800-810.

[13] Tong Gu-sheng, Yao Liang-fa, Xu Pan. Determination of type Ⅱ fracture toughness of C40 self-compacting concrete[J]. Chinese Quarterly of Mechanics, 2021, 42(4): 800-810.

[14] 王康, 吴佰建, 李兆霞. 损伤跨尺度演化导致的混凝土强度尺寸效应[J]. 东南大学学报:自然科学版, 2014, 44(6): 1230-1234.

[15] Wang Kang, Wu Bai-jian, Li Zhao-xia. Size effect on concrete strength caused by trans-scale damage evolution[J]. Journal of Southeast University(Natural Science Edition), 2014, 44(6): 1230-1234.

[16] Ngo D, Scordelis A C. Finite element analysis of reinforced concrete beams[J]. ACI Journal Proceedings, 1967, 64(3): 152-163.

[17] Rashid Y R. Ultimate strength analysis of prestressed concrete pressure vessels[J]. Nuclear Engineering and Design, 1968, 7(4): 334-344.

[18] 周正峰, 康玉峰, 罗君豪, 等. 基于黏聚区模型的混凝土细观断裂分析[J]. 东南大学学报: 自然科学版, 2021, 51(2): 270-277.

[19] Zhou Zheng-feng, Kang Yu-feng, Luo Jun-hao,et al. Mesoscale fracture analysis on concrete based on cohesive zone model[J]. Journal of Southeast University(Natrural Science Edition), 2021, 51(2): 270-277.

[20] Cope R J, Rao P V, Clark L A, et al. Modelling of reinforced concrete behaviour for finite element analysis of bridge slabs[C]∥Proceedings of the International Conference on Numerical Methods for Nonlinear Problems, Swansea, UK, 1981:457-470.

[21] Rots J G, Nauta P, Kuster G M A, et al. Smeared crack approach and fracture localization in concrete[J]. Heron, 1985, 30(1): 1-48.

[22] Yang Z J, Chen J F. Finite element modelling of multiple cohesive discrete crack propagation in reinforced concrete beams[J]. Engineering Fracture Mechanics, 2005, 72(14): 2280-2297.

[23] Bažant Z P, Oh B H. Crack band theory for fracture of concrete[J]. Matériaux Et Construction, 1983, 16(3): 155-177.

[24] Peerlings R H J, de Borst R, Brekelmans W A M, et al. Gradient enhanced damage for quasi-brittle materials [J]. International Journal for Numerical Methods in Engineering, 1996, 39(19): 3391-3403.

[25] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing[J]. International Journal for Numerical Methods in Engineering, 1999, 46(1): 131-150.

[26] Li W, Guo L. A mechanical-diffusive peridynamics coupling model for meso-scale simulation of chloride penetration in concrete under loadings[J]. Construction and Building Materials, 2020, 241: No.118021.

[27] Bourdin B, Francfort G A, Marigo J J. Numerical experiments in revisited brittle fracture[J]. Journal of the Mechanics and Physics of Solids, 2000, 48(4): 797-826.

[28] Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations[J]. International Journal for Numerical Methods in Engineering, 2010, 83(10): 1273-1311.

[29] Wu J Y, Nguyen V P, Nguyen C T, et al. Chapter One - Phase-field modeling of fracture[M]//Bordas S P A, Balint D S. Advances in Applied Mechanics. Amsterdam:Elsevier, 2020.

[30] Wu J Y. A unified phase-field theory for the mechanics of damage and quasi-brittle failure[J]. Journal of the Mechanics and Physics of Solids, 2017, 103: 72-99.

[31] Loew P J, Poh L H, Peters B, et al. Accelerating fatigue simulations of a phase-field damage model for rubber[J]. Computer Methods in Applied Mechanics and Engineering, 2020, 370:No.113247.

[32] Marbceuf A, Bennani L, Budinger M, et al. Electromechanical resonant ice protection systems: numerical investigation through a phase-field mixed adhesive/brittle fracture model[J]. Engineering Fracture Mechanics, 2020, 230:No. 106926.

[33] Wu J Y, Huang Y L, Nguyen V P, et al. Crack nucleation and propagation in the phase-field cohesive zone model with application to Hertzian indentation fracture[J]. International Journal of Solids and Structures, 2022, 241: No.111462.

[34] Yin B B, Zhang L W. Phase field method for simulating the brittle fracture of fiber reinforced composites[J]. Engineering Fracture Mechanics, 2019, 211: 321-340.

[35] Zhang P, Hu X F, Wang X Y, et al. An iteration scheme for phase field model for cohesive fracture and its implementation in Abaqus[J]. Engineering Fracture Mechanics, 2018, 204: 268-287.

[36] Zhang P, Hu X F, Yang S T, et al. Modelling progressive failure in multi-phase materials using a phase field method[J]. Engineering Fracture Mechanics, 2019, 209: 105-124.

[37] Miehe C, Welschinger F, Hofacker M. A phase field model of electromechanical fracture[J]. Journal of the Mechanics and Physics of Solids, 2010, 58(10): 1716-1740.

[38] Feng D C, Wu J Y. Phase-field regularized cohesive zone model (CZM) and size effect of concrete[J]. Engineering Fracture Mechanics, 2018, 197: 66-79.

[39] Wu J Y. A geometrically regularized gradient-damage model with energetic equivalence[J]. Computer Methods in Applied Mechanics and Engineering, 2018, 328: 612-637.

[40] Wu J Y, Mandal T K, Nguyen V P. A phase-field regularized cohesive zone model for hydrogen assisted cracking[J]. Computer Methods in Applied Mechanics and Engineering, 2020, 358:No. 112614.

[41] Wu J Y, Vinh Phu N. A length scale insensitive phase-field damage model for brittle fracture[J]. Journal of the Mechanics and Physics of Solids, 2018, 119: 20-42.

[42] Chen W X, Wu J Y. Phase-field cohesive zone modeling of multi-physical fracture in solids and the open-source implementation in Comsol Multiphysics[J]. Theoretical and Applied Fracture Mechanics, 2022, 117: No.103153.

[43] Ambati M, Gerasimov T, de Lorenzis L. A review on phase-field models of brittle fracture and a new fast hybrid formulation[J]. Computational Mechanics, 2015, 55(2): 383-405.

[44] Wu J Y. Unified analysis of enriched finite elements for modeling cohesive cracks[J]. Computer Methods in Applied Mechanics and Engineering, 2011, 200(45): 3031-3050.

[45] Gálvez J C, Elices M, Guinea G V, et al. Mixed mode fracture of concrete under proportional and nonproportional loading[J]. International Journal of Fracture, 1998, 94(3): 267-284.

[46] Yu K, Yang Z, Li H, et al. A mesoscale modelling approach coupling SBFEM, continuous damage phase-field model and discrete cohesive crack model for concrete fracture[J]. Engineering Fracture Mechanics, 2023, 278: No.109030.

[47] Dong Q, Zhao X, Chen X, et al. Numerical simulation of mesoscopic cracking of cement-treated base material based on random polygon aggregate model [J]. Journal of Traffic and Transportation Engineering (English Edition), 2023, 10(3): 454-468.

[48] 王国彤. 钢渣透水水泥混凝土力学性能的多尺度评价研究[D]. 南京: 东南大学交通学院, 2021.

[49] Wang Guo-tong. Multi-scale research on the mechanical properties of steel slag pervious cement concrete[D]. Nanjing: School of Transportation, Southeast University, 2021.

[50] 袁嘉伟. 基于内聚力模型的水稳碎石材料半圆弯曲开裂细观数值模拟研究[D]. 南京: 东南大学交通学院, 2020.

[51] Yuan Jia-wei. Research on meso-scale numerical simulation of cement treated base materials semi-circular bending cracking based on cohesive zone model[D]. Nanjing: School of Transportation, Southeast University, 2020.

[52] Xiao J, Li W, Corr D J, et al. Effects of interfacial transition zones on the stress–strain behavior of modeled recycled aggregate concrete[J]. Cement and Concrete Research, 2013, 52: 82-99.

[53] Gao J M, Qian C X, Liu H F, et al. ITZ microstructure of concrete containing GGBS[J]. Cement and Concrete Research, 2005, 35(7): 1299-1304.

[54] Huang Y, Yang Z, Ren W, et al. 3D meso-scale fracture modelling and validation of concrete based on in-situ X-ray computed tomography images using damage plasticity model[J]. International Journal of Solids and Structures, 2015, 67-68: 340-352.

[55] Mondal P. Nanomechanical Properties of Cementitious Materials[M]. Evanston: Northwestern University, 2008.

[56] Naderi S, Zhang M. A novel framework for modelling the 3D mesostructure of steel fibre reinforced concrete[J]. Computers & Structures, 2020, 234: No.106251.

[57] de Borst R. Computation of post-bifurcation and post-failure behavior of strain-softening solids[J]. Computers & Structures, 1987, 25(2): 211-224.

基金资助

国家自然科学基金项目(51978163)

国家自然科学基金项目(52208439)

江苏省自然科学基金项目(BK20200468)