半无限大弹性板裂纹尖端奇异应力场计算
时间：2023-01-17 17:20:11 来源：柠檬阅读网 本文已影响 人 
刘钧玉，张天禹，苏艳，宁宝宽
沈阳工业大学 建筑与土木工程学院，辽宁 沈阳 110870
在Williams[1]给出的裂纹尖端奇异应力场解析表达式中，奇异参数包括张开型（I型裂纹）应力强度因子（stress intensity factor）KI和剪切型（II型裂纹）应力强度因子KII，除此之外还包括被称为裂尖参数的T应力（T-stress）和高阶奇异参数。随着时代进步，计算机数值计算手段日趋发展，针对T应力和高阶项等奇异参数的相关研究内容逐渐增多[2-6]。裂纹尖端奇异场的参数受很多因素影响，包括结构的几何构型、材料属性、荷载类别以及裂纹倾角等，理论分析方法可以求解简单情况下的参数，但对于不同因素条件下的实际问题，一般使用数值方法来对参数进行计算。目前，计算裂纹尖端奇异应力场的数值方法主要有权函数法、有限差分法、有限元法、边界元法以及比例边界有限元法等。
有限元法需要对裂纹尖端部位进行离散，边界元法需要对裂纹面进行离散，且二者都是基于无法近似奇异点周围解析解的分片光滑函数[7]。比例边界有限元法具有有限元法和边界元法的优势。在计算裂纹尖端奇异应力场时，在裂纹面通过相似中心的情况下，可以不用离散裂纹面，这为前处理提供了便捷，使计算量减少。除此之外，比例边界有限元法相比于有限元法不需要基本解，且计算的数值结果在径向是精确的，同样也是环向收敛于有限元意义的解，裂纹尖端奇异场的这些参数由此可根据定义直接进行提取。因此比例边界有限元法在当下已在很多方面得到应用，例如对各向异性材料[7]、温度应力作用下的奇异应力场[8]的计算，以及对动应力强度因子[9]及不同材料高阶奇异性的求解[5]。文献[10]计算了劈拉试件裂纹尖端的奇异应力场，文献[11]对比例边界有限元法近年来的发展以及T应力和裂纹尖端奇异参数的计算进行了总结。
本文基于比例边界有限元法对半无限大弹性板裂纹尖端的奇异应力场进行计算，提取了裂纹尖端处的T应力以及高阶奇异项等参数。与数值求解结果进行对比验证了其精度与有效性。且对半无限大弹性板进行了断裂分析，通过改变裂纹长度和加载方式等影响因素，得出了裂纹尖端T应力以及高阶奇异参数的变化规律。在对比后可以看出T应力以及高阶奇异参数对于判断裂纹的断裂特性具有一定价值。
弹性力学平面问题的平衡方程：
式中：σ和u分别为应力和位移，ρ为质量密度，L为微分算子。比例边界坐标变换如图1所示。
图1 比例边界坐标变换
式(2)为比例边界坐标变换公式：
式中：N(η)为插值形函数，为相似中心横纵坐标，ξ、η为比例边界有限元坐标。将边界ξ=1上结点的坐标(xˆ,yˆ)用(x,y)表示。对式(1)按式(2)将坐标转换为比例边界坐标ξ、η，可得：
该方程组为二阶线性常微分方程组，沿 ξ方向的位移场{u(ξ)}可以解析求解。
系数矩阵E0、E1、E2和M0可以用数值积分方法直接计算，但是仅需在边界上进行离散计算[12-14]。
式中：M0项在动力荷载问题中存在，Nu(η)为位移插值函数。
1.1 基本方程求解
静力情况下，M0项为零。引入新变量：
式中Q(ξ)代表域内一点的等效节点力。
式(3)转化为一阶常微分方程：
式中Z为哈密顿系数矩阵，具体形式为
首先求解Z的特征值问题
式中Λ为对角矩阵。然后将Λ与Φ进行分块排列
可得：
式中λi、-λi均为Λ的一组特征值（λi的实部为负）。由此式(4)的解可设定为
式中c1、c2为积分常数。由边界处位移u(ξ=1)求解。
对于有限域来说，有限域刚度为
对于无限域来说，无限域刚度为
1.2 裂纹尖端奇异应力场
如图2所述，有裂纹的子块内保持比例边界有限元的特点。
图2 带裂纹子结构块（超单元）
式(1)中位移场的解为
式中：n为方阵Φ11的维数，Φ为特征向量矩阵，ϕi是矩阵Φ11的第i列，ci是积分常数矢量c1的第i个元素。插值函数Nu(η)可以计算单元位移场。
应力σ(ξ)为
式中：D为弹性矩阵，B1(η) 、B2(η)计算过程见文献[7]和文献[9] 。将位移式(5)代入到式(6)得到子结构内部点的应力：
式(5)和式(6)为子块内的位移和应力计算公式。为了计算裂纹尖端奇异应力场，将子块内的位移和应力的表达式表示为极坐标的形式比较方便。此时，径向坐标为
将式(8)代入到式(7)得到
将式(8) 和式(9)组合可形成Williams 以极坐标表达类似的应力场。Williams[1]裂纹尖端奇异应力场对应关系参见文献[5] 。
2.1 半无限大板裂纹边受均布拉伸荷载
图3为半无限大弹性板受均布拉伸荷载的计算模型。其中有限域第1块（相似中心O1）的宽度取W=1.5，高度H=4，裂纹长度a=1，均布力加载长度b分别取0.35、0.5、0.75和0.8。弹性模量E=1，泊松比ν=1.5。图4给出了比例边界有限元法离散网格的示意图，竖边共离散单元数N=16，第1块的相似中心坐标是(0,0)，无限域（第2块）的相似中心坐标是(-1,0)，第3块的相似中心坐标是(-1,2)。本文结果与解析解提取的高阶奇异参数a2、a3见表1。
图3 半无限大弹性板受均布拉伸荷载
图4 半无限大弹性板网格离散示意
表1 半无限大弹性板裂纹边受均布拉伸荷载（W=1.5，N=16）
由表1中数据可以看到比例边界有限元法的计算结果与解析解的最大误差为3%，这说明比例边界有限元法能够较为精确地计算应力强度因子。
2.2 半无限大弹性板裂纹边受集中剪切荷载
图5为半无限大弹性板受集中剪切荷载的计算模型。其中宽度W=2 ，高度H=4 ，裂纹长度a分别取0.7、0.8、0.9、1.0、1.1、1.2和1.3。弹性模量E=1，泊松比ν=0.25。共离散单元数N=19 。比例边界有限元法离散网格的示意图见图4，中心坐标取在裂纹尖端（0,0）。
图5 半无限大弹性板受集中剪切荷载
本文结果与文献[15]给出的解析解对比见表2。
表2 半无限大弹性板裂纹边受集中剪切荷载
由表2中数据可以看到本文计算结果与解析解的最大误差在2%以内，表明比例边界有限元法能够较为精确地计算应力强度因子。
2.3 半无限大弹性板裂纹边受均布剪切荷载
图6为半无限大弹性板受均布剪切荷载的计算模型。计算中宽度W=1.5，高度H=4，裂纹长度取a=1，均布力加载长度b分别取 0.65、 0.70、0.75、0.80、0.85、0.90和0.95 。弹性模量E=1，泊松比ν=0.25。共离散单元数取N=8。中心取在裂纹尖端坐标（0，0）。文献[15] 给出了解析解，本文结果与解析解对比见表3。
图6 半无限大弹性板受均布剪切荷载
表3 半无限大弹性板裂纹边受均布剪切荷载
由表3中数据可以看到本文计算结果与解析解的最大误差在4%以内，表明比例边界有限元法可以较为精确地计算应力强度因子。
2.4 半无限大弹性板裂纹边受集中拉伸荷载
图7为半无限大弹性板受集中拉伸荷载的计算模型。计算中宽度W=2，高度H=4，裂纹长度a=1，弹性模量E=1，泊松比ν=0.25。模型受集中力P=1。中心取在裂纹尖端坐标为（0，0）。文献[15]给出了解析解，本文结果与解析解和数值结果的对比见表4。
图7 半无限大弹性板裂纹边受集中拉伸荷载
表4 半无限大弹性板裂纹边受集中拉伸荷载计算结果
计算中离散单元数取N=11。由表4中数据可以看到本文结果a1、a2与解析解最大误差在2%以内，表明了比例边界有限元法能够较为精确地计算应力强度因子；
而高阶奇异参数a3的计算结果与文献中结果的最大误差在4%以内，表明本文方法能够较为精确地提取高阶奇异参数。
本文基于比例边界有限元法计算了半无限大弹性板裂纹尖端奇异应力场，给出了静力无限域刚度的计算方程，分析了比例边界元法位移场、应力场表达式与断裂力学裂纹尖端奇异场参数的对应关系。对半无限大弹性板在集中荷载和均布荷载作用下的裂纹尖端奇异场进行了研究，提取了应力强度因子、T 应力以及高阶项等奇异参数，将提取结果与解析解和文献中结果进行了对比，结果表明比例边界有限元法能够较精确地计算半无限大弹性板裂纹尖端应力场的奇异参数。对于拉伸荷载作用下，随着裂纹长度的增加，应力强度因子逐渐增大；
对于集中剪切荷载作用下，随着裂纹长度的增加，应力强度因子逐渐降低；
均布剪切荷载作用下应力强度因子逐渐增加。
