为什么我的模型计算出的应力强度因子与裂缝扩展的曲线

gradywzy

为什么我的模型计算出的弯曲受力的应力强度因子与裂缝扩展的曲线是增长型的   而不是抛物线型的呢？以下是我的命令流（静力分析）小弟是新手  希望各位前辈  帮我指点指点

!Fracture mechanistic ,dual wheel
*CREATE,FRACT,MAC
!  MACRO TO CREATE 3D SOLID95 CRACK TIP ELEMENTS FROM 3D SOLID45 ELEMENTS
!     MAKE A COMPONENT CONTAINING THE CRACK TIP NODES (CRACKTIP)
!     THE CRACK TIP IS BETWEEN NODES K AND O
!     SET ELEMENT TYPE TO POINT TO SOLID95
!     SET ARG1 TO  N  (THE TYPE OF THE ELEMENTS AROUND THE CRACK TIP)
!
/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             !all units are mm and MPa
*set,if_bond,1     ! the bonding condition between AC and PCC/old AC? if bonding fully, the value is 1; if separate,the value is 0
*set,if_has_c,1    !!!if the PCC or old AC has crack? if No, the value is 0; if Yes,the value is 1
*set,if_c_wid,1    !!!if the PCC or old AC has crack,then, does it has width? if No, the value is 0; if Yes,the value is 1
                   !!!Attention!if the if_has_c is 0, the if_c_wid must be 0.(namely, when the PCC has not crack,it must is 0)
*set,crack_w,1     !the width of PCC-crack for PCC or old AC layer
*set,crack_L,140  !the length of crack from the AC layer bottom to crack tip
*SET,N_AC_T,7      !!the division amount from crack tip to road surface, it should be changed according to the crack length and AC thickness
*SET,N_AC_B,1      !!the division amount from crack tip to road bottom,  it should be changed according to the crack length and AC thickness   
*set,AC_T,150     !the thickness of the AC layer
*set,PCC_T,150      !the thickness of the PCC or old AC layer
*set,BL_T,150     !the thickness of the old base layer
*set,SUB_T,5000    !the thickness of the subgrade
*SET,M1,1200    !the modulus of AC/MPa
*SET,M2,7000        !the modulus of PCC or old AC
*SET,M3,7000      !the modulus of Base
*SET,M4,50         !the modulus of subgrade
*set,P,0.689        !the load/MPa
*set,road_W,4572   !the half width of calculated road
*set,road_L,3000   !the half length of calculated road
*set,LOAD_L,198    !the length of the load contact area (along driving direction)
*set,LOAD_W,136    !the width of the load contact area
*set,LOAD_S,343    !the space of dual wheel
*set,LOAD_D,1828.8   !the distance between dual wheel
*set,a,5           !the amount of crack element belong to half tire contact area along road width direction
*set,b,10          !the amount of crack element for tire width
*set,c,10          !the amount of crack element between dual wheel (the outside edge and inside edge)
*set,d,10          !the amount of crack element between the outside edge of tire and road edge, it will not be used in this program
*set,N_ratio,15    !the ratio between last and first section, when NGEN the node and DRAG the volume from tire edge to road edge
*SET,N_Load,5      !!the division amount for part of load area length (along the driving direction,the other part is for crack element)
*SET,N_PCC,2       !the division amount of PCC or old AC                  
*SET,N_BL,2        !the division amount of Base                  
*SET,N_SUB,5       !the division amount of subgrade
*SET,U1,0.35       !the Poson's Ratio of AC                  
*SET,U2,0.20       !the Poson's Ratio of PCC or old AC                  
*SET,U3,0.20       !the Poson's Ratio of Base                 
*SET,U4,0.40       !the Poson's Ratio of subgrade
mp,ex,1,M1   
mp,nuxy,1,U1  
mp,ex,2,M2  
mp,nuxy,2,U2  
mp,ex,3,M3
mp,nuxy,3,U3  
mp,ex,4,M4   
mp,nuxy,4,U4   
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
ET,4,plane42
mat,1
CSYS,1                       ! CYLINDRICAL COORDINATE SYSTEM
N,1
NGEN,5,20,1
N,2,0.25
N,82,0.25,180
FILL,2,82,3,22,20
CSYS,0  
N,3,1
N,23,1,1
N,43,0,1
N,63,-1,1
N,83,-1,0            !the node of closest to crack tip element
NGEN,a+1,100,1,83,,,,(LOAD_D/2-(LOAD_S/2+LOAD_W/2))/a,0.05
NGEN,b+1,100,100*a+1,100*a+83,,,,LOAD_W/b
NGEN,c+1,100,100*(a+b)+1,100*(a+b)+83,,,,(LOAD_S-LOAD_W)/c
NGEN,b+1,100,100*(a+b+c)+1,100*(a+b+c)+83,,,,LOAD_W/b
NGEN,d+1,100,100*(a+2*b+c)+1,100*(a+2*b+c)+83,,,,(road_W-LOAD_D/2-LOAD_S/2-LOAD_W/2)/d,N_ratio   
*do,i,0,a+2*b+c+d-1
E,2+i*100,22+i*100,1+i*100,1+i*100,102+i*100,122+i*100,101+i*100,101+i*100
EGEN,4,20,-1
*enddo
*do,i,0,a+2*b+c+d-1
E,2+i*100,3+i*100,23+i*100,22+i*100,102+i*100,103+i*100,123+i*100,122+i*100
EGEN,4,20,-1
*enddo
TYPE,2                          
EMODIF,1
*REPEAT,4*(a+2*b+c+d),1
NUMMRG,NODE                 ! MERGE COINCIDENT NODES
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
/NERR,DEFA  
/OUTPUT
k,1,1,0
k,2,AC_T-crack_L,0  
k,3,AC_T-crack_L,1
k,4,1,1
k,5,AC_T-crack_L,LOAD_L/2
k,6,1,LOAD_L/2
k,7,0,LOAD_L/2
k,8,0,1
k,9,-1,LOAD_L/2
k,10,-1,1
k,11,-1,0
K,12,-crack_L,LOAD_L/2
K,13,-crack_L,crack_W
K,14,-crack_L,0        !Attention! this still is the AC layer (wearing course)
a,1,2,3,4
a,4,3,5,6
a,8,4,6,7
a,10,8,7,9      
a,13,10,9,12
a,14,11,10,13     !#1~6 areas is for the AC layer
lsel,s,,,1
lsel,a,,,3
lsel,a,,,6
lesize,all,,,N_AC_T   !the amount of element from crack element to road surface
lsel,s,,,14
lsel,a,,,17
lsel,a,,,15
lesize,all,,,N_AC_B   !the amount of element from crack element to road bottom  
lsel,s,,,5
lsel,a,,,7
lsel,a,,,10
lsel,a,,,13
lsel,a,,,16
lesize,all,,,N_Load   !the amount of element for part of load area (the other part have been for crack element)  
lsel,s,,,2
lsel,a,,,4
lsel,a,,,8
lsel,a,,,9
lsel,a,,,11
lsel,a,,,12
lsel,a,,,18
lsel,a,,,19
lesize,all,,,1   
type,4
asel,all
mshkey,1
amesh,all     !mesh area #1~6
k,15,-PCC_T-crack_L               !for PCC or old AC layer (along X axis,namely the road depth direction)
k,16,-BL_T-PCC_T-crack_L          !for base layer (along X axis,namely the road depth direction)
K,17,-SUB_T-BL_T-PCC_T-crack_L    !for subgrade (along X axis,namely the road depth direction)
K,18,0,road_L                     !along Y direction,namely the road length direction (driving direction)
l,14,15     !#20 line. for draging the PCC or old AC layer
l,15,16     !#21 line. for draging the Base layer thickness
l,16,17     !#22 line. for draging the Subgrade calculated thickness
            
l,7,18      !#23 line. for draging the residual length of calculated road
lesize,20,,,N_PCC       !PCC or old AC layer
lesize,21,,,N_BL        !Base layer
lesize,22,,,N_SUB       !subgrade
lesize,23,,,15,15        !it is residual length of calculated road length,and it is devided to 15 segments
lsel,s,loc,x,-crack_L
lgen,2,all              !creat the line #24,#25, for the separate node or the bonding condition between AC and PCC layers
lsel,a,,,20
lsel,a,,,21
lsel,a,,,22
adrag,24,25,,,,,20,21,22   !form the area for draging PCC or old AC layer, base layer and subgrade
asel,all
*if,if_c_wid,EQ,1,then
  adele,8                   !delete the area for the crack width of PCC/old-AC
*endif
lsel,all
lsel,s,loc,y,LOAD_L/2
lsel,r,loc,x,-crack_L,AC_T-crack_L
lsel,a,,,23
adrag,all,,,,,,23         !form the area for draging the AC layer (wearing course)
lsel,all
lsel,s,loc,y,LOAD_L/2
lsel,r,loc,x,-SUB_T-BL_T-PCC_T-crack_L,-crack_L
lsel,a,,,23
adrag,all,,,,,,23         !form the area for draging the PCC or old AC layer, base layer or subgrade
                          !Attention!! It is separate between the AC layer and PCC or old AC layer
asel,all
asel,u,,,1,6
mshkey,1
amesh,all
!following lines are for leading to form the volum
!Attention,the volume are draged along Z axis, namely, the road width direction
numstr,kp,100
numstr,line,100
k,100,0,0,0
k,101,0,0,LOAD_D/2-(LOAD_S/2+LOAD_W/2)
k,102,0,0,LOAD_D/2-(LOAD_S/2-LOAD_W/2)
K,103,0,0,LOAD_D/2+(LOAD_S-LOAD_W)/2
K,104,0,0,LOAD_D/2+LOAD_S/2+LOAD_W/2
K,105,0,0,road_W
l,100,101   !#100 line, the distance between central to edge of tire
l,101,102   !#101 line, the width of one tire contact area
l,102,103   !#102 line, the distance between dual wheel
l,103,104   !#103 line, the width of the other tire contact area
l,104,105   !#104 line, the distance between edge of tire to road edge
lesize,100,,,a,0.05         !folowing 5 lines are the line that will induce the volume
lesize,101,,,b
lesize,102,,,c
lesize,103,,,b
lesize,104,,,d,N_ratio     
asel,all
lsel,all
type,1
mat,1
asel,s,loc,x,-crack_L,AC_T-crack_L
asel,u,loc,x,-crack_L
asel,u,loc,x,AC_T-crack_L
vdrag,all,,,,,,100,101,102,103,104                       !from the AC layer
mat,2
asel,s,loc,x,-PCC_T-crack_L,-crack_L
asel,u,loc,x,-PCC_T-crack_L
asel,u,loc,x,-crack_L
vdrag,all,,,,,,100,101,102,103,104                       !from the PCC of old AC layer
mat,3
asel,s,loc,x,-BL_T-PCC_T-crack_L,-PCC_T-crack_L
asel,u,loc,x,-BL_T-PCC_T-crack_L
asel,u,loc,x,-PCC_T-crack_L
vdrag,all,,,,,,100,101,102,103,104                       !from the Base layer
mat,4
asel,s,loc,x,-SUB_T-BL_T-PCC_T-crack_L,-BL_T-PCC_T-crack_L
asel,u,loc,x,-SUB_T-BL_T-PCC_T-crack_L
asel,u,loc,x,-BL_T-PCC_T-crack_L
vdrag,all,,,,,,100,101,102,103,104                       !from the subgrade

nsel,all
esel,all
asel,all
aclear,all             !delete the area meshed
lsel,s,,,20,23
lsel,a,,,100,101
ldel,all               !delete the line
nsel,all
esel,all
lsel,all
*if,if_bond,EQ,1,then          !when the interface is bonding fully
  nsel,all
  nummrg,node
*else                          !when the interface is separate
  nsel,all
  nsel,u,loc,x,-crack_L
  nummrg,node
*endif
*if,if_has_c,EQ,0,then                    !!!if the PCC or old AC has not crack
  nsel,s,loc,y,0
  nsel,u,loc,x,-crack_L+0.001,-0.001
  dsym,symm,y   
      
  vsel,s,loc,x,-crack_L,-SUB_T-BL_T-PCC_T-crack_L
  nsel,s,loc,y,0
  nsel,r,loc,x,-crack_L
  nslv,r,1                             ! chose the node only belong to the PCC/old-AC layer
  dsym,symm,y                          ! symmetric B.C.'S at y=0 only the node of belonging to the PCC/old-AC layer
*else                                  !!!if the PCC or old AC has crack, has or hasnot width
  nsel,s,loc,y,0
  nsel,u,loc,x,-crack_L-PCC_T+0.001,-0.001  !!the value 0.001, is for excluding the node
  dsym,symm,y  
  
*endif
nsel,all
vsel,all
asel,all
finish
/solu
nsel,all
nsel,s,loc,x,-SUB_T-BL_T-PCC_T-crack_L
D,all,uy,0                                  !the bottom of subgrad is constrained fully
D,all,ux,0
nsel,all
nsel,s,loc,y,road_L
D,all,uy,0                                   ! one side of pavement is constrained partly
nsel,s,loc,z,0
D,all,all,0            !symmetric B.C.'S at z=0; namely,the half width of road will be calculated
nsel,all
nsel,s,loc,z,road_W
D,all,ux,0  
D,all,uy,0                                   ! one side of pavement is constrained partly
nsel,all
asel,all
asel,s,loc,x,AC_T-crack_L
asel,r,loc,y,0,LOAD_L/2
asel,r,loc,z,LOAD_D/2-(LOAD_S/2+LOAD_W/2),LOAD_D/2+LOAD_S/2+LOAD_W/2
asel,u,loc,z,LOAD_D/2-(LOAD_S/2-LOAD_W/2),LOAD_D/2+(LOAD_S-LOAD_W)/2
sfa,all,1,pres,P
dtran
/DELETE,FRACT,MAC
SOLVE
nsel,s,loc,y,0         
nsel,r,loc,z,LOAD_D/2-(LOAD_S/2-LOAD_w/2)
nsel,r,loc,x,-0.25,0
/post1  
PATH,KI,3,,20
PPATH,1,1501              
PPATH,2,4822
PPATH,3,1582
KCALC,,,1