ANSYS积分法和节点位移法求解应力强度因子附命令流

基于ANSYS数值求解I-II混合模型裂纹应力强度因子的计算方法目前对于求解断裂力学中应力强度因子的数值计算方法以交互积分法为主，常用的通用有限元软件ANSYS以及ABAQUS在此方面已能进行很好的求解，得出比较精确的结果。今天，作者介绍基于ANSYS的两种应力强度因子的数值计算方法。1.交互积分法计算应力强度因子1.1交互积分方程的表达首先，介绍一下目前基于ansys的应力强度因子计算的数值方法——交互积分法，其定义如下：I=−∫qi,j(σklεklauxσij−σkjauxuk,i−σkjuk,iaux)dV/∫δSqndsV其中：𝜎𝑖𝑗,𝜀𝑖𝑗,𝑢𝑖分别为应力，应变及位移；𝜎𝑖𝑗𝑎𝑢𝑥,𝜀𝑖𝑗𝑎𝑢𝑥,𝑢𝑖𝑎𝑢𝑥分别为辅助场的应力，应变及位移；𝑞𝑖为裂纹扩展矢量。交互积分和应力强度因子的关系如下：I=2E∗(K1K1aux+K2K2aux)+1μK3K3aux其中：Ki(i=1,2,3)为I、II、III型应力强度因子；Kiaux(i=1,2,3)为辅助场的I、II、III型应力强度因子；𝐸∗=𝐸/(1−𝜈2)，对于平面应变问题，μ为剪切模量，E为泊松比，ν为泊松比。1.2基于ansys的交互积分法应力强度因子的计算本计算算例模型为混合型裂纹的计算，模型材料参数如下图1-1所示：图1-1模型尺寸及材料参数Ansys有限元模型裂尖采用1/4节点的奇异单元，其网格划分及，模型以有限元建立如下图1-2、1-3所示：图1-2裂尖单元的网格划分图1-3有限元模型在建立有限元物理模型的过程中，值得注意的是，裂纹面的处理应该为上下两条重合的线，但各自隶属于不同的面。裂纹面的建立如下图1-4所示：图1-4裂纹面的处理Ansys计算裂尖应力强度因子采用交互积分法时，采用的命令为CINT，首先得在以裂尖为原点，X轴平行于裂纹面，Y轴垂直于裂纹面的局部坐标系。本算例该部分的命令流为：/SOLULOCAL,11,0,A*SIN(ALPH),A*COS(ALPH),,BETA,,,1,1ALLSEL,ALLCINT,NEW,1!DEFINECRACKIDCINT,TYPE,SIFS!CALCULATESTRESS-INTENSITYFACTORCINT,CTNC,CRACK_TIP!DEFINERIGHTCRACKTIPNODECOMPONENTCINT,SYMM,OFF!SYMMETRYOFFCINT,NORMAL,11,2!DEFINECRACKPLANENORMALCINT,NCON,6!NUMBEROFCOUNTOURSANTYPE,STATICSOLVE/POST1PRCINT,1,CRACK_TIP,K1PRCINT,1,CRACK_TIP,K2本算例围道积分次数为6次，在使用交互积分进行求解应力强度因子时，通常会设置好多次围线积分（NUMBEROFCOUNTOURS），除了第一道围线积分所得的值舍去外，计算出的值求平均值即可。图1-5、1-6所示，各次围线积分的结果及第一主应力云图。最终计算得到的应力强度因子为：KI=228.56Mpa*(mm)1/2，KII=221.95Mpa*(mm)1/2带斜裂纹的平板模型应力强度因子的计算公式为：KI=𝜎√𝜋𝑎𝑠𝑖𝑛2𝛽𝐾II=𝜎√𝜋𝑎𝑐𝑜𝑠2𝛽即解析解为：KI=217.08Mpa*(mm)1/2，KII=217.08Mpa*(mm)1/2计算误差分别为：5.29%，2.24%。得到的结果比较精确。图1-5交互积分计算结果图1-5第一主应力云图1.3基于ansys的节点位移法计算应力强度因子节点位移法为ansys计算应力强度因子的另一种数值方法，该方法首先在裂尖节点处建立局部坐标系，其X轴平行于裂纹面，Y轴垂直于裂纹面。最后计算裂纹面同一点的相对水平位移Δu及竖向位移Δv，其图解如图1-5所示。混合裂纹模型的应力强度因子计算表达式为：KI=√2πG1+κ∆v√rKII=√2πG1+κ∆u√r其中：G为剪切模量；κ为材料常数，对于平面应力问题，取3−𝜈1+𝜈；∆u为裂纹面在某点处的水平相对位移；∆v为裂纹面在某点处的垂直相对位移。图1-6位移法图解根据断裂力学对于三种裂纹的定义，当∆v0时，KI为正，裂纹上下面相对位移为顺时针为正，即顺时针时，∆u0，KII为正；反之为负。理论上，当取上下裂纹面同一位置的点，当该点趋向于裂尖时，结果更精确，本算例取奇异单元上1/4处的节点的位移进行计算，计算模型同上。首先，先对有限元模型进行求解，然后进入到后处理层，求出在局部坐标系系下，所处裂纹上下面的奇异单元上1/4处节点的水平及竖直位移ux，uy，然后求出裂纹面的相对位移∆u、∆v，最后代入上式即可。计算结果如图1-7所示：KI=223.84Mpa*(mm)1/2，KII=217.63Mpa*(mm)1/2。计算误差分别为：3.1%、0.25%。图1-71/4节点位移法计算得到应力强度因子1/4节点法的节点至裂尖的距离为：a/48，其中a为倾斜裂纹半长。若取距离裂尖a/6处的节点的位移值，计算得到的结果如下：KI=218.98Mpa*(mm)1/2，KII=212.48Mpa*(mm)1/2距离裂尖不同距离的节点位移值计算出来的应力强度因子如下表所示：距裂尖的距离（mm）a/48a/12a/8a/65a/24KI/Mpa*(mm)1/2223.84227.78221.36218.98216.47KII、Mpa*(mm)1/2217.63219.86214.99212.48210.04以上结果说明在1/4节点处计算应力强度因子比较精确，距裂尖越远，计算越不精确，当然和网格划分的精密度也有很大关系。Ansys中的KCALC命令是运用位移外推法来计算应力强度因子的，该方法在计算过程中需要定义一个由五个点确定的路径，这五个点为裂尖节点，上裂纹面两个节点以及下裂纹面相对应位置的两个节点，该过程由命令PATH来完成。位移外推法的理论计算公式为：KI=√2πG1+κ|∆v|√rKII=√2πG1+κ|∆u|√r与1/4节点位移法不同之处在于，位移外推法是运用取极限的方式求得的应力强度因子的计算表达式，所以无法判断应力强度因子的正负。当然，对于I型应力强度因子，不存在负值，其原因在于物质的相互重合叠加在物理上是不允许的。本算例取距离裂尖a/48、a/12两处四个节点和裂尖节点组成该路径用位移外推法来计算混合型裂纹的应力强度因子。得到的结果如下：KI=251.75Mpa*(mm)1/2，KII=242.42Mpa*(mm)1/2以下计算不同倾角下，分别采用交互积分法和1/4节点位移法事计算的结果跟理论解进行比较。图1-8不同倾角下应力强度因子的计算结果对比2．讨论交互积分法的原理理解起来相对困难，但计算应力强度因子对网格划分精度的要求不高，然而这种计算是基于积分与路径无关的前提下，对于和积分路径有关的计算，其积分表达式就不成立，计算的结果自然就不正确。对于位移外推法，其计算结果的精确度不高，这跟网格划分的精密度有很大关系，在此也说明了位移外推法对网格精密度的要求比1/4节点位移法所要求的更高。因为1/4节点位移法在满足和交互积分法同一网格划分精度的情况下，所得到的应力强度因子更加精确，并且不会有其它的限制条件，方法简单容易理解，精度较高，普适性较广。FINISH/CLEAR/TITIE,INTERACTIVEINTEGRATIONMETHODBYIDUTER-ANSYS/PREP7/RGB,INDEX,100,100,100,0/RGB,INDEX,80,80,80,13/RGB,INDEX,60,60,60,14/RGB,INDEX,0,0,0,15/REPLOT!-------------------!UNIFIEDUNIT(N,MM)PI=ACOS(-1)*SET,H,80!THEHEIGHTOFMODEL*SET,W,50!THEWEIGHTOFMODEL*SET,A,0.12*W!HALFLENGTHOFTHEANGLEDCRACK*SET,BETA,90-45!THEINCLINEDANGLEOFCRACK*SET,ALPH,(90-BETA)*PI/180!RADIANSYSTEM*SET,SIGMA,100!SIGMAR1=1!FIRSTROWOFELEMENTRADIUSR2=2!THIRDROWOFELEMENTRADIUSR3=3!SIXTHROWOFELEMENTRADIUS!----------------------------DEFINEMATERIALPARAMETERET,1,183KEYOPT,1,3,0MP,EX,1,10E9MP,PRXY,1,0.3MP,MU,1,COF!---------------------------CREATFINITEELEMENTMODELWPSTYLE,,,,,,,,1WPROTA,BETA,,CSWPLA,11,0RECTNG,-(A+R3+R2),A+R3+R2,0,R2+R3RECTNG,-(A+R3+R2),A+R3+R2,-(R2+R3),0CYL4,-A,0,R3,0,,180,,CYL4,-A,0,R2,0,,180,,ASBA,1,3,,DELE,KEEPASBA,3,4,,DELE,KEEPCYL4,-6,0,R3,0,,-180,,CYL4,-6,0,R2,0,,-180,,ASBA,2,3,,DELE,KEEPASBA,3,6,,DELE,KEEPKSEL,S,LOC,X,-(A+R3+R2),-ANUMMRG,KPALLSEL,ALL,ALLCYL4,6,0,R3,0,,180,,CYL4,6,0,R2,0,,180,,ASBA,5,3,,DELE,KEEPASBA,3,8,,DELE,KEEPCYL4,6,0,R3,0,,-180,,CYL4,6,0,R2,0,,-180,,ASBA,7,3,,DELE,KEEPASBA,3,10,,DELE,KEEPKSEL,S,LOC,X,A,A+R3+R2NUMMRG,KPALLSEL,ALL,ALL!-------------------WPROTA,-BETA,,CSYS,0RECTNG,-18,18,-18,18!MIDDLERECTANGLEALLSELASBA,3,ALL,,DELETE,KEEPRECTNG,-25,25,-40,40!BIGRECTANGLEALLSELASBA,3,ALL,,DELETE,KEEP!-------------------KSCON,14,0.5,1,8,1KSCON,20,0.5,1,8,1!-------------------TYPE,1LESIZE,12,,,8LESIZE,19,,,8AMESH,4AMESH,6LESIZE,1,,,2LESIZE,24,,,2LESIZE,7,,,2AMAP,1,13,10,12,9AMAP,2,13,10,18,15!-------------------LESIZE,21,,,8LESIZE,30,,,8AMESH,8AMESH,10LESIZE,16,,,2LESIZE,22,,,2LESIZE,28,,,2AMAP,5,11,19,8,17AMAP,7,16,23,8,17!-------------------AESIZE,9,1AMESH,9AESIZE,11,1AMESH,11AESIZE,12,2AMESH,12AESIZE,13,2AMESH,13!-------------------NSEL,S,LOC,X,-25NSEL,R,LOC,Y,-40D,ALL,ALL,0ALLSEL,ALLNSEL,S,LOC,Y,-40D,ALL,UY,0ALLSEL,ALLNSEL,S,LOC,Y,40SF,ALL,PRES,-SIGMAALLSEL,ALL!--------