应用三维有限单元法计算应力强度因子
1-3应力强度因子的求解
{ε } = [ N ′( x)]{uie }
T 由应变能密度 W = ε Cε / 2 可得单元刚度
[ K ]e = ∫ [ N ′( x)]T C[ N ′( x)]dV
Ve
将所有单元组装，可得 总刚度矩阵
[ K ]{u} = {F }
节点总位移矢量
广义节点力矢量
常规单元
• 裂纹尖端应力场的奇异性要求网格划分足够细，网格尺寸一般为裂纹 尺寸的1/1000~1/100 • 求解平衡方程，得到各节点位移，取裂纹附近节点位移，根据
ξ = −1+ 2 r / l
u = ∑ N i ui
i =1 3
将二次常规单元边 中点移到1/4处
对等参单元
应变
ε=
du 1 3 l l 1 l u1 − 4 − 2 u2 + 2 − u3 = 2 − dr l 2 r r 2 r
⊙⊙⊙⊙ σ 2a
K I = σ πa 1+ κ u 2 ( x, a ) = σ a2 − x2 4µ
K II = σ πa 1+ κ u1 ( x, a) = σ a2 − x2 4µ
K III = σ πa σ u3 ( x , a ) = a2 − x2 4µ
应力强度因子的求解
一般情况
K = Yσ πa
Y为与裂纹体几何有关的参数 a为裂纹几何尺寸相关参数
σ为外载
应力强度因子的求解方法
• • • • • • • 普遍形式的复变函数方法 积分变换法 应力集中系数法 位错连续分布法 权函数法 边界配置法和边界元法 有限元法
权函数法
p(x) P*(x) p(x) P*(x) + =
应用三维有限单元法计算应力强度因子
应用三维有限单元法计算应力强度因子应用三维有限单元法计算应力强度因子摘要描述了两种基于有限单元计算面形裂纹应力强度因子的方法，建议了一种创造三维有限单元网格的途径。
计算方法的精度通过和其它解析解或数值解的比较得到了说明。
关键词应力强度因子有限元损伤容限设计断裂评定无论在损伤容限设计还是在缺陷评定阶段，工程师们需要知道正在分析的构件中裂纹的应力强度因子，因为判断含裂纹构件的断裂，或者计算剩余疲劳寿命大多依赖于这一参量。
因此，在断裂力学发展中，如何求取应力强度因子一直是一个重要的课题。
当前已有许多方法可用来计算应力强度因子，较为典型的有解析法、边界配位法、有限单元法、边界元素法、体力法、权函数法和线弹簧模型。
利用这些方法，大量的应力强度因子解已经获得，已出版的应力强度手册［1］中收编了许多典型的解。
尽管如此，工程师们仍然会感到自己所需要的应力强度因子解很难找到，这是因为要解决的工程问题往往是一些受复杂载荷的构件，包含的裂纹也往往是一些不规则裂纹。
本文简单介绍了两种基于三维有限单元法计算面形裂纹应力强度因子的方法。
有限单元法已经成为工程设计分析领域中一个强有力的计算工具，它能模拟非常复杂的构件。
基于有限元的应力强度因子计算方法，自然也将具有卓越的工程能力。
除了计算方法的介绍以外，还将简单描述一种简化网格的生成方法。
最后提供了一些所得到的典型应力强度因子解，并和大家熟知的解进行了比较，以说明本文所描述的方法的可靠性。
1 计算方法简介图1表示了笔者［2，3］建议的裂纹尖端网格形图1 裂纹前沿单元网格式。
网格由3个半环共12个20节点三维等参单元组成，每个半环有4个单元，其中第一个半环单元的中节点被移至1/4点位置，以模拟裂纹尖端应力和应变场的奇异性［4，5］。
根据线弹性断裂理论［6～8］，裂纹尖端的位移场可表示为u t=0 (3)式中，n和t分别为裂纹前沿的法线方向和切线方向；z为垂直于裂纹平面方向(见图2)；μ为材料剪切模量。
计算应力强度因子
基于ANSYS的断裂参数的计算本文介绍了断裂参数的计算理论，并使用ANSYS进展了实例计算。
通过计算说明了ANSYS可以用于计算断裂问题并且可以取得很好的计算结果。
1 引言断裂事故在重型机械中是比拟常见的，我国每年因断裂造成的损失十分巨大。
一方面，由于传统的设计是以完整构件的静强度和疲劳强度为依据，并给以较大的安全系数，但是含裂纹在役设备还是常有断裂事故发生。
另一方面，对于一些关键设备，缺乏对不完整构件剩余强度的估算，让其提前退役，从而造成了不必要的浪费。
因此，有必要对含裂纹构件的断裂参量进展评定，如应力强度因了和J积分。
确定应力强度因了的方法较多，典型的有解析法、边界配位法、有限单元法等。
对于工程上常见的受复杂载荷并包含不规如此裂纹的构件，数值模拟分析是解决这些复杂问题的最有效方法。
本文以某一锻件中取出的一维断裂试样为计算模型，介绍了利用有限元软件ANSYS计算应力强度因子。
2 断裂参量数值模拟的理论根底对于线弹性材料裂纹尖端的应力场和应变场可以表述为：其中K是应力强度因子，r和θ是极坐标参量，可参见图1，(1)式可以应用到三个断裂模型的任意一种。
图1 裂纹尖端的极坐标系应力强度因子和能量释放率的关系：G=K/E" (3)其中：G为能量释放率。
平面应变：E"=E/(1-v2)平面应力：E=E"3 求解断裂力学问题断裂分析包括应力分析和计算断裂力学的参数。
应力分析是标准的ANSYS线弹性或非线性弹性问题分析。
因为在裂纹尖端存在高的应力梯度，所以包含裂纹的有限元模型要特别注意存在裂纹的区域。
如图2所示，图中给出了二维和三维裂纹的术语和表示方法。
图2 二维和三维裂纹的结构示意图3.1 裂纹尖端区域的建模裂纹尖端的应力和变形场通常具有很高的梯度值。
场值得准确度取决于材料，几何和其他因素。
为了捕获到迅速变化的应力和变形场，在裂纹尖端区域需要网格细化。
对于线弹性问题，裂纹尖端附近的位移场与成正比，其中r是到裂纹尖端的距离。
三维裂纹扩展时动态应力强度因子的有限元研究
56
数。 如果 t=0 时的初始位移和初始速度为{u0}与
{u 0}，则由式(8)求得初始加速度
{ } ⎡⎣M ⎤⎦ {u0} = F (t) − ⎡⎣C ⎤⎦ {u0} − ⎡⎣K ⎤⎦ {u0} （11）
再根据式(9)— (11)式，求出下一时刻△t 的
{u∆t}，{u ∆t}，{u ∆t}，由此即可得到所有时间
3 结束语
(1) 利用 ANSYS 建立三维模型的方法是多 种多样的，本文主要考虑的圆柱体中心圆裂纹的 结构。还可以利用相同的单元和建模方法建立其 他典型裂纹结构，如中心穿透裂纹，椭圆裂纹等。
(2) 通 过 与 算 例 的 比 较 ， 可 以 认 为 利 用 ANSYS 的 APDL 参数化语言编程建立三维裂纹 模型是可行的，有限元法是计算动态 K1(t)因子的 有效方法，本文提供的合理计算方案，也可为其 它结构件动态 K1(t)因子的计算作参考。
关键词：动态应力强度因子 三维裂纹扩展 有限元
中图分类号：TM464
文献标识码：A
文章编号：1003-4862(2010)06-0055-04
The Finite Element Method Study on Stress Intensity Factor of three Dimensional Dynamic Crack Propagation
尖端用退化三角形奇异等参元，其余地方用标准
等参元。
表 1 材料及模型参数
材 料 E(GPa)
υ
ρ（kg/m3）
参数 200
0.3
7800
模 型 Ro(mm)
h(mm)
a
参数 100
1000
不定
本文采用突加载荷，如图 4，其中 σ0（t）=100 MPa。对不同裂纹长度的动态应力强度因子的有 限元计算结果列入表 2 中。图 5 所示为裂纹尖端 应力场，典型的应力集中现象。
计算三维问题的应力强度因子的程序
计算三维问题的应力强度因子的程序（有点简单）/COM,ANSYS MEDIA REL. 60 (090601) REF. VERIF. MANUAL: REL. 60/VERIFY,VM143*CREATE,FRACT,MAC/NOPRNSEL,ALL*GET,N,NODE,,NUM,MAX ! CURRENT MAXIMUM NODE NUMBERCMSEL,S,CRACKTIP ! SELECT THE TIP NODESESLN ! ANY ELEMENTS ATTACHED*GET,ELMAX,ELEM,,NUM,MAX ! CURRENT MAXIMUM ELEMENT NUMBER*DO,IEL,1,ELMAX ! LOOP ON MAX ELEMENTELMI=IEL*IF,ELMI,LE,0,EXIT ! NO MORE SELECTED*GET,ELTYPE,ELEM,ELMI,ATTR,TYPE ! GET ELEMENT TYPE*IF,ELTYPE,NE,ARG1,CYCLE ! CHECK FOR SELECTED ELEMENTN3 = NELEM(ELMI,3) ! GET NODE 3 (K)*IF,NSEL(N3),LE,0,CYCLE ! IT MUST BE SELECTEDN7 = NELEM(ELMI,7) ! GET NODE 7 (L)*IF,NSEL(N7),LE,0,CYCLE ! IT MUST ALSO BE SELECTEDN1 = NELEM(ELMI,1) ! GET NODE 1 (I)N2 = NELEM(ELMI,2) ! GET NODE 2 (J)N5 = NELEM(ELMI,5) ! GET NODE 5 (M)N6 = NELEM(ELMI,6) ! GET NODE 6 (N)X3 = 0.75*NX(N3) ! WEIGHTED POSITION OF N3Y3 = 0.75*NY(N3)Z3 = 0.75*NZ(N3)X = 0.25*NX(N2) + X3 ! QUARTER POINT LOCATION ( NODE (R) )Y = 0.25*NY(N2) + Y3Z = 0.25*NZ(N2) + Z3N = N + 1 ! NEXT NODEN10 = NN,N10,X,Y,Z ! MIDSIDE NODE LOCATIONX = 0.25*NX(N1) + X3Y = 0.25*NY(N1) + Y3Z = 0.25*NZ(N1) + Z3N = N + 1N12= NN,N12,X,Y,ZX7 = 0.75*NX(N7)Y7 = 0.75*NY(N7)Z7 = 0.75*NZ(N7)X = 0.25*NX(N6) + X7 Y = 0.25*NY(N6) + Y7 Z = 0.25*NZ(N6) + Z7 N = N + 1N14 = NN,N14,X,Y,ZX = 0.25*NX(N5) + X7 Y = 0.25*NY(N5) + Y7 Z = 0.25*NZ(N5) + Z7 N = N + 1N16 = NN,N16,X,Y,ZN4=N3N8=N7NSEL,ALLTYPE,3EN,ELMI,N1,N2,N3,N4,N5,N6,N7,N8 ! REDEFINE THE ELEMENTEMORE,0,N10,0,N12,0,N14,0,N16EMORE,*ENDDOCMSEL,U,CRACKTIP ! UNSELECT THE TIP NODES NUMMRG,NODE ! MERGE MIDSIDE NODES NSEL,ALL ! SELECT ALL ELEMENTSESEL,ALL ! SELECT ALL ELEMENTS/GOPR*END/PREP7*afun,degInnerRadius=0.1 !InnerRadius为裂纹半径OuterRadius=1 !OuterRadius为圆柱半径Scaler=0.025 !Scaler为裂纹前沿单元范围，已经是最佳BaseHeight=0.51 !BaseHeight为基层高度,it may be the bestLayerHeight=0.18 !LayerHeight为扩展层高度,差别越大可能越好LayerAmount=16 !LayerAmount为层数,与精度关系不大RotationAngle=6 !RotationAngle为单元旋转的角度Rotationtimes=90/RotationAngle !Rotationtimes为旋转的次数Height=LayerAmount*LayerHeight+BaseHeight !Height为总高度SMRT,OFF/TITLE, VM143, FRACTURE MECHANICS STRESS INTENSITY - CRACK IN A FINITE WIDTH PLATEC*** BROWN AND SRAWLEY, ASTM SPECIAL TECHNICAL PUBLICATION NO. 410./COM, ****** CRACK IN 3-DIMENSIONS USING SOLID45 AND SOLID95ANTYPE,STATIC ! STATIC ANALYSISET,1,SOLID45ET,2,SOLID45 ! ELEMENTS AROUND THE CRACK TIP ET,3,SOLID95 ! CRACK TIP ELEMENTS CREATED USING MACRO FRACTMP,EX,1,2e4MP,NUXY,1,.3 ! CYLINDRICAL COORDINATE SYSTEM local,33,1,InnerRadiuscsys,33N,1NGEN,9,20,1N,11,ScalerN,171,Scaler,180FILL,11,171,7,31,20local,33,0,InnerRadiuscsys,33FILL,1,11,9,2,1,9,20,3N,15,OuterRadius-InnerRadiusN,75,OuterRadius-InnerRadius,BaseHeightFILL,15,75,2,35,20N,155,-InnerRadius,BaseHeightNGEN,2,200,155 !155点在轴线上，故而在此与355为同一点，目的在于产生轴线上的单元FILL,75,155,3,95,20N,172,-InnerRadiusNGEN,2,200,172 !且因为无论什么坐标对此无影响，所以可以直接产生，与155点类似FILL,155,172,5,177,-1,,,.15ngen,2,200,173,177,1 !将轴线上的节点依次产生其同一位置的节点，便于产生轴线单元csys,33FILL,11,15,3,,,7,20,3csys,5ngen,2,200,1,177,,,RotationAngle !在此只是产生了一次的节点，并不是21次（180度）csys,0E,2,22,1,1,202,222,201,201 !产生裂纹边沿的一个单元EGEN,8,20,-1 !旋转生成八个单元，均在裂纹边沿E,2,3,23,22,202,203,223,222 !裂纹次边沿EGEN,8,20,-1 !旋转生成八个单元EGEN,9,1,-8 !由此边沿向外扩展生成另外8层的单元TYPE,3EMODIF,1 ! MODIFY ELEMENTS 1 TO 8 FROM TYPE,1 TO TYPE,2*REPEAT,8,1NUMMRG,NODE ! MERGE COINCIDENT NODEScsys,33NSEL,S,LOC,X,0NSEL,R,LOC,Y,0CM,CRACKTIP,NODE/NERR,0 ! TEMPORARILY NO WARNINGS OR ERRORS PRINTOUT! (IN ORDER TO AVOID WARNING MESSAGES DUETO! MIDSIDE NODES LOCATION)FRACT,2 ! CONVERSION MACRO, TYPE 2 IS SOLID45 ! ELEMENTS AROUND THE CRACK TIPcsys,0type,1e,171,371,172,172,151,351,173,173 !产生左边的单元，轴线附近（最底层）e,151,351,173,173,131,331,174,174 !产生左边的单元，轴线附近（次底层）e,131,331,174,174,132,332,175,175 !产生左边的单元，轴线附近（第三层）egen,3,1,-1 !产生左边的单元，轴线附近（映射至第五层）e,134,334,177,177,135,335,155,155 !产生左边的单元，轴线附近（第六层最高层）e,11,12,32,31,211,212,232,231 !产生右边的单元，第一层第一列egen,4,1,-1 !产生右边的单元，第一层第二列至第五列e,31,32,52,51,231,232,252,251 !基本同上，为第二层，第一列egen,4,1,-1e,51,52,72,71,251,252,272,271 !基本同上，为第三层，第一列egen,4,1,-1e,71,72,92,91,271,272,292,291 !基本同上，为第四层，第一列egen,4,1,-1e,91,92,112,111,291,292,312,311 !基本同上，为第五层，第一列egen,4,1,-1e,111,112,132,131,311,312,332,331egen,4,1,-1csys,5egen,Rotationtimes,200,1,110,1,,,,,,,RotationAnglecsys,0nsel,s,loc,y,BaseHeightngen,2,5000,all,,,,LayerHeighte,135,335,155,155,5135,5335,5155,5155csys,5egen,Rotationtimes,200,-1e,135,115,315,335,5135,5115,5315,5335egen,Rotationtimes,200,-1e,115,95,295,315,5115,5095,5295,5315egen,Rotationtimes,200,-1e,95,75,275,295,5095,5075,5275,5295egen,Rotationtimes,200,-1csys,0nsel,s,loc,y,BaseHeight+LayerHeightesln,s,0,allegen,LayerAmount,5000,all,,,,,,,,,LayerHeightnsel,allesel,all/NERR,DEFA ! TURN ON THE WARNINGS OR ERRORS PRINTOUT/OUTPUTOUTPR,,ALLcsys,0NSEL,S,LOC,z,0DSYM,SYMM,z ! SYMMETRIC B.C.'S AT X = 0d,all,uznsel,allnsel,s,loc,y,0csys,5nsel,r,loc,x,InnerRadius,OuterRadiusDSYM,SYMM,Y ! SYMMETRIC B.C.'S AT Y = 0 EXCEPT CRACK NODESd,all,uynsel,allcsys,0nsel,s,loc,x,0dsym,symm,xd,all,uxNSEL,S,LOC,Y,HeightSF,ALL,PRES,-1000NSEL,ALLESEL,ALLFINISH!now is the solution/OUTPUT,SCRATCH/SOLUSOLVEFINISH/OUTPUT/POST1C*** IN POST1 DETERMINE KI (STRESS INTENSITY FACTOR) USING KCALC !**csys,33 !define the local systermPATH,KI1,3,,48 ! DEFINE PATH WITH NAME = "KI1" PPATH,1,1 ! DEFINE PATH POINTS BY NODE PPATH,2,170PPATH,3,171KCALC,,,1 ! COMPUTE KI FOR A HALF-MODEL WITH SYMM. B.C.*GET,KI1,KCALC,,K,1 ! GET KI AS PARAMETER KI1老大，你的裂纹处的单元怎么也是8节点的45号单元啊，这样算出来的应力强度因子对不对？呵呵，那里的单元是8节点的SOLID45单元,但是那里的单元已经经过奇异单元的处理,所以结果正确，我与标准的结果进行了比较,答案是一致的pengfmBlueblueday同志：您好！我看了您所给的例题很受启发，但还有个问题想向您请教。
T型焊接件的表面裂纹应力强度因子的计算
中图分类号：U661
文献标识码：A
文章编号：1006-7973（2021）02-0057-02
引言 T 型焊接接头是工程结构中常采用的典型结构，其疲劳 断裂性能严重影响着焊接结构的使用安全[1]。应力强度因子 是表征裂纹尖端应力场的重要参量[2]，研究含裂纹结构强度 时必须计算裂纹处应力强度因子大小，以判断结构工作时的 安全状态。随着有限元理论的发展和有限元软件功能的完善， 现广泛借助有限元软件求解应力强度因子[3]。 在实际工程结构中，由于疲劳、缺陷、应力集中等会出 现各种形式的裂纹，表面裂纹是其中常见的裂纹形式。本文 以 T 型焊接件为研究对象进行焊缝表面裂纹的热-结构耦合 应力强度因子的计算，利用稳态热分析模块和静力学分析模 块对焊接件进行稳态热应力分析，可为同类型应力强度因子 的计算提供借鉴，对工程的实际应用有重要的应用价值。 二、T 型焊接件有限元模型建立 1．T 型焊接件结构分析 本文中 T 型焊接件为结构钢直接焊接而成，T 型焊接件 所用材料为工程常用材料 Q235A 钢，其材料参数见表 1 所 示，结构模型如图 1 所示。
表 1 Q235A 钢材料参数
参数
值
参数
值
弹性模量/MPa
2.1
屈服强度/MPa
≥235
密度/（kg/m³）
7,900
泊松比
0.277
2．表面裂纹的插入 将用 Solidworks 软件建立的三维立体模型导入 ANSYS 软件中，工程上 3D 表面裂纹的形状一般采用半椭圆来描述。 在一侧的焊缝表面建立局部坐标系，此局部坐标系用来给裂 纹定位，插入表面半椭圆裂纹（工程中最常见的裂纹类型）， 裂纹的长半轴（裂纹长度）为 0.008m，裂纹的短半轴（裂 纹深度）为 0.003m，影响区域为 0.008m，裂纹前缘的径 向划分为 20 份，周向划分为 16 份（需为 8 的倍数），最大 绕线半径为 0.0035m。 3．网格划分 利用 Ansys 分析求解获得足够精度的应力强度因子的 关键是单元的选取和网格的划分。进行网格划分，由于焊接 件结构简单，直接用 workbench 进行四面体网格主导方法 进行划分。 利用有限元法求解应力强度因子，裂纹前沿周围网格的 划分是关键[3]。为了兼顾模拟结果的精度，裂纹所在面需要 更精密的网格划分，基体的网格划分太粗糙，所以用局部控 制 Sizing 设置该焊接面网格尺寸 0.001m，重新细化该面网 格，提高计算精度，T 型焊接件模型网格划分结果如图 2 所 示。
三维裂纹前缘应力强度因子数值计算方法
0 引 言
船舶与海洋工程结构在建造和营运过程中不可避 免地存在一些裂纹，这些裂纹在交变载荷作用下将会 发生扩展，最终可能会导致结构发生断裂失效，严重 影响结构的安全使用。采用数值方法预报含裂纹结构 裂纹扩展形貌及扩展寿命具有重要的工程应用价值， 可为结构裂纹检测周期的确定、疲劳寿命评估及抗疲 劳设计提供支撑。结构裂纹前缘应力强度因子（SIF） 是含裂纹结构疲劳裂纹扩展形貌及寿命数值预报方法
· 10 ·
舰船科学技术
第 41 卷
型应力强度因子的计算方法。采用平板表面裂纹研究 三维裂纹整个裂纹前缘各个位置应力强度因子随裂纹 网格划分参数变化情况，通过整个裂纹前缘应力强度 因子数值计算结果与 Newman-Raju 解析公式[6]计算结果 的对比分析，验证了数值计算方法的准确性，最终提 出裂纹网格划分参数取值范围。为了验证文中提出的 数值计算方法计算其他类型裂纹前缘应力强度因子的 准确性，采用数值计算方法计算分析了疲劳裂纹扩展 三点弯曲试样裂纹扩展过程中裂纹前缘应力强度因 子，并与解析解进行对比分析，验证了数值计算方法 计算穿透性裂纹前缘应力强度因子的准确性。
文中基于 20 节点奇异单元 1/4 节点位移法，提出 三维裂纹网格划分方法及裂纹前缘不同位置、不同类
收稿日期: 2018 – 06 – 30 基金项目: 国防预研项目资助项目 (3020503010101) 作者简介: 黄如旭 (1987 – )，男，工程师，研究方向为船舶与海洋工程结构疲劳与断裂理论、实验及数值模拟。
中的关键参数，一些学者提出了采用奇异单元求解裂 纹尖端应力强度因子的有限元模型，并以裂纹最深点 及裂纹端点处的应力强度因子作为计算参量模拟简单 焊接接头裂纹面内扩展形貌[1 – 5]。然而，由于工程结构 承受载荷复杂性及结构中裂纹分布随机性致使结构最 大主应力与裂纹平面并不垂直，裂纹扩展形貌将不再 保持平面。结构三维裂纹前缘各个位置处应力强度因 子准确计算将是模拟真实裂纹扩展形貌的前提。
基于ANSYS的应力强度因子计算
基于 ANS YS 的应力强度因子计算
王 伟
(湖北超高压输变电公司 ,武汉 430050)
摘 要 : 以平板表面裂纹为例 ,介绍和分析了运用有限元软件 ANSYS 计算应力强度因子的方法 。通过对求得的
应力强度因子值与解析解的比较 ,表明用有限元方法计算应力强度因子具有相当高的精度 ,并且操作简便 。
表 1 应力强度因子
裂纹表面
裂纹最深处
分析结果 公式结果
22 22. 4
36. 1 36. 7
由表 1 可知 ,ansys 计算结果与应力强度因子手
册公式计算结果相差均在 5 %以内 ,分析结果略大
于公式结果 ,能满足精度及安全储备要求 。
4 结 语
通过以上分析可以看出 ,应力强度因子在 AN2 S YS 中的计算是可行的 。与传统的利用断裂力学中 的公式直接求解相比 ,利用有限元软件 ANS YS 来 求解应力强度因子是一种简单而准确的方法 。
[ 4 ] 瞿伟廉 ,鲁丽君 ,李 明. 带三维穿透裂纹结构的有限 元实体建模方法[J ] . 武汉理工大学学报 ,2008 ,30 (1) : 87290 .
收稿日期 :2010203224. 作者简介 :王 伟 (19842) ,硕士. E2mail :wwwj132 @163. com
(上接第 72 页)
关键词 : ANSYS ; 表面裂纹 ; 应力强度因子
Calculate the Stress Intensity Factor by the Finite Element Sof t ware ANSYS
W A N G Wei ( Hubei EHV Transmission & Substation Company ,Wuhan 430050 ,China)
基于FRANC3D计算三维疲劳裂纹最大应力强度因子
Y U Shi ( Co lle g e
o
g
u a n
e r s
Z HU R
o
p
s
g
d A
te
s
f Me
e r
c
ha
ta
n
ic
a
l E
n
g in
a
e e r
in g
r e e
，
Na
—
n
j in
g Un iv
ity
m o
f A
e r o n a u
tic
t r o n a u t ic s
r e n
，
，
。
当前 已 有 许 多
。
种计 算应 力 强 度 因 子 的方 法 如 边 界 元 法 有 限元 法 等
，
、
裂纹前沿应力
，
随 着 边 界元 有 限 元 理 论 的发 展 以 边 界 元 法 有 限 元
、
，
、
强 度 因子 的计算方法 可
以直接利用边界元理 论
法求 解 应 力 强 度 因 子 得 到 了 广 泛 的应 用
r v
e
le
m e n
t m e th o
d ；i n
te n s
ity fa
c
to
r
( S I F )；F R A N C 3 D
工 程 结 构 经 常 存 在表 面 裂纹 和 缺 陷
，
例 如 齿 轮结 构 作
，
软件 相 继 涌 现 出来
一
为传 递 功 率 和 运 动 的 重 要 零 件 经 常会 发 生 齿 根 裂 纹 现 象
#有限元法在应力强度因子计算中应用96884
有限元法在应力强度因子计算中的应用本文构建了含裂纹平板的二位1/4模型及三维1/2模型，分别用于计算张开型、滑移型和撕开型裂纹尖端应力强度因子。
有限元分析的结果在误差范围内可以较好地与解读解吻合。
计算结果表明，裂纹周向单元的划分会严重影响有限元计算的结果，三维模型沿厚度方向提高划分密度可以有效提高计算精度。
1前言项目分析中，材料中的裂纹会对结构可靠性带来很大的影响。
历史上有很多航空航天事故、建筑事故都是因为裂纹引起的断裂导致结构失效的。
为了检验结构是否能够一般用于判断裂纹是否延伸的重要判据就是应力强度因子K<Stress Intensity Factor，SIF）。
对于任何材料，其应力强度因子极限K C只与材料本身的属性有关，而和裂纹尺寸、裂纹周围应力强度无关。
在具体的项目分析中，评估含裂纹结构稳定性，只需要计算含裂纹结构在要求的工况下的裂纹尖端应力强度因子K值，若K>Kc，则裂纹会发生扩展，导致结构失效。
具体工况下，应力强度因子K的计算可以通过多种方法完成。
弹性力学给出了三种基本断裂模式<分别为张开型、滑移型和撕开型，见图1）的应力强度因子解读解。
但是对于一般几何结构而言，求解读解的复杂程度会随结构的复杂程度成倍增加，很多情况甚至无法求出解读解。
有限单元法及边界元法可以用于具有复杂几何结构的含裂纹模型分析。
有限单元法中，经常使用的方法包括1/4节点位移法和J积分法。
本文使用有限元法建立裂纹的二维、三维模型，分别求解I型、II型和III型裂纹的应力强度因子，并对有限元法求解应力强度因子的精确性进行分析。
FFF<a） <b） <c）图1张开型<a）、滑移型<b）和撕开型<c）裂纹的受力状态2有限元法求解应力强度因子理论分析有限元法求解应力强度因子的方法包括1/4节点法、J积分法等方法，其中，J积分法计算应力强度因子时的计算过程较为繁杂，不便于多次重复求解分析，因此本文选择1/4节点法计算裂纹尖端应力强度因子K值。
三维表面裂纹应力强度因子参数敏感性有限元分析
三维表面裂纹应力强度因子参数敏感性有限元分析赵慧;赵玮;宋磊;刘伟【期刊名称】《郑州航空工业管理学院学报》【年(卷),期】2024(42)3【摘要】三维表面裂纹是许多工程结构中普遍存在的一种缺陷,无论是在损伤容限设计或缺陷评估阶段,人们都要求对裂纹的应力强度因子进行精确的估计,这一参量对裂纹的预测和剩余疲劳寿命的估计起到至关重要的作用。
借助有限元软件Abaqus和断裂力学软件Franc3d建立三维半椭圆表面裂纹实体模型,依据断裂力学相关知识,采用M-积分法,对含有表面裂纹的模型进行数值模拟,研究了模型在不同裂纹形状比(a/b=1.5、1.36、1.25和1.15)以及裂纹倾角(θ=0°、15°、30°和45°)下的裂纹扩展情况和应力强度因子参数敏感性。
主要研究内容与结论如下:裂纹形状比a/b越大,应力强度因子越小,应力强度因子分布形貌由“凹”向“凸”转变;裂纹倾角达到0°时,裂纹深处应力强度因子最大,裂纹更容易扩展,随着裂纹倾角(0°、15°、30°和45°)的增大,应力强度因子逐渐减小。
【总页数】8页(P18-25)【作者】赵慧;赵玮;宋磊;刘伟【作者单位】西安航空学院飞行器学院;太原国际机场有限责任公司通信导航室;山西太钢不锈钢股份有限公司炼钢二厂【正文语种】中文【中图分类】V215.6【相关文献】1.T型管节点表面裂纹应力强度因子的全三维有限元分析2.埋头紧固件椭圆形表面裂纹应力强度因子的三维有限元分析3.含表面半椭圆裂纹板裂纹尖端应力强度因子的三维光弹分析4.CD350套管吊卡应力场三维有限元分析及危险截面上裂纹应力强度因子求解5.SL450水龙头中心管螺纹齿根表面裂纹应力强度因子三维有限元分析因版权原因，仅展示原文概要，查看原文内容请购买。
应力强度因子的求解方法的综述
应力强度因子的求解方法的综述摘要：应力强度因子是结构断裂分析中的重要物理量，计算应力强度因子的方法主要有数学分析法、有限元法、边界配置法以及光弹性法。
本文分别介绍了上述几种方法求解的原理和过程，并概述了近几年来求解应力强度因子的新方法，广义参数有限元法，利用G*积分理论求解，单元初始应力法，区间分析方法，扩展有限元法，蒙特卡罗方法，样条虚边界元法，无网格—直接位移法，半解析有限元法等。
关键词：断裂力学；应力强度因子；断裂损伤；Solution Methods for Stress Intensity Factor of Fracture MechanicsShuanglin LU(HUANGSHI Power Survey&Design Ltd.)Abstract: The solution methods for stress intensity factor of fracture mechanics was reviewed, which include mathematical analysis method, finite element method, boundary collocation method and photo elastic method. The principles and processes of those methods were introduced, and the characteristics of each method were also simply analyzed in this paper.Key words: fracture mechanics; stress intensity factors0 引言断裂力学的基础理论最初起源于1920年Griffith的研究工作[1]。
Griffith在研究玻璃、陶瓷等脆性材料的断裂现象时，认为裂纹的存在及传播是造成断裂的原因。
应力强度因子
应力强度因子断裂与损伤力学应力强度因子数值计算方法综述2013年6月第一章应力强度因子求解方法概述含有裂纹的工程结构的断裂力学分析一直是一个重要问题，在断裂力学理论中应力强度因子是线弹性断裂力学中最重要的参量。
它是由构件的尺寸、形状和所受的载荷形式而确定。
由于裂尖应力场强度取决于应力强度因子，因此在计算各种构件或试件的应力强度因子是线弹性断裂力学的一项重要任务。
是边界条件。
这些边界条件则是通过在相应于真实裂纹体的假想边界上施加一系列的集中力和集中力矩来满足的，先把假想的边界离散化为一组线段，在每一段的中心，在离开假想边界处加上一对集中力和力矩，这些力和力矩的值可通过近似地满足边界条件得以确定。
与其他数值方法相比，边界力法有其明显的优点。
由于这一方法已精确地满足了裂纹面上的边界条件，所以它不需要像边界元法那样把裂纹面视为边界的一部分。
另外，它也克服了边界配位法中所需要的对每一类裂纹问题都要建立新的应力函数的缺点。
这种解法只要较小的自由度就能达到相当高的精度。
因此它在求解几何形体复杂的裂纹向题中有着明显的优点，但在处理复杂载荷的能力方面，则远非如权函数法那样灵活。
2.5权函数法权函数法是一种求解在任意受载条件下裂纹应力强度因子的高效方法。
这种解法的高效性在于它把影响应力强度因子的两个因素，即载荷与几何，作了变量分离。
权函数仅反映了裂纹体的几何特性，它可以根据一种受载情况下的已知解确定。
一经导出，它就能被用来不受限制地求解任意加载条件下的k 值，求解中只需作一个积分运算：dx x x a m K a )().,(0σ⎰=式中m(a,x)为权函数，)(x σ为无裂纹体中假想裂纹处的应力分布。
除了灵活通用，简单经济等特点外，这一方法所得的结果有高的可靠性。
2.6 有限元法有限元法在断裂力学中有着非常广泛的应用，它不受解析方法常遇到的因裂纹体几何或载荷的复杂性的限制。
这种方法的基本思路是用一系列离散化的，区段连续的场变量来对任何连续的场交量作逼近。
基于有限元法对裂纹尖端应力强度因子的计算
基于有限元法对裂纹尖端应力强度因子的计算杨巍;张宁;许良【摘要】基于ANSYS有限元软件通过相互作用积分法建立了求解三维穿透裂纹应力强度因子的有限元模型,将有限元法和解析法求得的应力强度因子进行比较验证了模型的准确性.研究了载荷、裂纹长度、试样宽度、厚度对裂纹尖端应力强度因子的影响,在对比结果的基础上分析了裂纹尖端应力强度因子的三维效应.结果表明:在不同条件下有限元模型都可以很好的模拟出应力强度因子的值,二维状态时应力强度因子的分布规律与三维状态时的分布规律有较大差异,出于安全的考虑不应忽略应力强度因子的三维效应,对三维应力强度因子的有限元求解有一定的指导意义.【期刊名称】《沈阳航空航天大学学报》【年(卷),期】2014(031)003【总页数】5页(P19-23)【关键词】相互作用积分方法;三维穿透裂纹;应力强度因子;ANSYS【作者】杨巍;张宁;许良【作者单位】中航工业沈阳飞机工业(集团)有限公司制造工程部,沈阳110136;沈阳航空航天大学机电工程学院,沈阳110136;沈阳航空航天大学机电工程学院,沈阳110136;沈阳航空航天大学航空制造工艺数字化国防重点实验室,沈阳110136【正文语种】中文【中图分类】O346.1断裂是工程构件最危险的一种失效方式，尤其是脆性断裂，它是突然发生破坏，断裂前没有明显的征兆，这就常常引起灾难性的事故。
大量断裂事故分析中发现，断裂皆与结构中存在缺陷或裂纹有关。
由缺陷或裂纹所引起的机械、结构的断裂失效,是工程中最重要、最常见的失效模式[1]。
传统强度理论是建立在假设材料无缺陷的基础上的,但工程实际中很多按传统强度理论设计的结构由于裂纹的产生和扩展,出现大量断裂失效。
为保证含裂纹构件的安全性和可靠性,必须预测裂纹的扩展速率和构件的断裂强度,在断裂力学的工程应用中,应力强度因子K是判断含裂纹结构的断裂和计算裂纹扩展速率的重要参数。
应力强度因子的计算方法有解析法、有限元法、边界元法、权函数法等[2]。
ansys计算应力强度因子
本文使用ANSYS13.0中的互动积分法（Interaction Integrals ）计算了三维贯穿裂纹的应力强度因子，计算结果表明该方法计算可靠，为计算更复杂的三维裂纹提供了一种途径。
据一些工业化国家统计，因材料和结构的破坏所造成的损失占国民经济生产总值的8% -12%多。
破坏事故所造成的人员伤亡的损失更不可估量。
我国作为一个发展中国家，在这方面的情况比西方发达国家更严重。
因此无论是为了减少破坏事故的损失还是研发满足现代工业所需要的新材料，都要求对材料的破断过程有科学的、全面的、定量化的认识。
三维裂纹作为工程中常见的裂纹形式，早在六十年代初就有不少研究者开始研究，到现在已有大量的文献资料论及这一问题，出现了一些有特点的分析方法。
工程上常见的表面裂纹的断裂分析，由于其实质是三维问题，也几乎同时开始被人们所关注。
三维裂纹问题的危害极大，断裂造成了大量的灾难性事故发生，这使得断裂力学在机械工程、海洋工程、核工程，特别是今天的航空航天工程中受到更广泛的重视和深入研究。
因此对含三维裂纹结构断裂特性尤其对三维裂纹体的应力强度因子的研究有重要的现实意义。
本文使用ANSYS成功的计算了三维贯穿裂纹的应力强度因子，为计算三维裂纹提供了一种便捷方式。
1.模型的建立图1 三维贯穿裂纹模型本文三维裂纹模型长度为L，高度为H，宽度为W，裂纹半长为a，裂纹位于模型的中心部位。
几何参数见表1。
模型的为线弹性材料，其弹性模量为2.1E11Pa，泊松比为0.3。
模型的边界条件为：底端固定，顶端承受拉应力σ为2E6Pa。
表1 模型的几何参数本文采用二维奇异单元PLANE183建立二维的裂纹模型，然后通过拉伸并使用三维奇异单元SOLID186来建立三维贯穿裂纹模型。
图2-图5给出了二维裂纹模型和三维裂纹模型。
在13.0中对应力强度因子的计算增加了一种计算方法即互动积分法(Interaction Integral s )，这种方法与计算J积分的主域积分法类似。
基于ANSYS的三维应力强度因子的计算
基于ANSYS的三维应力强度因子的计算摘要:本文在总结断裂力学各种行为研究方法的基础上，采用有限单元法思想，利用ANSYS软件建立裂纹体有限元模型。
通过计算得出I型裂纹尖端的应力强度因子，其计算结果与理论值吻合良好。
表明模型的选取和网格的划分是合理的，具有可靠的精度，其结果完全可以指导工程设计。
关键词：断裂力学；有限单元法；ANSYS；应力强度因子1 前言随着现代生产技术的高速发展，新材料、新工艺在航空、航天、压力容器、核反应堆、机械、土木等领域得到了广泛应用，结构在高速、高温、高压等环境中使用时，按照传统强度理论设计的结构在应用中却出现大量的断裂事故。
目前关于断裂力学的研究，理论上只能求解较简单的模型或做出较强的假设条件；通过实验探求其规律性的成本较高、周期较长。
因此对于断裂力学问题，尤其是三维裂纹问题目前大多借助数值方法进行研究。
裂纹问题严格来说都是三维问题，并且工程中最后发生事故的裂纹问题的物体(如机械构件，土建结构)总是有限弹性的。
因此对有限弹性裂纹问题进行三维分析在实际工程上有重要的现实意义[1]。
迄今在平面断裂力学中已形成了一整套相当成熟的计算方法，但在三维断裂问题中，尤其在表面裂纹方面，还有很多问题有待进一步探讨，本文正是在这方面进行了探索和研究。
通常板表面裂纹应力强度因子可以统一表示为：(1)只是不同的解给修正因子F赋予不同的表达式。
由Newman的研究成果可以看出，他主要考虑了拉伸或弯曲载荷下半椭圆表面裂纹的各种裂纹形状（a/c）和板宽、板厚对裂纹前缘应力强度因子的影响，而没有考虑板长的影响，也就是没有考虑平行边界对应力强度因子的影响。
2 裂纹类型与有限元法2.1 裂纹类型在断裂力学中，按裂纹受力情况和裂纹面的相对位移方向将裂纹分为三种基本类型[2]：即张开型（I型），滑移型（II型），撕裂型（III型），如图1。
两种或三种基本类型的组合称之为复合型裂纹问题。
图1中所示三种基本类型的裂纹模型的受力特点如下：I型裂纹受垂直于裂纹面的拉应力作用；II型裂纹受平行于裂纹面而垂直于裂纹前缘的剪应力作用，又称为面内剪切型；III型裂纹受既平行于裂纹面又平行于裂纹前缘的剪应力作用，对应于反平面剪切，又称为面外剪切或纵向剪切。
构造应力场三维数值模拟的有限单元法
构造应力场三维数值模拟的有限单元法的报告，800字
有限单元法（Finite Element Method, FEM）是一种使用数值方法来求解复杂的场计算问题的重要技术。
本文以构造应力场的三维数值模拟为例，详细讨论有限单元法技术的原理和实现步骤。
首先，基于地形，我们将复杂地形分割成简单的解析片，并计算物体表面上每个节点的位置，然后根据定义的有限元形状函数，把节点和元素之间的连接形成单元，其中物体表面上每个单元都可以表示为一个复合形式。
此外，基于对空气力学的几何观念，我们运用多种有限元形式（如旋转四面体、梯形和直角四面体），以实现物体的三维空间剖分，并计算三维空间内每个单元中的体积，所有单元的体积总和即为物体的总体积。
然后，根据应力场的数值模拟需求，采用有限单元法计算每个单元上每个节点的应力分布，即原有三维场集中在单元中的结构、状态，从而计算出表面上每个节点的应力分布。
通常情况下，我们会定义一套有限元函数，作为换算关系，以解决每个单元上各种节点的应力计算问题。
接着，采用非线性的差分迭代求解的方法，假定场的应力场状态发生变化，求解方程组，这套方案非常有效，可以大大减少计算量，并保证计算结果的准确性。
最后，得到精确的应力场三维数值结果。
有限单元法技术在构造应力场三维数值模拟中有着不可替代的
作用，本文正是对有限单元法的原理及实现步骤进行讨论，从而为应用者提供了一种有效的计算解决方案。
7085铝合金三维中心斜裂纹应力强度因子的数值计算
7085铝合金三维中心斜裂纹应力强度因子的数值计算宋彦琦;刘小珍;王卓;王石磊【摘要】To explore stress intensity factor of fatigue crack propagation on 7085-T7485 aluminum alloy,the fi-nite element model of aluminum alloy with three-dimensional center oblique crack was established by ANSYS.The stress intensity factor of the crack is calculated by the displacement extrapolation method and the interaction integral method.The results were compared with the theoretical results, and analysis the influence of length and angle of center crack on stress intensity factor.The results show that the two methods are both close to theoretical results in a certain range,but with the length and angle's increase of crack,the interaction integral is more precision.There-fore,the stress intensity factor of crack propagation is calculated by the interaction integral method.The results are significance to study the three-dimensional fatigue crack propagation mechanism of aluminum alloy.%为研究7085-T7485铝合金疲劳裂纹扩展的应力强度因子,用ANSYS软件建立了含三维中心斜裂纹铝合金的有限元模型.利用位移外推法和相互作用积分法计算了裂纹的应力强度因子;并将两种方法计算得到的结果与解析解对比,分析了中心斜裂纹的长度和倾斜角度对应力强度因子的影响.结果表明,在一定的范围内,两种方法与解析解都比较接近;但随着裂纹长度和角度的增加,相互作用积分凸显出更高的精度.因此,采用相互作用积分法对裂纹扩展的应力强度因子进行了计算,研究结果对铝合金三维复合型疲劳裂纹扩展机制研究有一定的指导意义.【期刊名称】《科学技术与工程》【年(卷),期】2017(017)028【总页数】5页(P251-255)【关键词】7085-T7452铝合金;应力强度因子;ANSYS模拟;位移外推法;相互作用积分法【作者】宋彦琦;刘小珍;王卓;王石磊【作者单位】中国矿业大学(北京)力学与建筑工程学院,北京100083;中国矿业大学(北京)力学与建筑工程学院,北京100083;中国矿业大学(北京)力学与建筑工程学院,北京100083;中国矿业大学(北京)力学与建筑工程学院,北京100083【正文语种】中文【中图分类】TU512.4航空航天需要采用高性能大型整体高强铝合金锻件实现轻量化和高可靠性[1,2]，新型7085-T7452高强铝合金锻件具有高强度、高淬透性、高损伤容限及耐腐蚀性;并已成功应用于波音787飞机和空客A380飞机的翼梁、翼肋等重要承力构件。
三维裂纹应力强度因子数值计算
摘要：利用无网格 Galerkin 方法对三维裂纹问题进行分析。基于最小二乘法的无网格方法仅需一系列节点信息就
可构成物体的离散模型。在无网格 Galerkin 方法中，位移边界条件不包含在内，因此，在边界已知位移方向上设
置一系列弹簧来满足无网格方法的本征边界条件。弹簧的一端与已知位移的物体表面固连，而另一端固定，把计
收稿日期：2006–01–06；修回日期：2006–02–08 基金项目：国家自然科学基金资助项目(40272120)；国家杰出青年科学基金项目(50229901) 作者简介：张敦福(1963–)，男，硕士，1993 年于清华大学工程力学系一般力学专业获硕士学位，现任副教授，主要从事计算力学、振动理论与应 用方面的教学与研究工作。E-mail：zhangdf@
• 3836 •
岩石力学omain of influence. This method significantly increases the efficiency of computing approximate functions by limiting the size of the least-square problem. Three-dimensional discontinuous displacement method is used to evaluate stress intensity factors along the 3D crack front. Applications of the method to the determination of stress intensity factors along single edge planar cracks and single through edge planar crack in 3D finite bodies are presented. The obtained stress intensity factors for both problems are found to be in good agreement with SIF values reported in previous studies. It can guarantee the success for trace propagation of three-dimensional crack. Key words：rock and soil mechanics；stress intensity factor；three-dimensional crack；visibility criterion； three-dimensional discontinuous displacement method；enriched element-free Galerkin method
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
应用三维有限单元法计算应力强度因子
摘要描述了两种基于有限单元计算面形裂纹应力强度因子的方法，建议了一种创造三维有限单元网格的途径。
计算方法的精度通过和其它解析解或数值解的比较得到了说明。
关键词应力强度因子有限元损伤容限设计断裂评定
无论在损伤容限设计还是在缺陷评定阶段，工程师们需要知道正在分析的构件中裂纹的应力强度因子，因为判断含裂纹构件的断裂，或者计算剩余疲劳寿命大多依赖于这一参量。
因此，在断裂力学发展中，如何求取应力强度因子一直是一个重要的课题。
当前已有许多方法可用来计算应力强度因子，较为典型的有解析法、边界配位法、有限单元法、边界元素法、体力法、权函数法和线弹簧模型。
利用这些方法，大量的应力强度因子解已经获得，已出版的应力强度手册［1］中收编了许多典型的解。
尽管如此，工程师们仍然会感到自己所需要的应力强度因子解很难找到，这是因为要解决的工程问题往往是一些受复杂载荷的构件，包含的裂纹也往往是一些不规则裂纹。
本文简单介绍了两种基于三维有限单元法计算面形裂纹应力强度因子的方法。
有限单元法已经成为工程设计分析领域中一个强有力的计算工具，它能模拟非常复杂的构件。
基于有限元的应力强度因子计算方法，自然也将具有卓越的工程能力。
除了计算方法的介绍以外，还将简单描述一种简化网格的生成方法。
最后提供了一些所得到的典型应力强度因子解，并和大家熟知的解进行了比较，以说明本文所描述的方法的可靠性。
1 计算方法简介
图1表示了笔者［2，3］建议的裂纹尖端网格形
图1 裂纹前沿单元网格
式。
网格由3个半环共12个20节点三维等参单元组成，每个半环有4个单元，其中第一个半环单元的中节点被移至1/4点位置，以模拟裂纹尖端应力和应变场的奇异性［4，5］。
根据线弹性断裂理论［6～8］，裂纹尖端的位移场可表示为
u t=0 (3)
式中，n和t分别为裂纹前沿的法线方向和切线方向；z为垂直于裂纹平面方向(见图2)；μ为材料剪切模量。
图2 曲线裂纹前沿坐标系
1.1 1/4点位移法
根据式(2)，若裂纹表面上(θ=180°)某一点的垂直于裂纹平面的位移已知，那么应力强度因子
(4)
用从有限单元法求出的1/4点位移u z(1／4)(见图1)求应力强度因子的方法即为1/4点位移法，即
(5)
显然这一方法很简单，但要求裂纹尖端附近的应力应变场能较好地被模拟。
1.2 三维J积分法
J积分法实际上是一种能量方法，近来被广泛用来计算应力强度因子，因为两者可按下式转换
用有限元计算J积分的方法通常有2种：回路积分法与虚拟裂纹扩展法。
后者最早由Parks［9］和Hellen［10］独立提出，主要是通过移动有限元模型的节点位置来模拟裂纹扩展。
最近deLorenzi
［11，12］已经根据连续介质力学成功地推导了面形裂纹的能量释放计算公式
式中符号请参见文献［11］。
商用有限元分析软件ABAQUS已经收编了上述基于“虚拟裂纹扩展原理”的三维J积分方法，可以直接获得面形裂纹前沿各角节点和中节点的J积分值。
2 网格生成
为了方便三维网格的生成，笔者建议了图3所示的分块方法。
整个网格由裂纹块和非裂纹块
图3 结构分块示意图
组成，裂纹块包含了裂纹，其网格需要进行精心划分，因为它对应力强度因子有相当大的影响。
相对于裂纹块，非裂纹块显得不是很重要，一般网格可以划分得粗一些，但必须使载荷和约束有一个比较合理的分配。
裂纹块和非裂纹块网格的组合可采用“多点位移约束法”，使其位移在其接触面上基本保持连续，从而使组合后的网格保持变形协调。
对于复杂含裂纹构件，只要合理地设计裂纹块和非裂纹块即可大大简化网格的生成。
这种网格形式也特别适合于疲劳裂纹扩展的模拟计算，疲劳裂纹扩展后的新网格重新生成可以只对裂纹进行，而保持非裂纹块网格不变［2，3］。
图4表示了一个笔者建议的三维裂纹块网格的生成过程。
首先生成一个8节点等参单元二维网格，在这个网格里定义开裂表面和裂纹前沿，然后将其扩展成一个三维网格，单元的扩展有2种，见图4。
一个非裂纹前沿单元扩展成一个20节点单元，而一个裂纹前沿单元扩展成6个20节点单元。
围绕裂纹前沿，网格呈放射形。
最后组合成的裂纹尖端网格具有图1所示的结构。
显然，这种网格扩展方法进一步简化了三维网格的建立。
图4 三维裂纹块网格的形成
3 计算程序
笔者最近发展了如图5所示的计算应力强度因子的计算机程序。
对于一些常见的工程裂纹，如平板中的半椭圆表面裂纹、压力容器中的半椭圆表面裂纹、平板孔边的角裂纹、圆棒中的表面裂纹等，只要输入一些关键的几何参数和载荷参数，程序将自动输出应力强度因子结果。
该程序目前主要和商用软件ABAQUS相连接，在UNIX操作系统下运行。
图5应力强度因子计算流程图
4 结果比较
4.1 无限大体中的椭圆形埋藏裂纹
图6表示了沿椭圆形埋藏裂纹前沿的应力强度因子结果。
椭圆裂纹短长轴之比为a/c=0.5。
结果表明，根据位移法和J积分法获得的结果非常相近。
两者的误差约为1%。
图中的虚线为Irwin ［6］的解析解，本文结果略小于Irwin的解，最大误差小于3%。
图6 埋藏椭圆裂纹前沿应力强度因子变化
4.2 有限厚度平板中的半椭圆表面裂纹
近20年来，这一裂纹得到了深入的研究，大量的应力强度因子结果已经发表［13］，其中，Newman_Raju［14］根据他们自己的有限单元分析结果提出的应力强度因子方程，被认为具有比较可靠的精度，并且易编程。
这一方程已经在工程中得到了广泛的应用。
图7比较了笔者和Newman_Raju的结果，笔者的结果由1/4点位移法获得。
裂纹形状为半圆形a/c=1。
结果表明，无论在裂纹深度点还是在表面点，笔者的结果和Newman—Raju的解误差约在5%以内，全面的比较见文献［2］。
图7 平板中半椭圆表面裂纹前沿应力强度因子比较
(a)拉伸载荷(b)弯曲载荷
4.3 压力容器中的半椭圆形内表面裂纹
图8所示为压力容器筒体的局部网格及应力强度因子结果。
筒体的壁厚和内径之比为
t/R i=0.1，载荷为均匀内压力。
Raju—Newman［15］的结果也示于图8中。
结果表明，笔者的结果
和Raju—Newman的结果基本一致，只是当裂纹深度达到4/5的壁厚时，Raju—Newman的结果相对低一些。
从图8a中可以看出，筒体网格由裂纹块和非裂纹块组成。
进一步的比较见文献［16］。
图8 受内压容器内表面裂纹
(a)网格(b)应力强度因子比较
4.4 孔边的半椭圆表面裂纹
图9比较了笔者的结果［17］和Shivarkumar—Newman［18］及Zhao等人［19，20］的结果。
Shivarkumar—Newman的结果基于有限元分析，而Zhao等人用的是权函数方法。
可见，结果同样符合得较好。
图9 孔边表面裂纹应力强度因子比较
5 结束语
本文描述了两种基于有限单元法计算面形应力强度因子的方法，建议了一种裂纹尖端的网格形式，提供了一种能方便创造包含裂纹构件三维网格的思路，最后比较了几个工程中常见裂纹的应力强度因子结果。
结果表明，笔者的结果和一些已被广泛接受的应力强度因子解符合得相当好，这说明本文所描述的方法具有可靠的计算精度，可用来计算复杂工程裂纹的应力强度因子。