【分享】计算三维问题的应力强度因子的程序

============================

/COM,ANSYS MEDIA REL. 60 (090601) REF. VERIF. MANUAL: REL. 60
/VERIFY,VM143
*CREATE,FRACT,MAC
/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

/PREP7
*afun,deg
InnerRadius=0.1 !InnerRadius为裂纹半径
OuterRadius=1 !OuterRadius为圆柱半径
Scaler=0.025 !Scaler为裂纹前沿单元范围，已经是最佳
BaseHeight=0.51 !BaseHeight为基层高度,it may be the best
LayerHeight=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 PLATE
C*** BROWN AND SRAWLEY, ASTM SPECIAL TECHNICAL PUBLICATION NO. 410.
/COM, ****** CRACK IN 3-DIMENSIONS USING SOLID45 AND SOLID95
ANTYPE,STATIC ! STATIC ANALYSIS
ET,1,SOLID45
ET,2,SOLID45 ! ELEMENTS AROUND THE CRACK TIP
ET,3,SOLID95 ! CRACK TIP ELEMENTS CREATED USING MACRO FRACT
MP,EX,1,2e4
MP,NUXY,1,.3 ! CYLINDRICAL COORDINATE SYSTEM
local,33,1,InnerRadius
csys,33
N,1
NGEN,9,20,1
N,11,Scaler
N,171,Scaler,180
FILL,11,171,7,31,20
local,33,0,InnerRadius
csys,33
FILL,1,11,9,2,1,9,20,3
N,15,OuterRadius-InnerRadius
N,75,OuterRadius-InnerRadius,BaseHeight
FILL,15,75,2,35,20
N,155,-InnerRadius,BaseHeight
NGEN,2,200,155 !155点在轴线上，故而在此与355为同一点，目的在于产生轴线上的单元
FILL,75,155,3,95,20
N,172,-InnerRadius
NGEN,2,200,172 !且因为无论什么坐标对此无影响，所以可以直接产生，与155点类似
FILL,155,172,5,177,-1,,,.15
ngen,2,200,173,177,1 !将轴线上的节点依次产生其同一位置的节点，便于产生轴线单元
csys,33
FILL,11,15,3,,,7,20,3
csys,5
ngen,2,200,1,177,,,RotationAngle !在此只是产生了一次的节点，并不是21次（180度）
csys,0
E,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,3
EMODIF,1 ! MODIFY ELEMENTS 1 TO 8 FROM TYPE,1 TO TYPE,2
*REPEAT,8,1
NUMMRG,NODE ! MERGE COINCIDENT NODES
csys,33
NSEL,S,LOC,X,0
NSEL,R,LOC,Y,0
CM,CRACKTIP,NODE
/NERR,0 ! TEMPORARILY NO WARNINGS OR ERRORS PRINTOUT
                            ! (IN ORDER TO AVOID WARNING MESSAGES DUE TO
                            ! MIDSIDE NODES LOCATION)
FRACT,2 ! CONVERSION MACRO, TYPE 2 IS SOLID45
                            ! ELEMENTS AROUND THE CRACK TIP

csys,0
type,1
e,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,-1
e,51,52,72,71,251,252,272,271 !基本同上，为第三层，第一列
egen,4,1,-1
e,71,72,92,91,271,272,292,291 !基本同上，为第四层，第一列
egen,4,1,-1
e,91,92,112,111,291,292,312,311 !基本同上，为第五层，第一列
egen,4,1,-1
e,111,112,132,131,311,312,332,331
egen,4,1,-1
csys,5
egen,Rotationtimes,200,1,110,1,,,,,,,RotationAngle
csys,0
nsel,s,loc,y,BaseHeight
ngen,2,5000,all,,,,LayerHeight
e,135,335,155,155,5135,5335,5155,5155
csys,5
egen,Rotationtimes,200,-1
e,135,115,315,335,5135,5115,5315,5335
egen,Rotationtimes,200,-1
e,115,95,295,315,5115,5095,5295,5315
egen,Rotationtimes,200,-1
e,95,75,275,295,5095,5075,5275,5295
egen,Rotationtimes,200,-1

csys,0
nsel,s,loc,y,BaseHeight+LayerHeight
esln,s,0,all
egen,LayerAmount,5000,all,,,,,,,,,LayerHeight
nsel,all
esel,all
/NERR,DEFA ! TURN ON THE WARNINGS OR ERRORS PRINTOUT
/OUTPUT
OUTPR,,ALL
csys,0
NSEL,S,LOC,z,0
DSYM,SYMM,z ! SYMMETRIC B.C.'S AT X = 0
d,all,uz
nsel,all
nsel,s,loc,y,0
csys,5
nsel,r,loc,x,InnerRadius,OuterRadius
DSYM,SYMM,Y ! SYMMETRIC B.C.'S AT Y = 0 EXCEPT CRACK NODES
d,all,uy
nsel,all
csys,0
nsel,s,loc,x,0
dsym,symm,x
d,all,ux
NSEL,S,LOC,Y,Height
SF,ALL,PRES,-1000
NSEL,ALL
ESEL,ALL
FINISH
!now is the solution
/OUTPUT,SCRATCH
/SOLU
SOLVE
FINISH
/OUTPUT
/POST1
C*** IN POST1 DETERMINE KI (STRESS INTENSITY FACTOR) USING KCALC !**
csys,33 !define the local systerm
PATH,KI1,3,,48 ! DEFINE PATH WITH NAME = "KI1"
PPATH,1,1 ! DEFINE PATH POINTS BY NODE
PPATH,2,170
PPATH,3,171
KCALC,,,1 ! COMPUTE KI FOR A HALF-MODEL WITH SYMM. B.C.
*GET,KI1,KCALC,,K,1 ! GET KI AS PARAMETER KI1