船用灰铸铁材料裂纹扩展的近场动力学研究

李 双，吕海宁，黄小华，杨建民

（1.上海交通大学海洋工程国家重点实验室，上海 200240；2.上海交通大学三亚崖州湾深海科技研究院，海南 三亚 572000；3.广西大学防灾与结构安全教育部重点实验室，南宁 530004）

0 引 言

在船舶与海洋结构物的设计制造及使用过程中，不可避免地会出现裂纹的萌生和扩展问题，这对结构的承载能力和安全性能起着明显的削弱作用[1]。综合性能优良的高强度灰铸铁作为广泛应用于船舶工程各种复杂机械设备构件（内燃机的机架、汽缸盖、燃油泵、船体水泵、高压箱体、基座、推进器、轴承、高速齿轮、飞轮、防护外壳等）的主要金属材料[2]，一旦产生裂纹扩展和脆性断裂问题，将导致船体发动机无法正常运行、船舶动力装置受损，带来安全隐患并影响船舶的日常工作。因此，对船用灰铸铁构件进行裂纹扩展行为模拟，研究试件从裂纹萌生、扩展、交汇直至断裂的整个服役过程，对掌握裂纹扩展规律、预防构件脆性断裂、评估船体灰铸铁结构的完整性和安全性具有重要的指导意义。

目前，用来模拟裂纹扩展的数值方法中，有限元法、内聚力法和扩展有限元法等得到广泛应用。余天堂等[3]采用改进的扩展有限元法对张开型和闭合型裂纹的裂纹扩展进行了模拟；唐志波等[4]基于船体表面微裂缝在应力和海水腐蚀耦合作用下的扩展机理，根据Charles-Hillig 模型的有限元方法模拟了裂缝的动态扩展过程；陈家旺等[5]基于有限元分析软件ANSYS 研究了带有初始表面裂纹的潜艇典型结构在不同下潜深度状态下的裂纹扩展情况；钱怡等[6]提出一种基于累积损伤理论和裂纹尖端弹塑性应力场的数值方法来预测裂纹的扩展行为；王敏等[7]用扩展有限元方法来分析含裂纹板的裂纹扩展，采用求解Reissner 板问题的混合插值单元，避免了刚度矩阵无穷大导致的单元刚度矩阵变为奇异的问题；苗婷等[8]采用内聚力方法对X65管线钢进行了断裂韧度评定及裂纹扩展行为研究；李亚政等[9]基于损伤力学和扩展有限元方法对舰船蒸汽轮机叶片进行了裂纹扩展研究；李卓等[10]采用内聚力和扩展有限元方法对含双边预制裂纹的灰铸铁金属材料进行了裂纹扩展的试验和模拟研究；Chang等[11]采用扩展有限元方法对含单边裂纹的灰铸铁材料进行了裂纹扩展和断裂分析。

上述基于连续介质力学的数值模拟方法，例如有限元虽然在结构变形分析中取得较大成功，但在解决破坏问题时，由于运动方程涉及位移对空间的微分，在不连续处微分型运动方程没有定义，出现奇异性，必须人为地对不连续处进行专门的处理；而且上述有限元法、内聚力法和扩展有限元法等在处理不连续问题时，都不可避免地存在网格依赖性、裂纹的萌生和扩展需要引入额外的断裂准则、涉及多裂纹相互作用时存在较大局限性等问题，这无疑增加了数值模拟的复杂性。目前，针对船用灰铸铁的裂纹扩展研究还较少，灰铸铁作为一种脆性金属材料，具有良好的铸造、减振、耐磨、切削加工性能和较低的缺口敏感性，被广泛应用于结构较为复杂的船海工程中，因此，有必要研究新的数值计算方法来对船用灰铸铁的裂纹扩展行为做进一步分析。

2000 年，新兴的近场动力学（Peridynamics，PD）理论[12]的提出弥补了这些不足。它采用位移的空间积分形式完全重构了传统连续介质力学理论，这种积分型方程在不连续处仍有定义，避免了奇异性问题，且不存在网格依赖性。PD 的断裂准则通过一个记忆函数嵌入本构方程中，能够模拟裂纹自发的萌生和扩展，实现了对连续介质和非连续介质力学行为的统一描述，是一种特别适合模拟材料破坏的力学理论[13]。Silling[14]采用PD 方法研究了冲击荷载作用下Kalthoff-Winkler 韧性钢试件的裂纹扩展和动态断裂行为；黄小华等[13]采用双参数PD 方法模拟了单轴拉伸荷载下含平行双裂纹薄板的裂纹扩展行为；Ha等[15]采用PD方法研究了脆性材料PMMA在单轴拉伸荷载下裂纹动态扩展和分叉的力学行为；Ran 等[16]采用改进的PD 方法进行了金属钢材的的裂纹扩展分析和由氢原子导致的氢脆断裂问题研究。

上述非局部PD 模型中，虽然可以应用黄小华等[13]和Silling[14]的PD 方法进行船用灰铸铁金属脆性材料的裂纹扩展分析，但是其微观参数键刚度系数都为常数，不能反映物质点邻域内长程力随距离增大而衰减的空间分布规律，影响计算精度。而且，上述PD 模型存在的“边界效应”问题也需要进一步研究，PD“边界效应”问题[17]，即物体内部的中心物质点具有完整的邻域，而边界处物质点的邻域却小于物体内部物质点的邻域，造成边界处物质点连接键的缺失，在计算时会产生较大误差。

因此本文基于近场动力学理论，通过校正PD的“边界效应”，考虑PD力非局部长程力的影响并推导出一种新的损伤断裂准则来研究船用灰铸铁金属脆性材料的裂纹扩展行为，预测灰铸铁构件的裂纹萌生、扩展和断裂的整个完整过程，采用近场动力学数值方法研究船用灰铸铁材料的裂纹扩展行为。本研究可为预防船用构件脆性断裂和评估船体灰铸铁结构的安全性提供参考。

1 键基PD概述

近场动力学理论[12]是使用节点力的积分来代替运动方程中的空间导数，物体被离散为一系列带质量的物质点，各物质点在有限的范围δ内彼此相互作用。假设空间域R任一物质点X与δ范围内的其他物质点X'∈R:‖‖X'-X≤δ在某一时刻t之间存在相互作用力f，则根据牛顿第二定律可得PD 方法的运动方程为

式中，ρ为材料密度，u为物质点的位移，b为外荷载密度，HX为X的近场范围，dVX'表示物质点X’处的体积微元。近场范围外物质点与X的作用力f为零。

物质点间的本构力函数f包含了材料的本构信息，也称为本构力函数，在微弹性材料[12]中，f=f(η,ξ)可用一个标量函数ω(η,ξ)表示,

式中，η为相对位移η=u'-u；ξ为相对位置ξ=X'-X；ω(η,ξ)为微观应变能密度，类似于分子动力学方法中的对势函数，用于描述近场范围内两物质点间相互作用的强弱[12]。对于线弹性材料，ω(η,ξ)可表达为

式中：c(ξ,δ)为微观模量函数，在不同维度下的表达式如表1所示，其中E为弹性模量；s为物质点对的相对伸长率，s=(|ξ+η|- |ξ|)/ |ξ|。由式（2）～（3），对点力函数f可表示为

表1 键基PD微模量函数表达式Tab.1 Micro-modulus functions of bond-based PD

式中，μ(t,ξ)是一个历史依赖的标量函数，表示物质点对的破坏情况，

式中，s0为临界伸长率。在PD理论中，通过定义标量函数φ(X,t)来反映物质点X的损伤，

式中，0≤φ(X,t)≤1。φ(X,t)=0 表示材料未损伤；φ(X,t)=1 表示该点不再与其他点相互作用，即完全损伤。

2 改进键基PD

2.1 考虑PD非局部长程力修正

Huang 等[18]和章青等[19]指出传统键基PD 模型中的微观模量c(ξ,δ)为常数，不能反映长程力随距离增大而衰减的空间分布规律，影响计算精度。本文引入长程力函数g(ξ,δ) =(1-(ξ/δ)4)，则PD 微观模量可重新表达为

方程中的本构力函数f(η,ξ)可改进为

式中，c(0,δ)为空间维数相关的材料微观模量；g(ξ,δ)为反映非局部长程力特性的空间分布函数，其强化了近距离物质点间的相互作用，弱化了远距离物质点的相互作用[19]。PD应变能密度表达式为

通过PD应变能密度与连续介质力学应变能密度相等的方法，可推导得到在不同维度下材料微观模量c(0,δ)的表达式，如表2所示。

表2 考虑非局部长程力的键基PD微模量函数表达式Tab.2 Micro-modulus functions of bond-based PD with long-range force correction

2.2“边界效应”修正

在传统键基PD 和2.1节考虑长程力特性的PD 模型中，定义PD 相互作用的参数是通过假设PD 点具有完整的邻域来计算的。然而，物体内部的中心物质点具有完整的邻域，而物体边界处物质点的邻域却小于物体内部物质点的邻域，边界点或近边界点比体内点有更少的键。因此，根据应变能密度等效的方法求得微观模量c(ξ,δ)，在物体内部和边界处均取相同的常数，弱化了边界处的微观模量值，导致计算结果不准确，即存在所谓的“边界效应”问题。本文在2.1节的PD修正模型基础上，借鉴Ganzenmuller等[20]的思想，重新推导物质点在体内和边界处的PD键刚度系数，对边界效应进行修正。

在计算边界区域的应变能密度时，仅考虑边界区域近场范围内物质点的贡献，则当发生空间或平面内各向均匀变形时，PD应变能密度为

式中，Hi为HX的离散形式，ηij为物质点Xi和Xj的相对位移，ξij为物质点Xi和Xj在参考构型中的相对位置。为了方便起见，物质点X和X'被替换为物质点Xi和Xj进行表示。因此，在2.1节考虑长程力修正的PD模型基础上，PD应变能密度可进一步表示为

令PD 应变能密度与连续介质力学应变能密度相等，可得长程力效应修正和边界效应修正后的PD微模量函数表达式，如表3所示。

表3 考虑边界效应和长程力修正后的键基PD微模量函数表达式Tab.3 Micro-modulus functions of bond-based PD with long-range force and surface effect correction

在传统PD 模型中，所有相互作用的粒子微模量常数都是相同的，因此满足fij=-fji。在当前改进的PD模型中，fij=-fji必然成立，但采用离散求和方式计算时，Hi和Hj、ci(0,δ)和cj(0,δ)不一定相同，为保证动量守恒，取

2.3 损伤断裂准则

PD 中材料的损伤是通过消除物质点之间的相互作用进行考虑的，假设当两个物质点键的拉伸s超过其临界值s0时，损伤发生。在基于键基PD 模拟断裂的方法上，键的临界拉伸值s0可以由临界能量释放率G0求得，临界能量释放率表示产生单位面积裂纹所需的能量，对于脆性材料是一个可采用试验测量的值。

由于考虑了长程力效应修正和边界效应修正，传统基于键基PD 的临界拉伸准则便不再适用，本文根据修正过的临界微观应变能密度，推导出一种新的求解物质点键临界拉伸表达式。键基PD模型临界微观应变能密度ω0(ξ,η)可表示为

在三维状态下，临界能量释放率G0的积分形式可通过ω0(ξ,η)表达为

二维状态下，临界能量释放率为

其中令板厚h=1，边界效应修正后的微模量cij替代常数微模量c进行求解，根据方程（14）和方程（15）可得不同维度下损伤断裂准则，即键的拉伸临界值sij0，如表4所示。

表4 改进键基PD的键拉伸临界值Tab.4 Critical stretch of the improved PD model

需要注意的是，物体边界或近边界点不具有完整邻域，这将导致边界或近边界点的临界拉伸小于体内点的临界拉伸值。在数值模拟中，边界可能会首先发生破坏，为了避免这种情况而更好地模拟试件真实的裂纹扩展行为，在边界处设置虚拟物质点来填充边界处不完整的邻域空间，而这些虚拟点并不参与相关的计算，使边界处物质点和内部物质的键临界伸长量保持一致。

3 算例验证

3.1 微模量系数和应变能密度验证

为了验证本文提出的微观模量计算方法的合理性，以及改进PD模型与原PD模型的应变能密度的计算精确性，考虑一边长为0.5 m、厚度为0.01 m 的正方形薄板，其弹性模量E=1 GPa，泊松比ν=1/3，取物质点间距Δ=0.01 m 进行均匀离散，采用邻域半径δ=3.015Δ进行计算。采用改进的键基PD 模型计算物质点（y=0 和y=0.25）的键刚度系数，如图1（a）所示；计算的薄板均匀变形时物质点（y=0）的应变能密度误差如图1（b）所示。

图1 微模量系数和应变能密度相对误差对比Fig.1 Comparison of micro-modulus coefficients and strain energy density relative errors

由图1（a）可知，传统键基PD模型和长程力修正的键基PD模型的微模量系数均为常数；而改进的键基PD模型计算的边界处微模量系数值不再与体内物质点的微模量系数值相同，边界处物质点的键刚度系数的数值得到强化，合理有效地削弱了“边界效应”。由图1（b）可知，传统键基PD 模型在边界处物质点的应变能密度与理论解（弹性应变能）有较大偏差；考虑长程力修正的键基PD模型减小了这一相对误差；而本文改进的键基PD 模型在此基础上进一步减小了误差，边界处物质点应变能密度的相对误差控制在0.03%左右，内部误差降至零。

3.2 非破坏误差分析

为了验证提出的改进键基PD模型的计算精度，通过数值解与解析解进行对比分析。采用长1 m、宽0.5 m的长方形薄板，弹性模量E=200 GPa，密度矢量ρ=7850 kg/m3，泊松比ν=1/3，几何参数和材料参数与黄小华等[13]采用的误差分析算例相同，薄板被离散为100×50 个物质点，物质点间距Δ=0.1 mm，邻域半径δ=3.015Δ，沿矩形短边施加的单轴应力荷载为p=200 MPa，应力边界条件的施加方式是从薄板边界区域上以体力密度的形式施加在板内物质点上。

根据弹性力学基本理论，二维平面x和y位移解析解计算公式为ux(x,y=0)=px/E，uy(x=0,y)=-νpy/E。x和y位移的数值解与解析解ux、uy的相对误差计算公式为ex=(-ux)/ux，ey=(-uy)/uy，其中，ex、ey分别表示x和y位移的相对误差。

图2（a）和2（b）分别展示了传统键基PD 模型和改进PD 模型计算的位移相对误差云图。图2（c）展示了改进键基PD模型相对误差绝对值 |ex|≥5%的物质点数量和比例。

图2 相对误差位移云Fig.2 Contours of relative errors

如图2 所示，传统键基PD 模型和改进键基PD 模型的最大相对误差分别为25.4%和6.26%，内部相对误差分别为16.3%和0.41%。采用改进PD模型得到的具有相对误差绝对值 |ex|≥5%的物质点个数和比例仅为24个和0.5%，表明了改进键基PD方法的有效性和准确性。

4 船用灰铸铁裂纹扩展模拟

4.1 含单边缺口灰铸铁裂纹扩展

为了证明改进的键基PD方法对船用灰铸铁的损伤预测和裂纹扩展模拟的能力，本文采用与李卓等[10]相同的单边缺口灰铸铁试件进行模拟计算。试件长宽厚度分别为0.2 m、0.03 m 和0.007 m，试件的左右两端0.05 m 范围内是机械夹持部位，如图3 所示，分别在单边预制长度为a=0.005 和0.01m、宽度为b=0.002 m 的裂纹。试件的材料参数分别为：弹性模量E=210 Gpa，泊松比ν=1/3，单轴抗拉强度ft=100 MPa，质量密度ρ=7150 kg/m3，断裂伸长率φ=0.3%。在PD 模拟中，物质点均匀离散尺寸为Δ=0.5 mm，物质点邻域半径δ=3.015Δ。在试件的两端均匀地施加位移控制荷载，每一步的荷载增量为Δu=1.0×10-8m，采用动态松弛法进行求解，时间步为Δt=1.0。

图3 不同初始长度单边裂纹试件的裂纹扩展示意图Fig.3 Crack propagation of single-edge specimens with different pre-crack lengths

图4 展示了采用改进键基PD 模型模拟的含5 mm 单边预制裂纹灰铸铁试件在t=0、1720、1900 和2090 时刻的裂纹扩展示意图。图5 展示了采用改进键基PD 模型模拟的含10 mm 单边预制裂纹灰铸铁试件在t=0、1350、1600和2010时刻的裂纹扩展示意图。由图4和图5可知，两种试件的破坏模式相似，裂纹首先从预切口尖端开始扩展，然后在单轴拉伸荷载下沿垂直于加载方向扩展，最后到达自由边界，两试件发生破坏。但是，含0.01 m 裂纹灰铸铁试件的起裂时刻和贯通时刻早于0.005 m 灰铸铁试件。可知，在一定条件下，预制裂纹的长度对裂纹起裂和贯通时间有较大影响，预制裂纹越大，裂纹萌生越早，试件越容易断裂。

图4 5 mm单边预制裂纹灰铸铁试件的裂纹扩展示意图Fig.4 Crack propagation of gray cast iron specimens with 5 mm unilateral notch

图5 10 mm单边预制裂纹灰铸铁试件的裂纹扩展示意图Fig.5 Crack propagation of gray cast iron specimens with 10 mm unilateral notch

图6 展示了两种不同预制单裂纹灰铸铁试件的试验结果和模拟结果对比情况，模拟结果和试验结果基本保持一致。由试验结果可知，裂纹虽然也沿垂直于加载方向扩展，但扩展呈锯齿状，而数值模拟的裂纹路径更光滑。可能的原因是，试验所用的灰铸铁材料并非均质材料，从拉断之后的断面可看出试件内部存在很多细微的孔洞[10]，而PD 模型将灰铸铁试件假设为理想均质材料，不考虑内部微观孔洞对裂纹扩展路径的影响，因此，裂纹形状对比呈现一定差异性。

图6 两种不同预制单裂纹灰铸铁试件的试验[10]和模拟破坏形态对比示意图（扩大十倍）Fig.6 Numerical results of gray cast iron and experiments of two different single-edge specimens(magnified ten times)

4.2 含双边缺口灰铸铁裂纹扩展

采用和4.1 节相同的灰铸铁材料属性，预制六种双边裂纹形式：双边裂纹长度各为5 mm，沿垂直于加载方向的对称轴进行对称分布，双裂纹尖端距离分别为5 mm、10 mm和15 mm，为方便描述，分别记为5-5 试件、5-10 试件和5-15 试件；双边裂纹长度各为10 mm，沿垂直于加载方向的对称轴进行对称分布，双裂纹尖端距离分别为5 mm、10 mm 和15 mm，分别记为10-5 试件、10-10 试件和10-15 试件。采用改进的PD 方法模拟的六种灰铸铁试件的最终破坏形态如图7所示。如图7（a）～（c）所示，5-5试件在荷载作用下，初始裂纹张开并向前扩展之后，在两条裂纹之间的区域出现了交汇串接现象，5-10试件和5-15试件的双裂纹各自独立扩展，扩展速率相同，且并无串接趋势；三种试件的破坏模式分别和10-5 试件、10-10 试件和10-15 试件的破坏模式相似，当双裂纹尖端距离为5 mm 时，裂纹串接，裂纹尖端距离是10 mm 和15 mm 时，裂纹独自沿水平方向扩展，直至破坏。这也证实了影响裂纹串接的并不是初始裂纹长度，而是裂尖纵向距离的大小，裂纹初始长度改变的是裂纹萌发和试件断裂的时间。

图7 不同预制双裂纹灰铸铁试件的数值模拟裂纹扩展示意图Fig.7 Numerical simulation of crack propagation of different double-edge gray cast iron specimens

图8 给出了六种不同预制双裂纹灰铸铁试件的裂纹扩展试验结果，与图7 数值模拟结果对比可知，试验结果和模拟结果具有很好的一致性。由于材料的非均匀性问题，导致两条裂纹的扩展速率并不相同，在同一时刻的裂纹扩展长度不同，但最终破坏模式一致，都是沿垂直于加载方向扩展，直至裂纹尖端到达试件边界处，试件破坏。

图8 不同预制双裂纹灰铸铁试件的裂纹扩展试验结果[10]Fig.8 Experimental results of crack propagation of different double-edge gray cast iron specimens

实例分析结果证明了改进的键基PD 模型能够很好地模拟灰铸铁试件从连续体到非连续体的整个破坏过程，验证了该方法用于船用灰铸铁类材料的脆性裂纹扩展和破坏分析是可行的。此外，由于采用空间积分方程求解来模拟裂纹扩展，在近场动力学模拟过程中不存在常规连续介质方法模拟时的网格二次剖分或裂尖应力奇异性等问题，随着荷载的施加，初始裂纹自然起裂、扩展和贯通。

5 结 论

本文基于近场动力学理论，提出一种改进键基PD 方法来模拟船用灰铸铁材料的裂纹扩展问题。该模型考虑了PD非局部长程力的影响，修正了PD边界效应，推导了对应的键刚度系数和断裂准则表达式，研究了含单边裂纹和双边裂纹的船用灰铸铁构件在单轴拉伸荷载下的裂纹扩展和断裂失效问题，得出以下结论：

（1）通过薄板的变形计算了改进PD 模型的键刚度值、应变能密度和物质点位移相对误差，结果显示边界处物质点的键刚度系数的数值得到强化，合理有效地削弱了PD 边界效应；边界处物质点应变能密度的相对误差控制在0.03%左右，内部误差降为零；物质点位移相对误差大大减小，有效地提高了PD模型的计算精度。

（2）改进的PD 数值模拟的船用灰铸铁构件裂纹开裂位置、扩展方向和最终破坏路径与试验结果吻合较好，证明了本文提出模型的适用性和准确性。由于试验试件内部存在微孔洞的原因，模拟和试验的部分裂纹形状虽存在一定差异，但对宏观裂纹扩展路径和最终破坏形态影响较小。

（3）预制裂纹的裂尖纵向间距是影响裂纹串接的重要因素，预制裂纹的初始长度对裂纹萌发和试件断裂的时间有较大影响。