一、 問題描述
如圖所示為一等厚度空心圓盤（參看圖1.1），厚度4mm，內(nèi)徑r=10+Δ，Δ=22mm，外徑R=500mm，材料屬性數(shù)據(jù)：雙線性（參看圖1.3）E=2.1×1011Pa，ET=6.0×109Pa，屈服極限σs=500MPa，μ=0.3，密度ρ=8500kg/m3，采用Mises屈服準(zhǔn)則，裂紋初始長度為a0=2.0mm。裂紋如圖1.2，裂紋位于0度90度180度和270度的位置，圖中粗短線表示裂紋。
要求用有限元解
（1） 載荷為均布拉力q=150MPa時的KI。（參看圖1.1）
（2） 載荷為均布拉力q=1200MPa時的J積分值。（參看圖1.1）
（3） 當(dāng)載荷為轉(zhuǎn)速n=600r/min的KI ，并計算沿裂紋尖端不同路徑的積分值，與KI比較。



二、 求解
求解中統(tǒng)一采用國際單位制，長度m，壓力、應(yīng)力與彈性模量Pa，密度Kg/m3，轉(zhuǎn)速rad/s。
對圓盤的1/4進(jìn)行ANSYS建模，網(wǎng)格劃分如圖2.1。單元類型為6節(jié)點三角形單元plane2。裂紋附近單元邊長為0.0002m。載荷施加如圖2.2，扇形兩條半徑（裂紋處除外）上施加對稱位移邊界條件，弧上加均布拉力。裂紋處無位移約束。


計算的應(yīng)力強(qiáng)度因子和J積分結(jié)果如表2.1。前面給出的J積分由于坐標(biāo)系錯誤，做法不對，現(xiàn)在改正1200Mpa下的J積分結(jié)果。但是，結(jié)果差別不大。因為，即使在局部坐標(biāo)系下，J積分中用到的 XG, YG, ZG的坐標(biāo)還是在全球坐標(biāo)系下的。坐標(biāo)系的改變只對J積分的第一項有一點影響（對Y坐標(biāo)積分這部分，Y坐標(biāo)變了）。
新加兩個附件：改正后求解1200MPa下J積分的命令流，及兩張路徑定義的圖片命令流里面涉及到路徑定義的命令是GUI方式進(jìn)行的。所以要分開看。
。
表2.1 應(yīng)力強(qiáng)度因子和J積分結(jié)果
載荷 項目 不同路徑下的結(jié)果
150MPa KI 3.1359e7 2.986e7 3.3068e7 2.7695e7
（此處J積分沒改） J-Integral 3193.41629 3182.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
三、一些需要討論的問題（后面補(bǔ)的帖解決了大部分問題，這些留在這里供參考）：
1、求解方法選擇
1200MPa下采用默認(rèn)的牛頓-拉普生法老是遇到收斂問題，經(jīng)常不收斂，或者單元高度扭曲。而且需要采用的載荷子步也要很多。（有人做的時候很少遇到這種情況，這與網(wǎng)格劃分質(zhì)量及網(wǎng)格多少有關(guān)，網(wǎng)格越少越容易收斂）
采用弧長法能解決這個問題。內(nèi)徑參數(shù)r=10+Δ，Δ=22mm 這種情況直接就求出來了。Δ參數(shù)不同時可能也不能直接就求出來，也不收斂。但是可以在450，550MPa的地方加兩個載荷步就能收斂。如Δ=6mm時，直接用一個載荷步不收斂，我分成150，450，550，1200MPa幾個載荷步(依次存為載荷步1，2，3，4)y ，再求解150到1200MPa就成功了（lssolve,1,4,1）。不收斂的原因在于500MPa是屈服極限。
后來驗證過，網(wǎng)格畫稀一點，不用弧長法直接用默認(rèn)設(shè)置就可以收斂。
2、J積分解法的疑問
我是按照help上的介紹做的，但是對其求解的做法有些懷疑。個人認(rèn)為求解δuy/δy偏導(dǎo)數(shù)的話，應(yīng)該是把路徑沿y方向移動才是對的。Help里面δux/δ x 的求法在數(shù)學(xué)上正確的。另外，ANSYS里面給了一個路徑項求導(dǎo)的操作：general postproc->path operation ->differentiate 。（differentiate是不是求導(dǎo)，請指教） 那這個東西用來求偏導(dǎo)數(shù)不行嗎？ 為什么help 里面要那么來求偏導(dǎo)數(shù)。
3、這個實例的建模
這個實例的建模，我是建的1/4模型。1/8模型也可以建出來，但是我對于1/8模型還能不能用對稱邊界條件有懷疑。1/4模型用對稱邊界條件是絕對正確的。另外，對稱邊界條件得到的約束條件在載荷步里面查看到，約束是發(fā)生在環(huán)向的（柱坐標(biāo)系下），環(huán)向位移約束為零。加約束的時候直接加環(huán)向位移為零，也是可以的。
命令流部分四、 命令流（log文件另附）
1、 建模和求解部分（這里的建模網(wǎng)格劃分比較密，可能不是很實用,這里的網(wǎng)格劃分不好，在裂紋尖端第一行單元沒有奇異性，最好還是用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積分部分（對于半邊裂紋，如果你已經(jīng)定義了路徑的話，直接把這部分命令流輸入進(jìn)去就可以了）
/post1
!local,11,0,0.034,0,0 !這里不應(yīng)該建立局部坐標(biāo)系。只有計算應(yīng)力強(qiáng)度因子才需要。這里只需要保證全局坐標(biāo)系的X方向與裂紋平行就是了。
csys,0
！這里應(yīng)該有一個定義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、應(yīng)力強(qiáng)度因子
（1）方法一、先建立局部坐標(biāo)系：原點在裂紋尖端，x方向與裂紋平行，Y與裂紋垂直，笛卡爾坐標(biāo)系。定義路徑，直接點一下菜單路徑就出來了，或者用kcalc就可以了。
（2）方法二、線彈性情況下。先算出J積分然后根據(jù)J積分與應(yīng)力強(qiáng)度因子的關(guān)系來求應(yīng)力強(qiáng)度因子。對于 平面應(yīng)變模型，J積分＝應(yīng)力強(qiáng)度因子的平方×（1－泊松比×泊松比）/彈性模量
全尺寸裂紋模型前面所建立的模型都只有裂紋的半邊。為了驗證模型的正確性，后面又建了一個全裂紋的模型。與前面模型的建立方式有一些不同：
1、單元尺寸控制不使用lesize，而是在裂紋尖端用了一個KSCON命令建立concentrate keypoints。單元尺寸很粗糙。
2、模型的建立是在柱坐標(biāo)系下進(jìn)行，通過建立直線L實現(xiàn)的笛卡爾坐標(biāo)下弧線的建立。
3、可能全尺寸裂紋模型的建立方式對大家有參考。主要是裂紋的上下表面在同一個位置，用不同的線/面來表示。
在命令流里面可以看到，直接用程序默認(rèn)設(shè)置求解不收斂。把載荷步設(shè)一下就得到了求解結(jié)果。結(jié)果與前面的模型得到的結(jié)果接近。1200Mpa下的應(yīng)力強(qiáng)度因子為 0.55352E+10，J積分為4657362.7。說明：這個模型中積分路徑是完全的。而前面的模型中是半邊路徑，計算中乘了2。
更新一下全尺寸模型的命令流（文件caenet_060715_002.rar,主要是J積分坐標(biāo)系的問題）， 這里只給出一個局部變形圖。

locdeform.jpg