分享一个断裂分析实例（J积分和应力强度因子）（最终版）
有些地方没有加注释。加起来比较麻烦。后面不再修改了。建议参考help关于fracture部分以及vm143 部分。
有更新以上附件没有删除。
一、 问题描述
如图所示为一等厚度空心圆盘（参看图1.1），厚度4mm，内径r=10+Δ，Δ=22mm，外径R=500mm，材料属性数据：双线性（参看图1.3）E=2.1×1011Pa，ET=6.0×109Pa，屈服极限σs=500MPa，μ=0.3，密度ρ=8500kg/m3，采用Mises屈服准则，裂纹初始长度为a0=2.0mm。裂纹如图1.2，裂纹位于0度90度180度和270度的位置，图中粗短线表示裂纹。
要求用有限元解
（1） 载荷为均布拉力q=150MPa时的KI。（参看图1.1）
（2） 载荷为均布拉力q=1200MPa时的J积分值。（参看图1.1）
（3） 当载荷为转速n=600r/min的KI ，并计算沿裂纹尖端不同路径的积分值，与KI比较。
二、 求解
求解中统一采用国际单位制，长度m，压力、应力与弹性模量Pa，密度Kg/m3，转速rad/s。
对圆盘的1/4进行ANSYS建模，网格划分如图2.1。单元类型为6节点三角形单元plane2。裂纹附近单元边长为0.0002m。载荷施加如图2.2，扇形两条半径（裂纹处除外）上施加对称位移边界条件，弧上加均布拉力。裂纹处无位移约束。
计算的应力强度因子和J积分结果如表2.1。前面给出的J积分由于坐标系错误，做法不对，现在改正1200Mpa下的J积分结果。但是，结果差别不大。因为，即使在局部坐标系下，J积分中用到的 XG, YG, ZG的坐标还是在全球坐标系下的。坐标系的改变只对J积分的第一项有一点影响（对Y坐标积分这部分，Y坐标变了）。
新加两个附件：改正后求解1200MPa下J积分的命令流，及两张路径定义的图片命令流里面涉及到路径定义的命令是GUI方式进行的。所以要分开看。。
表2.1 应力强度因子和J积分结果
载荷 项目 不同路径下的结果
150MPa KI 3.1359e7 2.986e7 3.3068e7 2.7695e7
（此处J积分没改） J-Integral 3193.416293182.07917 3245.22042 3167.23352
1200MPa KI 5.4207e9 5.3398e9 5.382e9 5.285e9
J-Integral 4701408.68 4777189.41 4674756.71 4709657.93
600r/min KI 55011 54901 54672 54463
（此处J积分没改） J-Integral 1.5649802 1.5719383 1.53405748 1.57113766
三、一些需要讨论的问题（后面补的帖解决了大部分问题，这些留在这里供参考）：
1、求解方法选择
1200MPa下采用默认的牛顿-拉普生法老是遇到收敛问题，经常不收敛，或者单元高度扭曲。而且需要采用的载荷子步也要很多。（有人做的时候很少遇到这种情况，这与网格划分质量及网格多少有关，网格越少越容易收敛）
采用弧长法能解决这个问题。内径参数r=10+Δ，Δ=22mm 这种情况直接就求出来了。Δ参数不同时可能也不能直接就求出来，也不收敛。但是可以在450，550MPa的地方加两个载荷步就能收敛。如Δ=6mm时，直接用一个载荷步不收敛，我分成150，450，550，1200MPa几个载荷步(依次存为载荷步1，2，3，4)y ，再求解150到1200MPa就成功了（lssolve,1,4,1）。不收敛的原因在于500MPa是屈服极限。
后来验证过，网格画稀一点，不用弧长法直接用默认设置就可以收敛。
2、J积分解法的疑问
我是按照help上的介绍做的，但是对其求解的做法有些怀疑。个人认为求解δuy/δy偏导数的话，应该是把路径沿y方向移动才是对的。Help里面δux/δ x 的求法在数学上正确的。另外，ANSYS里面给了一个路径项求导的操作：general postproc->path operation ->differentiate 。（differentiate是不是求导，请指教） 那这个东西用来求偏导数不行吗？ 为什么help 里面要那么来求偏导数。
3、这个实例的建模
这个实例的建模，我是建的1/4模型。1/8模型也可以建出来，但是我对于1/8模型还能不能用对称边界条件有怀疑。1/4模型用对称边界条件是绝对正确的。另外，对称边界条件得到的约束条件在载荷步里面查看到，约束是发生在环向的（柱坐标系下），环向位移约束为零。加约束的时候直接加环向位移为零，也是可以的。
[ 本帖最后由 wishangtian 于 2008-4-24 15:23 编辑 ] 非常感谢，我看看在于你讨论
命令流部分
四、 命令流（log文件另附）1、 建模和求解部分（这里的建模网格划分比较密，可能不是很实用,这里的网格划分不好，在裂纹尖端第一行单元没有奇异性，最好还是用kscon 来做）：
/prep7
/COM,Structural
ET,1,PLANE2
KEYOPT,1,3,3
R,1,0.004,
MPTEMP,,,,,,,,
MPTEMP,1,0
MPDATA,EX,1,,2.1+011
MPDATA,PRXY,1,,0.3
MPDATA,DENS,1,,8500
TB,BISO,1,1,2,
TBTEMP,0
TBDATA,,5+008,6+009,,,,
R,1,0.004,
!*
wpstyle,0.001,0.1,-1,1,0.003,0,2,,5
k,1,0,0,0
circle,1,0.032,,,90
/PNUM,KP,1
/PNUM,LINE,1
/PNUM,AREA,1
/NUMBER,0
/REPLOT
k,4,0.034,0,0
k,5,0,0.034,0
k,6,0.042,0,0
k,7,0,0.046,0
k,8,0.060,0,0
k,9,0,0.060,0
k,10,0.080,0,0
k,11,0,0.080,0
k,12,0.12,0,0
k,13,0,0.12,0
k,14,0.18,0,0
k,15,0,0.18,0
k,16,0.32,0,0
k,17,0,0.32
circle,1,0.5,,,90
l,2,4
l,4,6
l,6,8
l,8,10
l,10,12
l,12,14
l,14,16
l,16,18
l,3,5
l,5,7
l,7,9
l,9,11
l,11,13
l,13,15
l,15,17
l,17,19
SAVE
al,all
/COLOR,NUM,DGRA,1
/COLOR,NUM,BMAG,2
/COLOR,NUM,RED,3
/COLOR,NUM,CBLU,4
/COLOR,NUM,MRED,5
/COLOR,NUM,GREE,6
/COLOR,NUM,ORAN,7
/COLOR,NUM,MAGE,8
/COLOR,NUM,YGRE,9
/COLOR,NUM,BLUE,10
/COLOR,NUM,GCYA,11
/REPLOT
aplot
lesize,1,0.0004
lesize,2,0.04
lesize,3,0.0002
lesize,4,0.0002
lesize,5,0.0004
lesize,6,0.001
lesize,7,0.002
lesize,8,0.004
lesize,9,0.008
lesize,10,0.016
lesize,11,0.0002
lesize,12,0.0002
lesize,13,0.0004
lesize,14,0.001
lesize,15,0.002
lesize,16,0.004
lesize,17,0.008
lesize,18,0.016
MSHAPE,1,2D
MSHKEY,0
CM,_Y,AREA
CMSEL,S,_Y
AMESH,_Y
CMDELE,_Y
save
/SOLU
NSUBST,5,50,2
AUTOTS,1
sfl,2,pres,-1.5+008,-1.5+008
allsel,all
SFTRAN
lsel,,,,4,10,1
lsel,a,,,12,18,1
lplot
dl,all,,symm
allsel,all
sbctran
lswrite,1
!*
NLGEOM,0
lsel,,,,2
sfl,2,pres,-1.2+009,-1.2+009,
allsel,all
sftran
NSUBST,15,1000,5
ARCLEN,1,0,0
lswrite,2
!*
allsel,all
save
!*
NLGEOM,0
sfl,2,pres,0.0,0.0,
allsel,all
sftran
OMEGA,0,0,62.831852,0
NSUBST,1,1,1
ARCLEN,0,0,0
lswrite,3
/GST,1
lssolve,1,2,1
2、 J积分部分（对于半边裂纹，如果你已经定义了路径的话，直接把这部分命令流输入进去就可以了）
/post1
!local,11,0,0.034,0,0!这里不应该建立局部坐标系。只有计算应力强度因子才需要。这里只需要保证全局坐标系的X方向与裂纹平行就是了。
csys,0
！这里应该有一个定义path，这里没有写出。
etable,volu,volu,
etable,sene,sene,
sexp,wden,sene,volu,1,-1,
pdef,wden,etab,wden,avg
pcalc,intg,wint,wden,yg
pcalc,intg,wint,wden,yg
pdef,sx,s,x,avg
pdef,sy,s,y,avg
pdef,sxy,s,xy,avg
pvect,norm,nx,ny,nz
pcalc,mult,sxnx,sx,nx
pcalc,mult,sxyxy,sxy,ny
pcalc,mult,syny,sy,ny
pcalc,mult,sxynx,sxy,nx
pcalc,add,tx,sxnx,sxyny
pcalc,add,ty,syny,sxynx
*get,dx,path,,last,s
pcalc,add,xg,xg,,,,-dx/200
pdef,ux1,u,x,avg
pdef,uy1,u,y,avg
pcalc,add,xg,xg,,,,dx/100
pdef,ux2,u,x,avg
pdef,uy2,u,y,avg
pcalc,add,xg,xg,,,,-dx/200
pcalc,add,pux,ux2,ux1,100/dx,-100/dx
pcalc,add,puy,uy2,uy1,100/dx,-100/dx
pcalc,mult,tpux,tx,pux
pcalc,mult,tpuy,ty,puy
pcalc,add,tpu,tpux,tpuy,
pcalc,intg,jtpu,tpu,s
pcalc,add,jint,wint,jtpu,1,-1
pcalc,add,jint,jint,jint
*get,jint,path,,last,jint
3、应力强度因子
（1）方法一、先建立局部坐标系：原点在裂纹尖端，x方向与裂纹平行，Y与裂纹垂直，笛卡尔坐标系。定义路径，直接点一下菜单路径就出来了，或者用kcalc就可以了。
（2）方法二、线弹性情况下。先算出J积分然后根据J积分与应力强度因子的关系来求应力强度因子。对于平面应变模型，J积分＝应力强度因子的平方×（1－泊松比×泊松比）/弹性模量
[ 本帖最后由 wishangtian 于 2006-7-19 14:26 编辑 ]
全尺寸裂纹模型
前面所建立的模型都只有裂纹的半边。为了验证模型的正确性，后面又建了一个全裂纹的模型。与前面模型的建立方式有一些不同：1、单元尺寸控制不使用lesize，而是在裂纹尖端用了一个KSCON命令建立concentrate keypoints。单元尺寸很粗糙。
2、模型的建立是在柱坐标系下进行，通过建立直线L实现的笛卡尔坐标下弧线的建立。
3、可能全尺寸裂纹模型的建立方式对大家有参考。主要是裂纹的上下表面在同一个位置，用不同的线/面来表示。
在命令流里面可以看到，直接用程序默认设置求解不收敛。把载荷步设一下就得到了求解结果。结果与前面的模型得到的结果接近。1200Mpa下的应力强度因子为0.55352E+10，J积分为4657362.7。说明：这个模型中积分路径是完全的。而前面的模型中是半边路径，计算中乘了2。
更新一下全尺寸模型的命令流（文件simwe_060715_002.rar,主要是J积分坐标系的问题）， 这里只给出一个局部变形图。
[ 本帖最后由 wishangtian 于 2006-7-15 08:29 编辑 ] 小弟，在此先谢了！ 我是按照help上的介绍做的，但是对其求解的做法有些怀疑。个人认为求解δuy/δy偏导数的话，应该是把路径沿y方向移动才是对的。Help里面δux/δ x 的求法在数学上正确的。
我觉得你是对的。因为书上的裂纹长度方向是x向，因此微分选项是dux/dx，如果裂纹方向旋转90度的话，当然是要对y微分了。 wishangtian 你好。我看了你做的例子。但是我对于裂纹分析还是有不明白的地方，尤其是积分路径的定义，希望能够得到你的帮助。因为自己一个人摸索是很困难的事情。能联系我吗？谢谢。我的qq是40864114
邮箱是lele@mail.nwpu.edu.cn希望能有您的回复哦。
三维J积分
前面提出的几个讨论问题通过验证已经解决了一些。通过与鬼脸嘟嘟 讨论后发现了一些新的问题以及错误。已经解决的问题：
1、收敛问题：网格画稀一点。
2、求解器与求解方法的选取：网格稀的情况下直接使用默认设置。网格比较密的情况下打开弧长法。对网格质量要求最低的求解器是Frontal Solver ,它适用于问题规模不大的非线性问题，求解耗费的时间多。
3、模型的建立：在模型具有对称性的情况下，可以建立部分模型。如1/4，1/8模型。
4、根据 鬼脸嘟嘟 的观点，按照help中讲的，实际上有两种求解断裂强度因子的方法。一种是在裂纹尖点的单元上取三个连续的点，然后直接kcalc,一种就是取包含裂纹尖点的路径，用J积分计算。在线弹性条件下，J积分有明确的意义，就是应变能释放率。因此才能够在线弹性阶段使用J积分来计算应力强度因子
尚待讨论的问题：
求解δuy/δy偏导数到底该不该按照help里面的介绍做。
更正前面的一个错误：
J积分不应该在我建立的那个局部坐标系下进行（附件里面的不方便改，把前面贴出来的命令流改了）。计算应力强度因子必须在那样的一个局部坐标系下进行。J积分应该在全局坐标系下进行。这个全局坐标系的X方向必须与裂纹平行。
——————————————————————————————————————————————————
更新：
1、讨论对于J积分所在坐标系的问题（讨论）：
最初我做J积分的时候是在计算应力强度因子所用的局部坐标系下进行的。后面发现help里面说要在全局坐标系下进行。认为在局部坐标系下进行是不对的。但是发发现两种坐标系下计算结果差不多。
现在想来，只要是与全局坐标系（要求x轴平行于裂纹）平行的局部笛卡儿坐标系，都能用来计算J积分。可以看到J积分的第二部分的所有量都与局部坐标系无关只与全局坐标系有关，而第一部分（也就是对应变能密度沿Y轴积分）与采用的局部坐标系的Y轴有关。由于J 积分具有积分路径无关性，我们选取一条关于裂纹对称的积分路径（积分路径是完全的）。那么，以这条路径上的应变能密度为纵轴，以Y轴为横轴画图，可以看到，得到的应变能密度——Y曲线是一条封闭曲线。对于一条封闭的曲线，对横坐标做积分，在坐标系平移的情况下，积分值不变。所以，J积分可以在于全局坐标系平行的坐标系下进行。
这只是一个不成熟的看法。希望清楚这个问题的朋友指正。
2、网格划分的问题
在参考书上和help中都要求靠近裂纹尖端的单元必须具有奇异性（也就是单元的中间节点要靠近裂纹尖端，这样才能更好地描述裂纹附近的应力场）。实际上，裂纹尖端附近单元不具有奇异性也无大碍。计算出来的J积分结果没有太大差别。
3、三维J积分
三维J积分的相对于二维J积分的难点在于：裂纹尖端附近单元没有KSCON这样的命令来直接生成，需要自己处理。Help里面的vm143例子有三维裂纹的详细命令流。
这里贴出help中三维J积分裂纹尖端附近单元生成的命令流，其余部分采用二维一样的做法。
下面这段命令流的使用说明：先划分网格，然后把裂纹尖端对应的那个节点定义为部件，名称为CRACKTIP，然后输入如下部分命令流，裂纹尖端附近第一行单元便调整成奇异单元。之后的加载和J积分部分无特别之处。
附件中是整理过的vm143中的三维裂纹J积分部分。
/NOPR
NSEL,ALL
*GET,N,NODE,,NUM,MAX ! CURRENT MAXIMUM NODE NUMBER
CMSEL,S,CRACKTIP ! SELECT THE TIP NODES
ESLN ! ANY ELEMENTS ATTACHED
*GET,ELMAX,ELEM,,NUM,MAX ! CURRENT MAXIMUM ELEMENT NUMBER
*DO,IEL,1,ELMAX ! LOOP ON MAX ELEMENT
ELMI=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 ELEMENT
N3 = NELEM(ELMI,3) ! GET NODE 3 (K)
*IF,NSEL(N3),LE,0,CYCLE ! IT MUST BE SELECTED
N7 = NELEM(ELMI,7) ! GET NODE 7 (L)
*IF,NSEL(N7),LE,0,CYCLE ! IT MUST ALSO BE SELECTED
N1 = 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 N3
Y3 = 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) + Y3
Z= 0.25*NZ(N2) + Z3
N= N + 1 ! NEXT NODE
N10 = N
N,N10,X,Y,Z ! MIDSIDE NODE LOCATION
X= 0.25*NX(N1) + X3
Y= 0.25*NY(N1) + Y3
Z= 0.25*NZ(N1) + Z3
N= N + 1
N12= N
N,N12,X,Y,Z
X7 = 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 + 1
N14 = N
N,N14,X,Y,Z
X= 0.25*NX(N5) + X7
Y= 0.25*NY(N5) + Y7
Z= 0.25*NZ(N5) + Z7
N= N + 1
N16 = N
N,N16,X,Y,Z
N4=N3
N8=N7
NSEL,ALL
TYPE,3
EN,ELMI,N1,N2,N3,N4,N5,N6,N7,N8 ! REDEFINE THE ELEMENT
EMORE,0,N10,0,N12,0,N14,0,N16
EMORE,
*ENDDO
CMSEL,U,CRACKTIP ! UNSELECT THE TIP NODES
NUMMRG,NODE ! MERGE MIDSIDE NODES
NSEL,ALL ! SELECT ALL ELEMENTS
ESEL,ALL ! SELECT ALL ELEMENTS
/GOPR
*END
[ 本帖最后由 wishangtian 于 2006-7-24 15:40 编辑 ]
全尺寸裂纹模型J积分路径的选取
由于全尺寸裂纹模型中，裂纹上下面在一块，同一个位置可能有两个点，造成J积分路径选取困难（如图1）。一个处理办法就是：求解后，画出模型的变形图（画应力图，位移图或者是应变图都行）。在变形图上选取积分路径就比较轻松。如图（如图2）。这样就得到想要的路径了。后面先plot以下单元，再plot path就如下图（图3）所示了。 我要找的原来在这里，谢谢搂住 多谢楼主分享 按照断裂力学的基本理论，积分路径应该是定义在弹性区就可以了。对吗？ 高手，谢谢了！！！！！！！！！！！！！！
求教
此模型可以用来进行模态分析么？在裂纹位置，可否检测出损伤的存在？哈哈，不错！！
不错不错不错不错不错不错不错不错不错不错!!! 原帖由 bewinner01 于 2007-2-5 14:48 发表此模型可以用来进行模态分析么？在裂纹位置，可否检测出损伤的存在？
你试试看吧。
我对振动、模态方法的东西没有兴趣也没有任何了解。
我算出的KI很小 KI = 0.11344E+06
使用对称模型,1200MPa情况下 ,你算出的"1200Mpa下的应力强度因子为0.55352E+10"
而我的是
****KI =0.11344E+06, KII = 0.0000 , KIII = 0.0000 ****