1、分享一个断裂分析实例（J 积分和应力强度因子）（最终版）有些地方没有加注释。 加起来比较麻烦。 后面不再修改了。建议参考 help 关于 fracture 部分以及 vm143 部分。有更新attach92985/attach以上附件没有删除。一、 问题描述如图所示为一等厚度空心圆盘（参看图 1. 1）， 厚度 4mm， 内径 r=10+ ， =22mm，外径 R=500mm， 材料属性数据： 双线性（参看图 1. 3） E=2. 11011Pa， ET=6. 0109Pa， 屈服极限 s=500MPa， =0. 3， 密度 =8500kg/m3， 采用 Mises 屈服准则， 裂纹初始长度为
2、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 比较。attach93106/attach attach93107/attachattach93108/attach二、 求解求解中统一采用国际单位制， 长度 m， 压力、 应力与弹性模量 Pa， 密度
3、 Kg/m3， 转速rad/s。对圆盘的 1/4 进行 ANSYS 建模， 网格划分如图 2. 1。 单元类型为 6 节点三角形单元 plane2。裂纹附近单元边长为 0. 0002m。 载荷施加如图 2. 2， 扇形两条半径（裂纹处除外） 上施加对称位移边界条件， 弧上加均布拉力。 裂纹处无位移约束。attach93109/attachattach93110/attach计算的应力强度因子和 J 积分结果如表 2. 1。 前面给出的 J 积分由于坐标系错误， 做法不对， 现在改正 1200Mpa 下的 J 积分结果。 但是， 结果差别不大。 因为， 即使在局部坐标系下， J 积分中用到的 X
4、G, YG, ZG 的坐标还是在全球坐标系下的。 坐标系的改变只对 J 积分的第一项有一点影响（对 Y 坐标积分这部分， Y 坐标变了）。新加两个附件： 改正后求解 1200MPa 下 J 积分的命令流， 及两张路径定义的图片命令流里面涉及到路径定义的命令是 GUI 方式进行的。所以要分开看。attach93704/attach。表 2. 1 应力强度因子和 J 积分结果载荷 项目 不同路径下的结果150MPa KI 3. 1359e7 2. 986e73. 3068e7 2. 7695e7（此处 J 积分没改） J-Integral 3193. 41629 3182. 07917 3245.
5、 220423167. 233521200MPa KI 5. 4207e9 5. 3398e95. 382e9 5. 285e9J-Integral 4701408. 68 4777189. 41 4674756. 714709657. 93600r/minKI550115490154672 54463（此处 J 积分没改） J-Integral 1. 5649802 1. 5719383 1. 534057481. 57113766三、 一些需要讨论的问题（后面补的帖解决了大部分问题， 这些留在这里供参考）：1、 求解方法选择1200MPa 下采用默认的牛顿-拉普生法老是遇到收敛问题， 经常
6、不收敛， 或者单元高度扭曲。 而且需要采用的载荷子步也要很多。（有人做的时候很少遇到这种情况， 这与网格划分质量及网格多少有关， 网格越少越容易收敛）采用弧长法能解决这个问题。 内径参数 r=10+ ， =22mm 这种情况直接就求出来了。 参数不同时可能也不能直接就求出来， 也不收敛。 但是可以在 450， 550MPa 的地方加两个载荷步就能收敛。 如 =6mm 时， 直接用一个载荷步不收敛， 我分成 150， 450， 550， 1200MPa几个载荷步(依次存为载荷步 1，2， 3， 4) y ， 再求解 150 到 1200MPa 就成功了（lssolve, 1, 4, 1）。不收敛
7、的原因在于 500MPa 是屈服极限。后来验证过， 网格画稀一点， 不用弧长法直接用默认设置就可以收敛。2、 J 积分解法的疑问我是按照 help 上的介绍做的， 但是对其求解的做法有些怀疑。 个人认为求解 uy/ y 偏导数的话， 应该是把路径沿 y 方向移动才是对的。 Help 里面 ux/ x 的求法在数学上正确的。 另外， ANSYS 里面给了一个路径项求导的操作： general postproc-pathoperation -differentiate 。（differentiate 是不是求导， 请指教） 那这个东西用来求偏导数不行吗？ 为什么 help 里面要那么来求偏导数。3
8、、 这个实例的建模这个实例的建模， 我是建的 1/4 模型。 1/8 模型也可以建出来， 但是我对于 1/8 模型还能不能用对称边界条件有怀疑。 1/4 模型用对称边界条件是绝对正确的。 另外， 对称边界条件得到的约束条件在载荷步里面查看到， 约束是发生在环向的（柱坐标系下）， 环向位移约束为零。 加约束的时候直接加环向位移为零， 也是可以的。 attach92985/attach命令流部分四、 命令流（log 文件另附） attach93111/attach1、 建模和求解部分（这里的建模网格划分比较密， 可能不是很实用, 这里的网格划分不好， 在裂纹尖端第一行单元没有奇异性， 最好还是用
9、kscon 来做）：/prep7/COM, StructuralET, 1, PLANE2KEYOPT, 1, 3, 3R, 1, 0. 004,MPTEMP, , , , , , , ,MPTEMP, 1, 0MPDATA, EX, 1, , 2. 1+011MPDATA, PRXY, 1, , 0. 3MPDATA, DENS, 1, , 8500TB, BISO, 1, 1, 2,TBTEMP, 0TBDATA, , 5+008, 6+009, , , ,R, 1, 0. 004,!*wpstyle, 0. 001, 0. 1, -1, 1, 0. 003, 0, 2, , 5k, 1,
10、 0, 0, 0circle, 1, 0. 032, , , 90/PNUM, KP, 1/PNUM, LINE, 1/PNUM, AREA, 1/NUMBER, 0/REPLOTk, 4, 0. 034, 0, 0k, 5, 0, 0. 034, 0k, 6, 0. 042, 0, 0k, 7, 0, 0. 046, 0k, 8, 0. 060, 0, 0k, 9, 0, 0. 060, 0k, 10, 0. 080, 0, 0k, 11, 0, 0. 080, 0k, 12, 0. 12, 0, 0k, 13, 0, 0. 12, 0k, 14, 0. 18, 0, 0k, 15, 0,
11、0. 18, 0k, 16, 0. 32, 0, 0k, 17, 0, 0. 32circle, 1, 0. 5, , , 90l, 2, 4l, 4, 6l, 6, 8l, 8, 10l, 10, 12l, 12, 14l, 14, 16l, 16, 18l, 3, 5l, 5, 7l, 7, 9l, 9, 11l, 11, 13l, 13, 15l, 15, 17l, 17, 19SAVEal, all/COLOR, NUM, DGRA, 1/COLOR, NUM, BMAG, 2/COLOR, NUM, RED, 3/COLOR, NUM, CBLU, 4/COLOR, NUM, MRE
12、D, 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/REPLOTaplotlesize, 1, 0. 0004lesize, 2, 0. 04lesize, 3, 0. 0002lesize, 4, 0. 0002lesize, 5, 0. 0004lesize, 6, 0. 001lesize, 7, 0. 002lesize, 8, 0. 004lesize, 9, 0. 008lesize
13、, 10, 0. 016lesize, 11, 0. 0002lesize, 12, 0. 0002lesize, 13, 0. 0004lesize, 14, 0. 001lesize, 15, 0. 002lesize, 16, 0. 004lesize, 17, 0. 008lesize, 18, 0. 016MSHAPE, 1, 2DMSHKEY, 0CM, _Y, AREACMSEL, S, _YAMESH, _YCMDELE, _Ysave/SOLUNSUBST, 5, 50, 2AUTOTS, 1sfl, 2, pres, -1. 5+008, -1. 5+008allsel,
14、allSFTRANlsel, , , , 4, 10, 1lsel, a, , , 12, 18, 1lplotdl, all, , symmallsel, allsbctranlswrite, 1!*NLGEOM, 0lsel, , , , 2sfl, 2, pres, -1. 2+009, -1. 2+009,allsel, allsftranNSUBST, 15, 1000, 5ARCLEN, 1, 0, 0lswrite, 2!*allsel, allsave!*NLGEOM, 0sfl, 2, pres, 0. 0, 0. 0,allsel, allsftranOMEGA, 0, 0
15、, 62. 831852, 0NSUBST, 1, 1, 1ARCLEN, 0, 0, 0lswrite, 3/GST, 1lssolve, 1, 2, 12、 J 积分部分（对于半边裂纹， 如果你已经定义了路径的话， 直接把这部分命令流输入进去就可以了）/post1!local, 11, 0, 0. 034, 0, 0 !这里不应该建立局部坐标系。 只有计算应力强度因子才需要。 这里只需要保证全局坐标系的 X 方向与裂纹平行就是了。csys, 0！ 这里应该有一个定义 path， 这里没有写出。etable, volu, volu,etable, sene, sene,sexp, wden,
16、 sene, volu, 1, -1,pdef, wden, etab, wden, avgpcalc, intg, wint, wden, ygpcalc, intg, wint, wden, ygpdef, sx, s, x, avgpdef, sy, s, y, avgpdef, sxy, s, xy, avgpvect, norm, nx, ny, nzpcalc, mult, sxnx, sx, nxpcalc, mult, sxyxy, sxy, nypcalc, mult, syny, sy, nypcalc, mult, sxynx, sxy, nxpcalc, add, tx
17、, sxnx, sxynypcalc, add, ty, syny, sxynx*get, dx, path, , last, spcalc, add, xg, xg, , , , -dx/200pdef, ux1, u, x, avgpdef, uy1, u, y, avgpcalc, add, xg, xg, , , , dx/100pdef, ux2, u, x, avgpdef, uy2, u, y, avgpcalc, add, xg, xg, , , , -dx/200pcalc, add, pux, ux2, ux1, 100/dx, -100/dxpcalc, add, puy
18、, uy2, uy1, 100/dx, -100/dxpcalc, mult, tpux, tx, puxpcalc, mult, tpuy, ty, puypcalc, add, tpu, tpux, tpuy,pcalc, intg, jtpu, tpu, spcalc, add, jint, wint, jtpu, 1, -1pcalc, add, jint, jint, jint*get, jint, path, , last, jint3、 应力强度因子（1） 方法一、 先建立局部坐标系： 原点在裂纹尖端， x 方向与裂纹平行， Y 与裂纹垂直，笛卡尔坐标系。 定义路径， 直接点一下
19、菜单路径就出来了， 或者用 kcalc 就可以了。attach93111/attach（2） 方法二、 线弹性情况下。 先算出 J 积分然后根据 J 积分与应力强度因子的关系来求应力强度因子。 对于平面应变模型， J 积分应力强度因子的平方（1泊松比泊松比）/弹性模量全尺寸裂纹模型前面所建立的模型都只有裂纹的半边。 为了验证模型的正确性， 后面又建了一个全裂纹的模型。 与前面模型的建立方式有一些不同：1、 单元尺寸控制不使用 lesize， 而是在裂纹尖端用了一个 KSCON 命令建立 concentratekeypoints。 单元尺寸很粗糙。2、 模型的建立是在柱坐标系下进行， 通过建立直
20、线 L 实现的笛卡尔坐标下弧线的建立。3、 可能全尺寸裂纹模型的建立方式对大家有参考。 主要是裂纹的上下表面在同一个位置， 用不同的线/面来表示。在命令流里面可以看到， 直接用程序默认设置求解不收敛。 把载荷步设一下就得到了求解结果。 结果与前面的模型得到的结果接近。 1200Mpa 下的应力强度因子为 0. 55352E+10，J 积分为 4657362. 7。 说明： 这个模型中积分路径是完全的。 而前面的模型中是半边路径， 计算中乘了 2。更新一下全尺寸模型的命令流（文件 simwe_060715_002. rar, 主要是 J 积分坐标系的问题）， attach93705/attach
21、 这里只给出一个局部变形图。 attach93340/attach三维 J 积分前面提出的几个讨论问题通过验证已经解决了一些。通过与 鬼脸嘟嘟 讨论后发现了一些新的问题以及错误。已经解决的问题：1、 收敛问题： 网格画稀一点。2、 求解器与求解方法的选取： 网格稀的情况下直接使用默认设置。 网格比较密的情况下打开弧长法。 对网格质量要求最低的求解器是 Frontal Solver , 它适用于问题规模不大的非线性问题， 求解耗费的时间多。3、 模型的建立： 在模型具有对称性的情况下， 可以建立部分模型。 如 1/4， 1/8 模型。4、 根据 鬼脸嘟嘟 的观点， 按照 help 中讲的， 实际
22、上有两种求解断裂强度因子的方法。 一种是在裂纹尖点的单元上取三个连续的点， 然后直接 kcalc, 一种就是取包含裂纹尖点的路径， 用 J 积分计算。 在线弹性条件下， J 积分有明确的意义， 就是应变能释放率。 因此才能够在线弹性阶段使用 J 积分来计算应力强度因子尚待讨论的问题：求解 uy/ y 偏导数到底该不该按照 help 里面的介绍做。更正前面的一个错误：J 积分不应该在我建立的那个局部坐标系下进行（附件里面的不方便改， 把前面贴出来的命令流改了）。 计算应力强度因子必须在那样的一个局部坐标系下进行。 J 积分应该在全局坐标系下进行。 这个全局坐标系的 X 方向必须与裂纹平行。更新：
23、1、 讨论对于 J 积分所在坐标系的问题（讨论）：最初我做 J 积分的时候是在计算应力强度因子所用的局部坐标系下进行的。 后面发现 help 里面说要在全局坐标系下进行。 认为在局部坐标系下进行是不对的。 但是发发现两种坐标系下计算结果差不多。现在想来， 只要是与全局坐标系（要求 x 轴平行于裂纹） 平行的局部笛卡儿坐标系，都能用来计算 J 积分。 可以看到 J 积分的第二部分的所有量都与局部坐标系无关只与全局坐标系有关， 而第一部分（也就是对应变能密度沿 Y 轴积分） 与采用的局部坐标系的 Y 轴有关。由于 J 积分具有积分路径无关性， 我们选取一条关于裂纹对称的积分路径（积分路径是完全的）
24、。 那么， 以这条路径上的应变能密度为纵轴， 以 Y 轴为横轴画图， 可以看到， 得到的应变能密度Y 曲线是一条封闭曲线。 对于一条封闭的曲线， 对横坐标做积分， 在坐标系平移的情况下， 积分值不变。 所以， J 积分可以在于全局坐标系平行的坐标系下进行。这只是一个不成熟的看法。 希望清楚这个问题的朋友指正。2、 网格划分的问题在参考书上和 help 中都要求靠近裂纹尖端的单元必须具有奇异性（也就是单元的中间节点要靠近裂纹尖端， 这样才能更好地描述裂纹附近的应力场）。 实际上， 裂纹尖端附近单元不具有奇异性也无大碍。 计算出来的 J 积分结果没有太大差别。3、 三维 J 积分三维 J 积分的相
25、对于二维 J 积分的难点在于： 裂纹尖端附近单元没有 KSCON 这样的命令来直接生成， 需要自己处理。 Help 里面的 vm143 例子有三维裂纹的详细命令流。这里贴出 help中三维J积分裂纹尖端附近单元生成的命令流， 其余部分采用二维一样的做法。下面这段命令流的使用说明： 先划分网格， 然后把裂纹尖端对应的那个节点定义为部件， 名称为 CRACKTIP， 然后输入如下部分命令流， 裂纹尖端附近第一行单元便调整成奇异单元。 之后的加载和 J 积分部分无特别之处。附件中是整理过的 vm143 中的三维裂纹 J 积分部分。 attach94253/attach/NOPRNSEL, ALL*G
26、ET, 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, A
27、TTR, 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
28、 = 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) + Y3 Z = 0. 25*NZ(N2) + Z3N = N + 1 ! NEX
29、T 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) + X7Y = 0. 25*NY(N6) + Y7Z = 0. 25*NZ(N6) + Z7N = N + 1N14 = NN, N14, X, Y, ZX = 0. 25*
30、NX(N5) + X7Y = 0. 25*NY(N5) + Y7Z = 0. 25*NZ(N5) + Z7N = 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 NODESNUMMRG, NODE ! MERGE MIDSIDE NODESNSEL, ALL ! SELECT ALL ELEMENTSESEL, ALL ! SELECT ALL ELEMENTS/GOPR*END由于全尺寸裂纹模型中， 裂纹上下面在一块， 同一个位置可能有两个点， 造成 J 积分路径选取困难（如图 1）。 一个处理办法就是： 求解后， 画出模型的变形图（画应力图， 位移图或者是应变图都行）。 在变形图上选
