独立变量的热力学约束：多孔弹塑性断裂相场模型流固耦合项的重构

一、引言

相场断裂建模通过引入连续标量损伤场将尖锐裂纹面规则化为有限宽度，从而摆脱了传统有限元方法中显式裂纹面追踪与扩展判据的束缚，使Griffith断裂问题的变分表述得以与各类标准及新兴有限元方法相衔接。经过多年发展，该框架已从脆性固体相继拓展至内聚力断裂、动态断裂、各向异性扩展、疲劳循环加载及多物理场耦合等方向，并被广泛植入开源软件生态。然而，在具有分辨微结构的非弹性材料断裂模拟方面，研究大多局限于二维表示，全三维应用因捕捉微结构特征与裂纹两侧陡峭损伤梯度所需的极细网格而计算代价高昂，长期难以推进。

该文指出，当采用线性形函数时，需要足够多的单元横跨裂纹半宽，即单元尺寸须远小于相场长度尺度；而高阶单元可在单元尺寸与长度尺度相当的条件下给出精确解且避免线性单元的虚假振荡。无矩阵方法避免了稀疏全局矩阵的显式组装与存储，直接在单元与积分点层面施加算子，提高算术强度、降低访存流量，与GPU架构天然契合，从而使高阶单元在每自由度意义上比低阶更高效。将无矩阵表示与多重网格预条件子结合，已在超弹性与脆性相场断裂中获得近优扩展性。

Jeremy Thompson、Jed Brown、Fabio Di Gioacchino等人在“Computational Physics"期刊发表了文章“Matrix-free phase-field modeling of fracture in micromechanical testing simulations of inelastic materials"，该文的贡献有二：其一，将相场断裂与Perić-Dettmer非弹性本构框架相结合，提出串联组装的流变断裂单元，使损伤累积在不改变非弹性性质的前提下使整个本构块失活，且损伤仅由储存于本构块中的纯弹性应变能驱动，从而实现微结构尺度下损伤区与微孔洞、微裂纹位置的对应；其二，将该框架实现于原生支持无矩阵算子与GPU加速的开源固体力学库Ratel中，采用p-多重网格预条件，并在El Capitan高性能计算原型机上完成算例验证。据作者所知，这是包含有限应变非弹性并使用GPU的无矩阵相场断裂研究。

二、Perić-Dettmer流变本构框架

该方法首先从长度尺度的视角厘清塑性区与损伤区的关系（图1）。在均匀化长度尺度上，损伤区代表含微孔洞与微裂纹的区域；而在微结构尺度上，损伤区与微结构中微孔洞形核、微裂纹扩展的位置直接重合，损伤的区域不再能够累积非弹性变形。这一区分在流变学表示上体现为：将断裂赋予与非弹性过程控制单元相串联的专用流变单元。

图1  不同长度尺度下塑性区与损伤区的示意区分：(a)均匀化尺度；(b)微结构尺度

Perić-Dettmer框架通过Hooke（弹性）、Maxwell（粘弹性）与Prandtl（弹塑性）流变分支的并联组装构造各向同性材料的非弹性本构模型（图2）。该文在其基础上串联一个流变断裂单元，其中储存于各分支中的弹性能通过退化函数打折后共同贡献于断裂驱动力。损伤累积使整个Perić-Dettmer块失活而不影响其内在非弹性性质，串联组装还保证非弹性与断裂性质仅在均匀化长度尺度上相互影响。

图2  非弹性与断裂力学的流变单元组装：Perić-Dettmer块中Hooke(H)、Maxwell(M)、Prandtl(P)分支并联，断裂单元(F)串联

该框架专用于Hencky材料（图3）：采用乘法分解将总变形梯度分为非弹性与弹性部分，通过Hencky应变能泛函与指数映射时间积分，当前构型中的应力更新可方便地复用小应变回映射算法，仅需在欧氏应变空间与对数应变空间之间增加前、后处理映射；再通过在基于主应力的格式中将更新限制于偏量分量，即可用极少的材料参数高效实现涵盖广泛各向同性力学响应的模型。乘法基还保证了非共轴载荷下的物理响应，这是大应变加法形式所不具备的。

图3  Perić-Dettmer流变块的总变形梯度乘法分解

三、断裂相场与流变断裂单元

断裂单元的能量泛函采用Ambrosio-Tortorelli形式的相场断裂模型。AT2模型中损伤自加载伊始即开始累积；AT1模型则引入能量阈值，损伤仅在弹性能密度超过阈值后才开始积累，适合裂纹形核发生在变形后期的情形。该文根据各算例的物理特征选择相应模型：脆性响应采用AT2，孔洞形核位于变形后期的颗粒增强材料拉伸算例则采用AT1。作用在断裂单元上的驱动力取为各流变分支弹性应变能贡献之和；为处理压缩应力状态下的裂纹闭合与摩擦等行为，能量分解采用仅由正（拉伸）特征应变能驱动损伤的谱分解格式。

该方案与既有将非弹性应变能项纳入损伤驱动力或使断裂韧性随非弹性应变累积而退化的方案形成对照：在流变串联组装下，纯弹性应变能是损伤的驱动力，损伤不通过屈服面收缩诱发塑性软化，而是在裂纹形成的区域抑制进一步的塑性累积，导致承载能力的快速丧失而非渐进颈缩至零载荷——这与实验中延性材料拉伸测试以灾难性断裂与回跳响应告终的特征相符。

数值稳定性方面，该文引入残余刚度因子与损伤黏性两项正则化：前者保留少量未损伤刚度以避免失强，后者对损伤演化施以黏性松弛，共同防止裂纹局域化时牛顿迭代收敛性的丧失，并使准静态框架能够数值持续地再现快速卸载乃至回跳行为。

四、无矩阵数值实现与计算平台

该框架实现于基于libCEED与PETSc构建的Ratel开源固体力学库。无矩阵实现通过单元局部算子在飞行中评估耦合残差与雅可比的作用，避免显式组装：全局自由度经聚集算子映射至单元局部自由度，基插值与母单元梯度算子在积分点上插值出场值及其梯度，残差与雅可比向量积复用同一套限制、基、积分与转置操作序列。四分量解向量包含三个位移分量与损伤场，通过引入与弹性模量、长度尺度及断裂能相关的相场残差缩放因子，使损伤残差幅值与力学残差可比，降低耦合雅可比的病态程度。

非线性问题采用带回溯线搜索的牛顿-拉夫逊法与自适应时间步进求解；每次牛顿迭代的线性化系统用GMRES求解，预条件采用以代数多重网格粗层求解器与Chebyshev光滑的p-多重网格。由于特征长度尺度与网格尺寸相近，损伤场被排除在提供给AMG的近零空间之外，降低了网格复杂度。模拟均在劳伦斯利弗莫尔国家实验室Tioga集群（El Capitan原型系统，AMD架构）上运行，经PETSc与Kokkos以HIP后端编译执行于MI250X GPU；微结构分辨算例的问题规模达约2000万自由度，在3节点、24逻辑GPU与12小时分配窗内约500个时间步完成。该文还给出微结构生成流程：借助microstructpy控制颗粒形状、尺寸与体积分数，经trimesh转换STL并过滤重叠后由Gmsh剖分、Netgen优化网格质量。

五、数值算例

5.1  尖锐缺口板剪切试验

第一个算例为尖锐缺口板的I-II混合型剪切（图4a为其流变配置）：底面夹持、顶面施加侧向位移的平面应变边界条件被修改为简支以允许裂纹在边界张开，剪切变形由两侧面的相对位移施加。预测的裂纹自缺口处形核，形核时刻与力-位移曲线上特征凸峰的第一个峰值重合；随后裂纹以与板边成一定角度在拉伸区内扩展，最终弯折向下并以直角抵达底面（图5）。与文献中夹持边界下裂纹贴底面传播的路径不同，简支边界下的主裂纹路径明显下弯。

图4  本研究采用的流变单元配置：(a)缺口钢板（算例1）及算例2、3的颗粒相；(b)粘弹性基体（算例2）；(c)粘弹塑性基体（算例3）

图5  尖锐缺口板剪切试验：(a)有限元网格与边界条件；(b)裂纹路径灰度图；(c)力-位移曲线及收敛率报告的阶段标注

非线性求解器在损伤累积与缺口裂纹形核阶段保持二次收敛，裂纹扩展阶段收敛率有所退化，最终韧带阶段则先经历缓慢的残差下降再进入局部收敛域，对应裂纹趋于失稳扩展的状态。裂尖应力演化表明切断最终韧带的回跳倾向：扩展阶段裂尖应力水平相当，而到最终韧带时约增大一倍（图6）。残余刚度因子与损伤黏性的引入有效维持了求解器收敛。

图6  裂纹扩展三个连续加载阶段的σxx演化（损伤值高于0.95的区域被消除以显示裂纹张开）

5.2  粘弹性颗粒材料的压缩

第二个算例模拟含较高体积分数硬椭球颗粒的粘弹性基体圆柱的准静态单轴压缩，代表聚合物粘结颗粒材料的微力学测试。1498个随机分布的椭球颗粒平均尺寸0.02 mm，约占柱体体积的0.2（图7）。颗粒相用脆性断裂流变模型描述，基体相用粘弹性模型描述，两相均采用AT2模型；借助Ratel已实现的接触算法考察润滑与侧向约束两类加载条件。

图7  压缩下的粘弹性颗粒材料：(a)颗粒俯视渲染；(b)侧视渲染；(c)试样整体有限元网格

三种基体模型（纯弹性E、有限偏量黏性VE、附加损伤黏性DV）在两种约束条件下的力-位移响应表明：提高偏量黏性降低流变应力，提高损伤黏性增大峰值力并使峰后软化更平缓；两类约束条件下峰值力与全断裂时的位移相当，但润滑条件下峰值力在更大应变处出现且峰后载荷下降更陡（图8）。

图8  润滑与侧向约束条件下三种基体模型的压缩力-位移响应：(a)侧向约束；(b)润滑压板

损伤演化与裂纹扩展模式解释了上述差异。侧向约束时，压板下方形成穹顶状无损区，损伤优先在偏应力集中的试样边缘形核，同时试样中心出现形核点并随变形增加聚合成贯通的长裂纹网络，最终中心与边缘裂纹汇合，三维上呈现脆性颗粒复合材料压缩的X形断裂模式（图9、图11a）。提高损伤黏性使损伤累积扩散化，侧向区域碎化程度更高，与PBX实验观察更为接近；而粘弹性仅延缓损伤累积、基本不改变裂纹模式，印证压缩下损伤由偏应力主导。润滑压板消除了边缘偏应力集中与裂纹起始，损伤分布更均匀，不形成X形模式，而是出现局限于试样下半部的主导裂纹（图10、图11b），断裂面范围较小、消耗断裂能较少，对应润滑条件下更快的峰后载荷跌落。

图9  侧向约束条件下试样中央截面上相场损伤随变形的演化：(a)峰值载荷I；(b)峰后75%峰值载荷II；(c)50%峰值载荷III

图10  润滑压板条件下相场损伤随变形的演化：(a)峰值载荷I；(b)75%峰值载荷II；(c)50%峰值载荷III

图11  压缩下的三维断裂形貌：(a)侧向约束（X形断裂模式）；(b)润滑压板

5.3  颗粒增强粘弹塑性材料的拉伸

第三个算例为低体积分数硬球形颗粒增强粘弹塑性基体的单轴拉伸，对应颗粒增强合金的微力学测试。约390个直径15 μm的颗粒置于平均网格尺寸5 μm的狗骨试样中；通过改变随机数生成种子获得三个统计等效微结构M1-M3，保持力学参数不变以单独考察微结构变异性的影响。两相均采用AT1模型，颗粒相为脆性断裂模型，基体为粘弹塑性模型。

图12  拉伸下颗粒增强粘弹塑性材料：(a)狗骨试样网格与三个统计等效微结构M1-M3；(b)预测力-位移曲线（虚线为抑制损伤的参考响应，I-IV标注峰值与峰后阶段）

抑制损伤的参考响应中，Hooke与Maxwell分支的应变硬化与颗粒强化共同阻止颈缩、维持载荷上升；损伤激活后，曲线平滑趋近峰值力并呈渐进峰后软化，且试样发生颈缩——损伤促进几何失稳，放大了微结构变异性对峰后响应的影响，体现为三条曲线在该阶段的轻微分叉（图12）。

图13  峰后加载中三个微结构的演化：(a)累积塑性应变；(b)相场损伤（列I-IV对应图12b标注的变形阶段）

截面场演化再现了延性断裂的经典序列（图13）：加载中塑性应变在硬颗粒周围累积形成特征性变形羽流；峰值应力时损伤在相同区域局域化并优先在相邻颗粒间的基体通道内发展——由于大部分偏应变能已通过塑性流动耗散，损伤主要由体积应变能驱动，与孔洞形核及长大机制一致；颈缩区内孔洞开始聚合，形成垂直于加载轴的盘状裂纹；随应变增加，侧向韧带中发展出与加载方向成45°的剪切带，裂纹沿剪切带向外扩展，完成从I型断裂到剪切主导II型断裂的转变，对应常规延性材料拉伸的杯锥断口。裂纹还优先沿颗粒-基体界面形核扩展、使颗粒裸露，与断口分析中的常见观察一致。由于损伤不诱发塑性软化，裂纹形成的区域不再出现盘状塑性局域化图案，承载能力快速丧失，再现了灾难性断裂与回跳收尾的不稳定失效（图14）。

图14  最终施加应变下三个微结构的累积塑性应变截面图（损伤值高于0.95的区域被遮蔽，虚线标出M1的杯锥断裂轮廓）

六、结论与展望

该文将Perić-Dettmer非弹性本构与相场断裂相结合，提出串联组装的流变断裂单元：弹性应变能作为损伤的驱动力，损伤使整个本构块失活而不改变其非弹性性质。算例表明该框架能够再现剪切带、表面断裂特征与定义延性断裂的裂纹模式，支持了流变框架特别适合微结构分辨模拟、损伤形核演化可与孔洞形核及微裂纹扩展相对应这一初始假设。

在计算层面，无矩阵实现借助Ratel库在El Capitan原型GPU系统上稳健求解高度非线性问题：位移-损伤弱耦合、高阶离散与损伤黏性的正则化共同保证了含p-多重网格预条件的四分量单块格式的稳健收敛，问题规模达两千万自由度。该实现已支撑了粘弹性断裂替代模型研发所需的数百次高保真模拟，展现出服务不确定性量化等大规模研究流程的潜力。

展望方面，作者计划将框架拓展到基于晶体塑性有限元模型的各向异性非弹性：塑性流动将取决于晶粒的Schmid因子而变得空间非均匀，弹性-塑性转变较各向同性情形更平缓，进一步佐证理想塑性Prandtl分支的合理性；重点将放在滑移传递与晶格曲率介导等塑性应变传递机制对抑制晶界等微结构界面处损伤形核的作用，并将微结构分辨模拟与微柱压缩等微力学测试实验相对接。