基于近场动力学的脆性材料热冲击损伤研究

廖 鑫 邓 强 张 超

（成都飞机工业（集团）有限责任公司，四川 成都 610092）

脆性材料（例如陶瓷、岩石等）在日常生活中随处可见，在温度载荷作用下，产生的内应力会导致断裂破坏，如何准确模拟脆性材料的复杂裂纹萌生和扩展一直是国内外研究的难题。

PD 理论假设连续体中粒子之间存在相互作用，不需要指定额外的失效准则，材料的断裂损伤自然发生，常用于研究裂纹扩展问题。Silling S A 等[1]根据PD 理论成功地捕捉了带有初始缺陷厚板的裂纹萌生和扩展现象，证明了近场动力学可以用来模拟裂纹扩展问题。秦洪远等[2]构建了可以消除泊松比限制的双参数微弹脆性本构模型，模拟了含双裂纹巴西圆盘的劈裂破坏过程。孙杰等[3]同样使用双参数微弹脆性本构模型研究了初始裂纹角度对多裂纹扩展路径和临界破坏载荷的影响。Huang D 等[4]分析了混凝土结构的静态弹性变形和开裂行为。针对脆性材料的热冲击损伤，郑州大学的王振宇[5]采用PD 方法研究了岩石材料的动态断裂与瞬时热传导行为。Zhang H 等[6]提出了一种考虑热力载荷的常规态基PD 模型，以预测金属和陶瓷双材料结构的断裂损伤问题。PD 理论适合处理复杂的裂纹萌生和扩展问题，根据该理论可以较好地研究脆性材料的热冲击损伤问题。

该文基于PD 理论建立一种分析脆性材料热冲击损伤问题的方法，推导近场动力学的传热方程和热力耦合本构方程，模拟脆性材料在热冲击载荷作用下发生多裂纹随机萌生和扩展的现象。

1 键基PD 热力耦合模型

键基PD 的运动方程由虚功原理推导而来，满足动量守恒定律和角动量守恒定律，如公式（1）所示。

式中：ρ为质点的密度；ü为加速度矢量；H为视界范围；b为体积力密度矢量；f为键力密度矢量，与质点的材料属性、质点间的拉伸情况和质点的视界半径有统计学意义；x为质点的坐标矢量；t为时间；u为质点的位移矢量；u'为x视界内其他质点的位移矢量；x'为视界内其他质点的坐标矢量。

键力密度矢量f与键的伸长率s呈线性关系，如公式（2）所示。

式中：ξ、η分别为2 个质点之间的相对位置、相对位移；c为键合常数；s为键的伸长率；μ为质点的位移矢量。

对各向同性材料来说，可以通过WPD=WCM（WPD为基于近场动力学的应变能；WCM为基于连续介质力学的应变能）来获得PD 参数与连续介质力学参数间的关系[3]，键合常数c在平面应变状态下如公式（3）所示。

式中：h为厚度；δ为近场半径；E为弹性模量。

引入一个标量值函数g，以控制质点间键力的存在和消失，如公式（4）所示。

式中：sc为键失效时的临界拉伸；α为热膨胀系数；T为温差。

令消除所有穿过新裂纹表面的键所做的功等于临界能量释放率G，可以推导二维情况临界拉伸与临界能量释放率间的关系，如公式（5）所示。

式中：Gc为临界能力释放率。

温度是标量，但是热传导具有方向性，质点之间的温差是热传导的驱动力，由傅里叶热传导定律可知，热流密度q如公式（6）所示。

式中：k为材料的热导率；T为温度梯度。

键基PD 传热方程可以用温差、热流密度等参数来表示，如公式（7）所示。

式中：ρ为质点的密度；cv为质点的比热容；T为 温度对时间的一阶导数；T'为质点x'的温度；x'为x视界内其他质点的坐标矢量；Vx'为x视界内所有质点；sb为热源；fh为热响应函数。

热响应函数与质点的温差、热导率和质点的相对位置有统计学意义，因此热响应函数如公式（8）所示。

式中：τ为质点的温差；κ为微热导率，可以通过假设简单线性温度场下PD 热势能和经典连续介质力学热势能是等价的来获得。

二维情况下近场动力学的微热导率如公式（9）所示。

式中：s为厚度。

对流边界条件下的单位体积产热率如公式（10）所示。

式中：n为垂直于边界Sf的法线方向单位矢量；Sf为实体热流边界。

包括标量函数g的传热方程如公式（11）所示。

式中：hs为单位体积的产热率。

通过时间的积分可以获得所有质点的当前温度，然后在键力密度矢量中引入热膨胀项c'T，以表征温度变化与键力间的关系，并计算键力密度矢量，如公式（12）所示。

式中：c'为PD 热力耦合参数。

c'与键合参数c的关系如公式（13）所示。

式中：α为热膨胀系数。

获得特定温度场影响下的所有质点键力密度矢量后，将其代入键基PD 运动方程，如公式（14）所示。

该方程可计算所有质点的加速度ü，再积分求得所有质点的位移结果。

2 数值模拟

2.1 PD 模型设置

PD模型取直径为l3 mm、厚度为1 mm的圆形薄板建立二维模型，薄板初始温度为300 ℃~500 ℃，温度分布均匀，并在圆盘外周施加1 圈对流换热边界，环境温度设置为T0=15 ℃。材料参数如下：E=370 GPa，υ=1/3，p=3 980 kg/m3，Gc=12.16 J·m2。

初步设置模型的质点间距Δ=50 μm，视界尺寸为δ=3.015Δ。热冲击温度为500℃，温度场计算步长设置为dt0=1×10-5s，位移场计算步长为dt=5×10-3s，为了保证ADR收敛，设置最大收敛步数为200、收敛安全系数为5。

2.2 PD 参数

进行δ收敛分析，以确定视界最优值，在较细的网格尺寸（Δ=50 μm）下进行，δ如公式（15）所示。

式中：δ为视界尺寸；m为视界尺寸与网格大小的比值，m分别取2、3、5、7、10 和15。

保持其他条件完全一致，给出水平中线上的质点，在不同m值下的PD 解与FEM 解的相对误差欧几里德范数的均值如公式（16）所示。

式中：ai为质点的待分析参数值；为对照组的参数值；h为质点数量。

由表1 可知，当m=3 时，温度场和位移场的相对误差最小分别为0.052 3%和0.547 5%，计算精度较高。

表1 不同m 值的PD 解相对误差

分析网格大小的最优值，将m值固定为3，以保证结果精度，令Δ为13 μm、25 μm、50 μm、100 μm、200 μm 和500 μm，其他条件完全一致，不同Δ值的PD 解与FEM 解的相对误差欧几里德范数的均值见表2。当Δ分别为50 μm和100 μm 时，相对误差较小，分别为0.547 5%和0.156 9%。考虑计算精度和效率，认为Δ为50 μm 或100 μm 是较好的选择。

表2 不同Δ 值下PD 解与FEM 解的相对误差

综上所述，当T=500 ℃时，m=3、Δ=50 μm 的结果精度较高。

3 计算结果分析

采用近场动力学方法计算不同温度下的热冲击裂纹扩展情况，在裂纹数量和长度方面，与试验结果[7]进行对比，如图1 所示。在热冲击载荷的作用下，陶瓷薄板发生了垂直于周向的裂纹萌生和扩展现象，裂纹具有周期性和层次性。热冲击裂纹沿陶瓷薄板四周均匀分布，周期性地出现长裂纹、中裂纹和短裂纹，裂纹层次鲜明，均由圆板四周萌生并向圆心处扩展。当温度为300 ℃时，PD 模拟结果出现4 条长裂纹，分布在0°、90°、180°和270°方向，分布较均匀，而试验结果同样明显出现4 条长裂纹，数量一致，大致分布在0°、70°、190°和270°方向。当温度为350 ℃时，PD 模拟结果出现8 条长裂纹，分布在0°、30°、60° 、90°、120°、150°、180°和270°方向，而试验结果明显出现6 条长裂纹，分布在0°、60° 、120°、150°、210°和270°方向，相比之下，PD模拟结果多产生了2 条长裂纹，二者有5 条长裂纹方向一致。当温度为400 ℃时，PD 模拟结果出现8 条长裂纹，分布在0°、30°、60° 、90°、120°、150°、180°和270°方向，而试验结果明显出现6 条长裂纹，分布在30°、60° 、120°、190°、270°和330°，相比之下，PD 模拟结果多产生了2 条长裂纹，二者有4 条长裂纹方向一致。当温度为500 ℃时，PD 模拟结果出现8 条长裂纹，分布在0°、30°、60°、90°、180°、210°、240°和270°方向，而试验结果同样明显出现8 条长裂纹，数量一致，分布在30°、60° 、140°、190°、210°、240°、270°和340°方向，二者有5 条长裂纹方向一致。综上所述，PD 模拟结果与试验结果的热冲击裂纹趋势较一致，PD 方法成功捕捉到了垂直于圆形薄板周向的裂纹萌生和扩展，裂纹的周期性和层次性规律也十分相似，主裂纹的长度也比较接近。但是裂纹分布情况存在一定差别，推测可能是由试验所用陶瓷板材料分布不均匀或存在初始缺陷导致，试验误差引起的淬火不均匀性也会影响试验结果。

图1 热冲击裂纹对比（左为PD 模拟结果，右为试验结果）

图2 统计了300 ℃~500 ℃，PD 模拟结果与试验结果的裂纹数量对比，根据裂纹长度与圆盘半径R的比值归类，将裂纹长度在0.1R~0.5R（一种定量的对比手段，将裂纹长度与圆盘半径进行对比，当裂纹长度等于1R时，表示裂纹贯穿；当裂纹长度0.5R~1.0R时，表示裂纹较长接近贯穿；当裂纹长度0.1R~0.5R时，表示裂纹较短；当裂纹长度小于0.1R时，属于微损伤，裂纹小且数量庞大，统计难度较大，因此忽略）范围内视为短裂纹，大于0.5R为长裂纹，忽略小于0. 1R的微小裂纹。模拟结果的长裂纹和短裂纹交错出现，规律性较强。当热冲击温度为300 ℃时，PD 模拟结果和试验结果的长裂纹数量均为4 条，短裂纹数量均为20 条，因此总裂纹数量相等。当热冲击温度为350 ℃时，PD 模拟结果和试验结果的长裂纹数量相差2 条，短裂纹数量均为20条，总裂纹数量较接近。同样，当热冲击温度为400 ℃时，PD 模拟结果和试验结果的长裂纹数量相差2 条，短裂纹数量均为23 条，总裂纹数量较接近。当热冲击温度为500 ℃时，PD 模拟结果和试验结果的长裂纹数量均为8 条，短裂纹数量均为28 条，因此总裂纹数量相等。综上所述，在不同热冲击温度下，对陶瓷薄板热冲击裂纹进行模拟，PD 方法较准确，不仅短裂纹数量高度一致，而且长裂纹数量偏差较小。

图2 试验与PD 模拟裂纹数量

最后，对比不同温度下的PD 模拟结果发现，长裂纹数量由300 ℃的4 条逐渐增至500 ℃的8 条，而短裂纹数量由300 ℃的20 条逐渐增至500 ℃的28 条，随着热冲击温度的升高，长裂纹和短裂纹的数量逐渐增加，主裂纹长度也逐渐增加，陶瓷板的损伤程度加剧。这一现象说明，热冲击载荷不仅会影响裂纹的数量，而且还会影响裂纹的整体长度。

综上所述，该文的PD 方法较好地模拟了陶瓷圆形薄板的热冲击裂纹扩展问题，仅使用少量计算资源就得到与试验较一致的现象，证明了PD 方法对该问题有较高的适用性和经济性。

4 结语

该文基于键基PD 理论建立了热冲击损伤模型，并对陶瓷薄板的热冲击过程进行模拟以及参数收敛性分析，保证了计算精度和效率，最后将模拟结果与试验结果进行比较。

模型的视界尺寸δ和网格尺寸Δ对计算精度和效率的影响较大。将m值分别取为2、3、5、7、10 和15，以改变视界尺寸。当m=3 时，PD 方法与FEM 方法的温度和位移结果相对误差均较小，计算时间短。网格尺寸Δ值分别为13 μm、25 μm、50 μm、100 μm、200 μm 和500 μm，当Δ为50 μm或100 μm 时，相对误差较小。因此认为m=3、Δ=50 μm 是该模型精度高、效率高的组合。

当热冲击温度为300 ℃~500 ℃时，分别对比了PD 结果与试验结果的裂纹情况，发现无论是裂纹分布规律，还是裂纹数量、长度和形式都较为一致，现象吻合较好，说明PD 模型较好地模拟了陶瓷薄板的热冲击损伤问题。此外，还发现随着温度升高，长裂纹和短裂纹的数量均逐渐增加，长裂纹在较高的热冲击温度下变得更长，陶瓷板的损伤程度加剧。