基于非协调广义混合元的应力强度因子计算

何雨轩 卿光辉

（中国民航大学航空工程学院，天津 300300）

0 引言

裂纹损伤是引发灾难事故的主要原因之一，因此裂纹尖端附近的应力、应变和裂纹的扩展规律一直是国内外学者的研究热点问题。应力强度因子是判断已有裂纹是否扩展的物理量，高精度的应力强度因子计算结果是裂纹扩展规律分析的重要保证。

在弹性范围内断裂力学进展中最重要的成就之一是1968年Rice提出的计算应力强度因子的J积分方法［1-2］。文献［3］采用位移外推法、J 积分法和虚拟闭合技术计算裂尖应力强度因子，比较了三者的优劣。主要结论是：在有限元网格相同的条件下，J 积分和虚拟闭合技术的强度因子结果精度高。另一方面，J 积分法最大的优点在于解决了位移外推法和虚拟闭合技术高度依赖裂纹尖端区域网格细分的问题，且具有积分路径的无关性，是当前求解应力强度因子值最常用的计算方法之一。

基于位移有限元法模型研究复杂结构中裂纹或结构中复杂裂纹裂尖附近的应力、应变和裂纹的扩展规律最为普遍［3-5］。一般情况下，基于有限元模型的J 积分算法的精度取决于位移结果的精度。事实上，基于最小势能原理的位移元法用有限的结点代替了真实结构中无穷多的点，减少了结构的自由度数，因而这种模型较真实模型偏硬，给出的位移数值解是真解的下界，所以通常需要非常细的有限元网格模型才能得到比较理想的位移结果。

HERMANN［6］最早提出混合有限元法分析板壳问题。混合元模型可以同时引入位移和应力边界条件［7-8］，较位移元模型增加了已知的应力边界条件约束，使得模型的刚度更符合实际情况，进而可以提高位移变量和应力变量的精度。最近，Qing 和Tian［9］结合最小势能原理和H-R 变分原理建立了分析静力学问题的非协调的混合元。文献［10］扩展了这种混合单元的应用范畴。非协调混合元模型的结果精度明显优于非协调位移元。

在商用软件中，对平面静态裂纹的应力强度因子进行数值计算时可提供的单元类型有Abaqus 中的非协调线性元CPS4I，协调非线性元CPS8［11］等，Ansys中的平面单元PLANE42，PLANE82 等。这些单元计算的应力强度因子的方式也有多种，归纳成两大类：基于节点位移求解（外推法、J 积分）和基于节点应力求解（虚拟闭合），但由于其计算精度不高，不能很好地表现裂纹尖端的奇异性或需要额外重新设置成奇异单元［12］。本文针对应力强度因子的精度问题，构建了在商用软件中针对二维问题还未曾应用过的非协调混合元用于计算应力强度因子。首先根据广义H-R 变分原理建立含参数的非协调广义混合有限元法的数学模型，然后简要地给出J积分算法的基本理论和过程，最后通过算例验证基于非协调的广义混合有限元模型的J积分法的精确性和可靠性。

1 非协调四边形广义混合元法

假设线性弹性体的位移边界条件事先满足u-=0，则不考虑体积力的H-R 变分原理可表示为［13］

式中：σ表示应力向量；C表示材料的刚度系数矩阵；∇表示微分算子；u表示位移向量；V表示连续体的体积是作用在边界表面上的载荷。

含参数α的广义混合变分原理：

式中，参数0≤α≤1。

由式（1）和式（2）可得：

非协调的四边形单元（4 节点）的应力场和应变场可近似表示为［13-14］：

式中，N和Nr分别为单元的协调部分和非协调部分（即单元内部节点）的形函数矩阵。

将式（4）代入式（3）中，可得：

对式（5）中pe,qe,re分别进行变分，消去re后可得非协调的广义混合元的列式：

根据式（6）组装各个单元后可得到模型的控制方程。通常情况下，控制方程中的参数α取0.75（参见文献［15］）。从理论上讲，非协调广义混合元模型与混合有限元模型一样，可以同时引入位移和应力边界条件。更重要的是模型中参数α的不同取值可以起到调节模型刚度的作用，因而可以提高位移结果精确度［8］。

2 J积分理论与应力强度因子计算

2.1 J积分理论

如图1 所示，考虑围绕裂纹尖端的逆时针回路Γ，关于J积分的表达式如下：

图1 围绕裂纹尖端的逆时针积分回路Fig.1 Anti clockwise integral circuit around crack tip

式中：ui为位移矢量；ds为积分路径Γ上的微小增量；w为应变能密度因子，其定义为

式中，σij和εij分别为应力张量和应变张量。

Ti为垂直于回路的应力矢量：

式中，nj为Γ的单位法矢量。

基于有限元模型的J 积分数值算法参见文献［3，16］。

2.2 应力强度因子的计算

式中：J为J积分值；E为材料的杨氏模量；μ为泊松比。

3 数值算例

3.1 算例1:平板中心水平裂纹问题[3]

如图2所示，给定平面内宽度为2b、长度为2h的均质平板，板厚1.0 mm，有一长度为2a的中心裂缝，受到均匀对称拉伸作用，平面应力状态下，材料的弹性模量E=200×103MPa，泊松比为μ=0.3，各向同性、均匀、线弹性，该裂纹板承受均匀应力σ=30 MPa。

图2 平板中心水平裂纹Fig.2 Plate with a central crack

几何参数取a=20 mm，b=100 mm，h=200 mm。这是一个典型的Ι 型直裂纹问题，其应力强度因子KΙ的计算公式为［15］

因此，

由式（10）可得到解析解：

四分之一对称的有限元模型，如图3（a）所示。有限元网格划分数为100×200。

如图3（b）所示，对称边界左侧水平方向位移约束为零，垂直方向位移自由；底部非裂纹区域垂直方向位移约束为零，水平方向位移自由。应力边界条件如图3（c），除了已知均布应力σ，对模型右侧给出水平方向应力约束为零，垂直方向待求；底部裂纹区域竖直方向约束为零，水平方向待求。

图3 初始边界条件Fig.3 Initial boundary conditions

积分路径与文献［3］的相同，如图4 所示。本文的计算结果为243.59，与解析解的结果243.6 MPa ∙相比，误差仅为-0.27%，该误差远小于文献［3］中的-1.3%。显然，基于非协调的广义混合元模型的精度明显高于非协调位移元模型的精度。

图4 1/4对称模型中等效积分区域Fig.4 The equivalent integral region in 1/4 symmetric model

3.2 算例2:三点弯曲问题[3]

三点弯曲模型是测试材料断裂韧度的经典模型。其几何构型和加载方式如图5 所示，其中，L、W、a、S、B和P分别表示模型的长度、宽度、裂长、跨度、厚度和载荷。

图5 三点弯曲模型Fig.5 Three point bending model

取L=220 mm，W=50 mm，S=200 mm，P=30 MPa，B=1 mm，弹性模量200×103MPa，泊松比为0.3。

应力强度因子的表达为［17］

式中，形状因子函数为

有限元模型的网格划分数为220×50，非协调位移元和非协调广义混合元模型的结果与解析解的误差列于表1。

表1 非协调位移元和非协调广义混合元结果与解析解的误差（SIF）Table 1 Errors between the results and analytical solutions of nonconforming displacement element and nonconforming generalized mixed element model（SIF）

从表1 中不难看出，非协调广义混合元的精度明显高于非协调位移元。

4 结论

（1）基于广义H-R 变分原理建立了含参数的非协调广义混合元的列式。混合元模型可以同时引入位移和应力边界条件，模型更加符合实际情况。另一方面，非协调广义混合元模型中的参数α的不同取值可以调节模型的刚度，参数α合理的取值可以提高位移结果的精度。

（2）在非协调的广义混合元模型的基础上，采用J 积分法对平板中心水平裂纹和三点弯曲垂直裂纹两个典型的实例进行了一类应力强度因子分析。数值结果表明，在有限元网格模型相同的情况下，非协调广义混合元模型的应力强度因子精度明显高于非协调位移元模型。