ANSYS命令流编写三维问题应力强度因子的计算程序

/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 = %26quot;KI1%26quot;

　　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