有限元法在应力强度因子计算中的应用

1/7FF有限元法在应力强度因子计算中的应用马源[e-mail]（原创）摘要：本文构建了含裂纹平板的二位1/4模型及三维1/2模型，分别用于计算张开型、滑移型和撕开型裂纹尖端应力强度因子。有限元分析的结果在误差范围内可以较好地与解析解吻合。计算结果表明，裂纹周向单元的划分会严重影响有限元计算的结果，三维模型沿厚度方向提高划分密度可以有效提高计算精度。关键词：1有限元应力强度因子精度分析1前言工程分析中，材料中的裂纹会对结构可靠性带来很大的影响。历史上有很多航空航天事故、建筑事故都是由于裂纹引起的断裂导致结构失效的。为了检验结构是否能够一般用于判断裂纹是否延伸的重要判据就是应力强度因子K（StressIntensityFactor，SIF）。对于任何材料，其应力强度因子极限KC只与材料本身的属性有关，而和裂纹尺寸、裂纹周围应力强度无关。在具体的工程分析中，评估含裂纹结构稳定性，只需要计算含裂纹结构在要求的工况下的裂纹尖端应力强度因子K值，若KKc，则裂纹会发生扩展，导致结构失效。具体工况下，应力强度因子K的计算可以通过多种方法完成。弹性力学给出了三种基本断裂模式（分别为张开型、滑移型和撕开型，见图1）的应力强度因子解析解。但是对于一般几何结构而言，求解析解的复杂程度会随结构的复杂程度成倍增加，很多情况甚至无法求出解析解。有限单元法及边界元法可以用于具有复杂几何结构的含裂纹模型分析。有限单元法中，经常使用的方法包括1/4节点位移法和J积分法。本文使用有限元法建立裂纹的二维、三维模型，分别求解I型、II型和III型裂纹的应力强度因子，并对有限元法求解应力强度因子的精确性进行分析。（a）（b）（c）图1张开型（a）、滑移型（b）和撕开型（c）裂纹的受力状态1马源机械71班2007010423Email：[e-mail]2/72有限元法求解应力强度因子理论分析有限元法求解应力强度因子的方法包括1/4节点法、J积分法等方法，其中，J积分法计算应力强度因子时的计算过程较为繁杂，不便于多次重复求解分析，因此本文选择1/4节点法计算裂纹尖端应力强度因子K值。使用ANSYS有限元分析软件计算K值时，只需要将裂纹尖端进行1/4节点划分（图2），进行求解，再沿着裂纹设定相应的路径，即可自动求解出该模型、边界条件下的裂纹尖端应力强度因子值。弹性力学给出裂纹尖端应力强度因子的解析解与1/4节点位移的关系如（1）式。(1/4)22(1)uKr(1)其中，(1/4)u为1/4节点的位移。图2二维、三维单元1/4节点3有限元分析模型的构建为了对张开型、滑移型和撕开型裂纹应力强度因子进行有限元求解，本文分别构建了二维1/4模型（PLANE82）和三维1/2模型（SOLID45及SOLID95）。其中，PLANE82为平面8节点单元，SOLID95为三维20节点单元。（图3）。图3平面1/4模型和三维1/2模型3/7图中，网格最密集的中心即为裂纹尖端（如箭头所指），裂纹尖端向左侧即为裂纹。为了使用1/4节点法计算应力强度因子，需要对裂纹尖端周围的单元进行1/4节点设置（见图2）。二维单元设置1/4节点可以使用KSCON命令实现。该命令可以控制裂纹尖端周围单元的划分。包括周向单元数以及第一周单元的尺寸。三维模型无法使用KSCON命令对单元进行自动划分，因此需要通过直接设定节点坐标再构建单元的方法进行建模。考虑到裂纹两侧表面在分析中设为距离为0，则两层节点会拥有相同的坐标。层与层之间的节点也会出现坐标重合的现象。为了保证计算精度，需要将层与层见的重合节点进行融合，同时还要保证裂纹两侧的节点不被融合。这个过程是通过节点的选择完成的。4裂纹尖端应力强度因子的计算4.1张开型裂纹应力强度因子的计算张开型裂纹常被记为I型裂纹，弹性力学给出了无限大平面上张开型裂纹尖端应力强度因子的计算公式IKa(2)实际问题中，大部分模型不可能具有无限大的外形，因此有限宽度板件上的裂纹尖端应力强度因子往往更能够说明问题。有限宽度板件的裂纹尖端应力强度因子修正公式是sec2IaKab(3)此情况下的受力状态如图4图4有限宽度平板张开型裂纹计算模型简图为以上模型进行如下赋值计算：一块宽度为2b,高度为2h,其中间沿宽度方向有长度为2a的裂纹，平板沿高度方向两侧承受拉力。材料：E=63010psi,υ=0.3几何：a=1in,5bin,h=5inσσ2a2b4/7载荷：0.5641895psi分别使用上文所述的平面模型和三维模型进行分析计算时，可以得到KI计算结果如表1表1平面、三维模型求解张开型裂纹尖端应力强度因子结果KI理论解数值解误差二维模型1.02491.05873.3%三维模型1.02490.99223.2%4.2滑移型裂纹应力强度因子的计算滑移型裂纹常被记为II型裂纹，弹性力学给出了无限大平面上张开型裂纹尖端应力强度因子的计算公式IIKa(2)现有文献并没有针对有限宽度平板滑移型裂纹尖端应力强度因子的计算。从4.1的张开型裂纹修正值与原始解相差不大，可以使用无限大平板的解析解来近似替代。本问题中的受力状态如图5图5有限宽度平板滑移型裂纹计算模型简图为以上模型进行如下赋值计算：一块宽度为2b,高度为2h,其中间沿宽度方向有长度为2a的裂纹，平板沿裂纹方向两侧承受剪切应力。材料：E=63010psi,υ=0.3几何：a=1in,5bin,h=5in载荷：1.08psi分别使用上文所述的平面模型和三维模型进行分析计算时，可以得到KI计算结果如表25/7表2平面模型求解滑移型裂纹尖端应力强度因子结果KII理论解数值解误差二维模型1.9141.9361.2%4.2滑移型裂纹应力强度因子的计算撕开型裂纹常被记为III型裂纹，弹性力学给出了无限大平面上张开型裂纹尖端应力强度因子的计算公式IIIKa(2)现有文献并没有针对有限宽度平板滑移型裂纹尖端应力强度因子的计算。从4.1的张开型裂纹修正值与原始解相差不大，可以使用无限大平板的解析解来近似替代。本问题中的受力状态如图6图5有限宽度平板滑移型裂纹计算模型简图为以上模型进行如下赋值计算：一块宽度为2b，高度为2h，厚度为t，其中间沿宽度方向有长度为2a的裂纹，平板沿裂垂直于平板方向承受剪切应力。材料：E=63010psi,υ=0.3几何：a=1in,5bin,h=5in载荷：1.00psi分别使用上文所述的平面模型和三维模型进行分析计算时，可以得到KI计算结果如表3表3平面模型求解滑移型裂纹尖端应力强度因子结果KIII理论解数值解误差三维模型1.77251.76613.6%6/75有限元划分方式对分析精度的影响由于裂纹尖端对于线弹性断裂力学而言属于奇点，而且要求尖端周围的单元进行特殊的1/4节点划分，因此容易推知裂纹尖端的划分情况会对有限元分析精度带来很大的影响。通过改变裂纹尖端的划分方式，可以推出数值分析精度随划分的改变的改变结果。根据图6可以看出，改变周向单元划分可以从以下两个方面讨论：1、周向单元划分数2、周向划分半径图6裂纹尖端局部的网格划分因此，在同样的载荷条件下，使用平面滑移型裂纹模型，分别改变裂纹尖端周向划分数和周向单元半径，得到结果如图7、图8。图7周向单元数对结果的影响图8周向划分半径对结果的影响周向单元数周向划分半径7/7从图7、图8中可以看出，随着周向划分单元的减少，误差呈增大趋势。同时，若周向划分半径过小，则计算精度大幅下降。从两条曲线的趋势也可以发现，对于平面滑开型裂纹，存在着可以消除的系统误差。周向单元数与周向划分半径对三种裂纹的影响方式是相同的。对于撕开型裂纹，平面应变模型不再使用，因而平板的厚度必然会对结果分析产生影响。因此，在选用厚度合适的模型下，改变厚度方向的划分密度，可以得到图9所示结果。图9厚度方向划分密度对分析精度的影响从图9可以看出，随着划分密度的提高，有限元分析精度越来越靠近解析解。6结论分别构建了裂纹的二维1/4模型及三维1/2模型，用于分析张开型、滑移型和撕裂型裂纹的裂纹尖端应力强度因子分析。从有限元分析的结果可以看出，这些模型的构建是有效的。使用有限元模拟的方法，得到的大部分结果精度与理论值的误差在4%以下。随着网格划分的不同，数值计算的结果会出现大幅度的波动，模型对网格附近的划分参数敏感度很大。对于裂纹周向划分，一般情况下需要将划分半径控制在合理范围内，不可过小，周向划分单元数可以尽量多。对撕开型裂纹应力强度因子计算，厚度方向划分密度对裂纹应力强度因子有影响，其中，提高划分密度可以有效提高计算精度。参考文献1陆毅中,工程断裂力学,西安交通大学出版社,19872ANSYSHelp周向划分半径