基于相场方法的多场耦合模型与早期大体积混凝土热裂纹演化机制

《Results in Engineering》：Multi-field Coupling Model and Evolution Mechanism of Thermal Cracking in Early-age Mass Concrete Based on Phase-field Method

本研究提出了一种热-力-相场耦合模型，用于模拟大体积混凝土施工期的温度裂缝。模型结合了 hydration heat、creep和autogenous shrinkage等时间依赖性参数，通过增量迭代算法求解多场耦合方程，并通过基准案例和实际工程验证了其准确性，为裂缝预测和防治提供了理论支持。

　　在大型水利工程中，早期混凝土结构在施工过程中容易发生热裂缝，这会显著影响结构的安全性和耐久性。本研究建立了一个基于热力学理论框架和相场方法的热-机械-相场多场耦合模型，该模型考虑了水化热、蠕变特性以及自生体积变形等随时间变化的材料属性。模型通过增量迭代算法求解，并通过经典基准案例验证了其准确性。利用该模型，对常见的水工结构如水闸和渡槽的热裂缝演变过程进行了模拟，揭示了热应力的分布特征以及裂缝起始和扩展的机制。模拟结果与工程实际观测到的裂缝分布高度一致，验证了模型的可靠性。本研究的成果为水工大体积混凝土结构的裂缝预测和预防提供了理论基础和技术支持。

在工程实践中，混凝土结构的热应力计算方法包括解析方法和数值模拟方法。解析方法通过数学推导或解析表达式来计算早期混凝土的热应力，这些方法通常基于已有的理论模型，将材料参数和结构几何特性代入相应公式进行计算。解析方法的优势在于其相对简单直观的计算过程，使结果易于理解和应用。常见的解析方法包括热应力公式中的变形系数和温度梯度计算方法。然而，对于具有复杂结构特征和边界条件的水工混凝土结构，解析方法的适用性受到一定限制。

数值模拟方法则通过数值计算技术对混凝土结构的热应力进行模拟和分析。常见的数值模拟方法包括有限元分析和离散元方法。这些方法能够更好地考虑混凝土材料的非线性特征和结构的复杂性。它们可以模拟不同形状和尺寸的结构，同时考虑温度场分布和边界条件等影响因素。数值模拟方法的优势在于能够提供更准确和全面的结果，使其适用于分析复杂结构和大规模的热应力问题。因此，目前数值模拟方法被广泛应用。

在早期混凝土热应力计算的理论模型中，研究涵盖了多个方面，包括线性热弹性模型、非线性热弹性模型以及考虑异质性的模型。线性热弹性模型是最常用的模型，用于计算大体积混凝土的热应力。该模型基于线弹性理论，假设混凝土的应变与温度变化之间存在线性关系。模型使用混凝土的热膨胀系数和材料弹性模量来计算温度变化引起的应力。在这一领域，已建立了系统的理论，即由Bofang Zhu提出的大型体积混凝土热应力与温度控制理论。该理论考虑了早期混凝土材料属性随时间的变化，包括蠕变、自生收缩和干燥收缩，并已被广泛应用于工程实践。

尽管该模型在一定范围内能够提供合理的计算结果，但由于其忽略了混凝土的非线性特征和结构的复杂性，可能导致某些情况下计算结果不够准确。为了更精确地描述混凝土的温度行为，研究人员提出了考虑混凝土非线性特性的模型。这些模型可以包括温度对混凝土弹性模量和热膨胀系数的非线性影响，以及温度梯度引起的温度差异。通过考虑这些因素，非线性热弹性模型能够更准确地预测早期混凝土的热应力。

早期混凝土结构通常具有复杂的几何形状和结构特征，因此仅基于均质性假设的模型可能无法满足实际工程需求。一些研究人员考虑了混凝土结构的异质性，如裂缝的影响，并将这些因素纳入模型中进行分析。此类模型可以更好地反映混凝土结构的实际行为，提高计算结果的准确性。

自20世纪90年代以来，关于早期混凝土裂缝的理论模型和数值模拟方法的研究开始出现。最初，大多数研究人员使用塑性模型来描述早期混凝土在裂缝后的性能。随后，随着“弥散裂缝模型”的兴起，研究人员提出了适用于早期混凝土的弥散裂缝模型来表征其性能。随着混凝土损伤力学的发展，经典的损伤模型被扩展到早期混凝土的数值模拟中。同时，一些学者提出了微力学模型，以考虑混凝土的介观结构，特别是粗骨料对早期混凝土抗裂性的影响。

虽然上述研究为深入探索早期混凝土的力学特性奠定了坚实的基础，但这些研究采用了相对简单的裂缝模型，未能充分考虑裂缝演变过程与热传导等过程之间的耦合关系。因此，准确描述早期混凝土的裂缝演变过程以及裂缝后材料力学性能的定量评估仍是需要进一步研究的领域。裂缝传播作为典型的非连续力学问题，可以通过两种主要方法进行数值模拟：基于连续介质的方法和基于离散介质的方法。在有限元框架内，传统的裂缝表示技术如元素删除法和粘聚力区模型具有较高的计算效率。然而，它们存在显著的网格依赖性，并且难以处理复杂的裂缝形态。在非几何描述方法中，扩展有限元方法允许裂缝在元素内传播，但仍需借助辅助技术如水平集方法来追踪裂缝表面，这在处理三维复杂裂缝网络时存在挑战。

相比之下，相场方法引入一个连续的相场变量来描述裂缝，无需预先定义裂缝路径或使用专门的追踪技术。这种方法能够自然地捕捉裂缝的复杂行为，如起始、分支和合并，并在近年来被广泛采用。基于能量最小化原理，相场方法能够直接耦合热传导、机械变形和时间依赖的材料属性，提供了一个统一的多物理场耦合问题的理论框架。

基于脆性相场模型，Nguyen等人提出了水化反应、热传导和裂缝演变的耦合分析方法，有效描述了早期混凝土裂缝演变过程。Li等人进一步发展了一个水热诱导的相场模型，成功模拟了热-孔隙弹性介质中的裂缝传播，为多物理场耦合条件下的断裂演化提供了关键的见解。然而，由于混凝土是一种准脆性材料，脆性相场模型未能合理反映混凝土的损伤和失效行为，而混凝土表现出应变软化特性。在Wu Jianying等人提出的统一相场理论基础上，研究考虑了裂缝过程、水化反应和热传导之间的相互作用。他们建立了裂缝相场演变特征与混凝土水化程度和温度之间的定量关系，提出了一个水化-热-机械耦合的相场粘聚裂缝模型。然而，他们的研究忽略了由于水分蒸发引起的混凝土蠕变变形和自生收缩对材料力学性能的影响。这些因素在考虑混凝土热应力时至关重要，不能被忽视。

本研究针对大体积混凝土热应力分析的理论局限性，提出了一个创新的解决方案。通过将传统的温度场和应力场计算模型与相场断裂理论相结合，建立了考虑蠕变和自生收缩影响的多场耦合模型。模型全面考虑了水化热释放、蠕变变形、自生收缩以及热机械性能演变等关键因素。耦合方程通过增量迭代方法求解，使得混凝土材料的时间依赖行为得到精确描述。

为实现这一目标，本研究采用有限元分析方法进行计算。在空间离散化后，温度场及其梯度、位移场和应变场，以及裂缝相场及其梯度分别表示为：

其中，是温度场节点自由度，是温度场的形状函数，是温度梯度协调矩阵。是位移场节点自由度，是位移场的形状函数，是位移-应变协调矩阵。是裂缝相场节点自由度，是裂缝相场的形状函数，是裂缝相场梯度协调矩阵。

通过将控制方程的弱形式进行求解，可以得到平衡方程。在本研究中，增量迭代算法采用分步求解策略。基于总计算时间和混凝土材料的年龄依赖特性，将计算过程划分为子步骤，其中时间步长根据混凝土年龄动态调整——在早期阶段采用较小的步长，后期逐渐增加。

在每个计算子步骤中，采用迭代求解方法以确保计算结果的准确性。该求解过程可以分为以下四个步骤，以某一增量步骤为例。每个载荷步骤包含四个关键阶段。

第一步：固定裂缝相场值，求解温度场。

第二步：固定温度场和裂缝相场自由度，求解位移场。

第三步：固定温度场和位移场自由度，求解裂缝相场。

第四步：进行收敛性评估。

在上述公式中，代表当前载荷步骤中混凝土的年龄，而代表混凝土年龄的当前载荷步骤。这一处理实际上放大了初始载荷步骤的裂缝历史场变量，确保裂缝不会闭合。

在该步骤中，若裂缝的演变过程满足条件，则迭代终止，输出当前载荷步骤下的温度场、位移场和裂缝相场的计算结果；若条件不满足，则重复前三步直到结果收敛。收敛容忍度可以选取为。

在之前开发的基于MATLAB的相场程序包（PFM）中，已添加了热-机械-相场耦合模型程序。该程序可以通过对应作者获取。

为验证所提出的热-机械-相场耦合模型，采用了基准案例进行测试。如图1a所示，模型底部固定，加载开始时温度场均匀分布为300 K。左、右边界为绝热边界条件。采用自适应网格细化技术进行模拟，初始网格配置如图1b所示。采用AT2相场断裂模型，其热力学参数如表1所示。图1a中的A-A线作为研究裂缝和温度场节点演变的参考线。

图1显示了单边缺口板的拉伸断裂和热传导模型。

表1显示了单边缺口板的热力学参数。

在该研究中，采用的自适应网格细化技术基于本团队之前的工作。主要实施步骤包括根据相场阈值触发局部细化，实现裂缝传播区域的预测性自动细化。使用过渡元素连接细化元素与原始粗网格，确保兼容性并避免悬挂节点。网格重构后，通过映射函数将温度、位移和历史变量转移，以继续计算。这种方法实现了沿裂缝传播路径的自适应细化，从而提高了计算效率和准确性。

如图1a所示，模型底部施加400 K的热载荷，顶部边界采用Dirichlet边界条件（第一类边界条件）维持300 K。在模拟过程中，温度边界条件保持不变。同时，施加了位移载荷，增量步长为1×10^-5 mm。为确保内部温度分布有足够时间，每个载荷步骤设置为100秒。

图2展示了在位移载荷为u=6.11×10^-3 mm、6.18×10^-3 mm和6.24×10^-3 mm时的相场和温度场分布。在裂缝传播开始时，温度场保持均匀分布的稳态。然而，随着裂缝向右传播，温度场逐渐发生变化。当裂缝传播接近完成时，温度场被明显分为上下两部分。上部分稳定在300 K，下部分稳定在400 K。这种现象源于热传导降解函数的定义。当结构完全断裂时，降解函数w(φ)达到零，表示裂缝不导热，从而隔离上下区域的温度传递，导致观察到的温度分割。

图2显示了不同载荷步骤下的裂缝传播和温度场分布。

图3展示了A-A线节点在三个载荷步骤中的温度和相场值演变。如图3所示，在裂缝起始前，温度场保持均匀连续。一旦裂缝开始传播，温度场发生变化，当裂缝穿过A-A线时，由于裂缝的热绝缘作用，温度场出现急剧转变。当裂缝传播接近完成时，上下区域表现出不同的温度值，裂缝区域的温度突然从400 K变为300 K。这表明损伤和断裂的发展改变了材料的热导率，影响了温度场的演变。这些结果验证了所提出的模型在处理热-机械耦合条件下的断裂问题方面的有效性。

为验证所提出的热-机械-相场耦合模型，采用了一个具有方形孔的早期混凝土圆盘进行实验。如图4a所示，圆盘的半径为0.5米，包含一个边长为0.4米的方形孔。为了模拟，圆盘的厚度被设置为0.05米。模型模拟了混凝土圆盘在浇筑后由于水化热引起的热裂缝。利用模型的对称性，选择了一个四分之一模型进行模拟。初始网格配置如图4b所示。在模型的左侧施加了x方向的位移约束，在底部施加了y方向的位移约束。

混凝土圆盘的浇筑温度为20°C，外部环境温度也为20°C。圆盘的外边界受到对流边界条件（第三类边界条件）的影响，而其余边界则被视为绝热边界。采用AT2模型进行模拟，其热力学参数参考了先前的研究，如表2和表3所示。

表2和表3分别显示了早期混凝土圆盘的固定热力学参数和随年龄变化的热力学参数。

图4显示了早期混凝土圆盘模型和网格。

图5展示了在不同时间间隔下，具有方形孔的混凝土圆盘的相场和温度场分布。如图5a所示，在浇筑后的前4.5小时内，由于水化反应释放的热量，圆盘的温度持续上升，中心点温度达到65.7°C。此时，圆盘的拉应力相对较低，除了与地板连接处，拉应力较高，超过混凝土在该年龄的抗拉强度，导致局部损伤。随后，结构的温度逐渐下降，圆盘内部的压应力逐渐减少，导致在4天后出现拉应力。随着圆盘温度的进一步下降，拉应力逐渐增加，在6天时达到拉应力的极限值，导致圆盘内部多个区域逐渐出现损伤。此时，结构底部的拉应力相对较高，而上部的拉应力较低，受损区域的应力相应降低。

图5显示了不同时间下的裂缝和温度场分布。

图6显示了早期混凝土圆盘的裂缝和温度场分布。

在整个温度诱导裂缝起始和传播过程中，有两个值得注意的点。首先，在拉应力达到峰值之前，相场值没有显著变化；然而，一旦应力达到峰值，它就不再增加。随后，在相应位置出现损伤，相场值增加，应力值减少，如图10b和图10c所示。其次，两条主要裂缝对称分布，将结构分为三等分，这符合裂缝传播路径使系统总能量最小化的原理。如果结构温度继续下降，新的裂缝仍会将完整区域分为三等分，如图9c所示。值得注意的是，两条主要裂缝对称分布，并精确地将结构分为三等分，这符合裂缝传播路径使系统总能量最小化的原理。如果结构温度继续下降，新的裂缝仍会将完整区域分为三等分。

为验证时间依赖参数的影响，建立了三个不同的计算案例，考虑蠕变和自生收缩对应力和相场值的影响。案例1考虑了蠕变和自生体积收缩，案例2基于案例1忽略了蠕变的影响，案例3基于案例1忽略了自生体积收缩的影响。关键点在浇筑后的温度、应力和相场值的时间历史曲线如图10所示。

从图中可以看出，所有三个案例都显示了相同的温度发展曲线，而蠕变和自生收缩仅影响应力场和相场值。在温度下降阶段，混凝土的内部应力迅速增加，相场值在拉应力达到峰值前没有显著变化。拉应力达到峰值后，它不再增加，损伤出现在相应位置，相场值增加，应力值减少，如图10b和图10c所示。比较三个案例中应力和相场值的变化，可以看出，当忽略蠕变时，拉应力在早期冷却阶段增长更快。因此，蠕变效应可以降低拉应力水平，从而降低结构裂缝的风险。因此，在实际工程应用中，当考虑蠕变时，可以适当放松温度控制措施。当忽略自生收缩时，拉应力在早期冷却阶段增加稍慢，表明自生收缩导致拉应力增加并提高了结构裂缝的风险。

值得注意的是，对于常规的水利工程，蠕变参数通常使用Bofang Zhu推荐的公式（公式18）进行计算，而无需进行专门的蠕变实验。自生收缩参数则基于实验或经验数据使用双指数函数进行拟合。当自生收缩的影响显著时，工程实践中通常会添加膨胀剂来抵消收缩效应。因此，本文选择了这三个典型的工况，未进一步对特定参数进行敏感性分析。

此外，图12展示了渡槽的模拟温度裂缝传播，其中在墙墩上可以观察到四条温度裂缝。与图13相比，模拟的裂缝分布与实际工程中观察到的温度裂缝高度一致。需要指出的是，模拟的裂缝宽度比实际物理裂缝更大，这是由于长度尺度参数（l）控制了裂缝扩散区的范围。这使得数值模拟的裂缝显示为具有一定宽度的损伤带，反映了实际结构中损伤区域的分布范围，而不是物理裂缝的几何尺寸。重要的是，这种表示方式不会影响裂缝传播路径或裂缝倾向的评估。

渡槽温度裂缝的原因是，当墙墩被浇筑时，其下方的地板已经相对成熟，对上部施加了较强的约束。在早期阶段，墙墩的混凝土强度相对较低，随着结构温度的升高，内部区域产生较小的压应力。由于结构相对较薄，热传导较快，当水化热产生的热量小于结构散热时，墙墩的温度会下降。此时，混凝土的强度逐渐增加。随着温度的下降，结构内部的压应力也逐渐减少，导致拉应力的出现。当拉应力增加并超过混凝土的抗拉强度时，墙墩将经历损伤，逐渐形成温度裂缝。

研究提供了一个可靠的理论基础和实用的数值模拟工具，用于早期混凝土的热裂缝预测和预防，为工程应用做出了贡献。