各向同性材料损伤本构模型
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, + RPL,DDSDDT,DRPLDE,DRPLDT,
+ STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, + NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, + CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
INCLUDE 'ABA_PARAM.INC' CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV), + DDSDDE(NTENS,NTENS),DDSDDT(NTENS),
+ DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS), + TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS), + COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
DIMENSION STRANT(6),TSTRANT(4),PT(1) DIMENSION OLD_STRESS(6)
DIMENSION DOLD_STRESS(6),D_STRESS(6)
DIMENSION C(6,6),CD(6,6),DSTRESS(6),BSTRESS(6),ROOT(3), + DFMNDE(6),DDMDE(6),DCDDM(6,6),ATEMP1(6), ATEMP2(6)
PARAMETER (ZERO=0.D0,ONE=1.D0,TWO=2.D0,FOUR=4.D0,HALF = 0.5D0) C start
C IF (NPROPS.LT.2) THEN
C WRITE(7,*) '** ERROR: UMAT REQUIRES *NPROPS=2' C STOP C END IF
E11 =PROPS(1) V12 =PROPS(2)
G12 =PROPS(1)/TWO/(ONE+PROPS(2))
C Critical values of stresses XT=PROPS(3) XC=PROPS(4) XS=PROPS(5)
GX=PROPS(6) !Fracture energy in matrix ETA=0.001 C Current strain
DO I = 1, NTENS
STRANT(I) = STRAN(I) + DSTRAN(I) END DO C Stiffness
DO I = 1, 6 DO J = 1, 6 C(I,J)=ZERO
END DO END DO
ATEMP = (1+V12)*(1-TWO*V12) C(1,1) = E11*(1-V12)/ATEMP C(2,2) = E11*(1-V12)/ATEMP C(3,3) = E11*(1-V12)/ATEMP C(1,2) = E11*V12/ATEMP C(1,3) = E11*V12/ATEMP C(2,3) = E11*V12/ATEMP C(4,4) = G12 C(5,5) = G12 C(6,6) = G12
DO I = 2, 6
DO J = 1, I-1 C(I,J) = C(J,I)
END DO END DO
C Critical values of strains
XET=XT/(C(1,1)-2*V12*C(1,2)) XEC=XC/(C(1,1)-2*V12*C(1,2)) XES=XS/C(4,4)
DMOLD = STATEV(1)
C Strain initiation criterion
A11 = STRANT(1)**TWO+STRANT(2)**TWO+STRANT(3)**TWO A12 = A11 / XET / XEC
A21 = STRANT(1)+STRANT(2)+STRANT(3) A22 = (XEC - XET) / XEC / XET * A21
A31 = STRANT(4)**TWO+STRANT(5)**TWO+STRANT(6)**TWO A32 = A31 / XES**TWO
A1= A12 + A22 + A32
C B11 = STRANT(2)**TWO C B12 = B11 / XET / XEC C B21 = STRANT(2)
C B22 = (XEC - XET) / XEC / XET * B21 C B31 = STRANT(5)**TWO C B32 = B31 / XES**TWO
C B1= B12 + B22 + B32
C C11 = STRANT(3)**TWO C C12 = C11 / XET / XEC C C21 = STRANT(3)
C C22 = (XEC - XET) / XEC / XET * C21 C C31 = STRANT(6)**TWO C C32 = C31 / XES**TWO
C C1= C12 + C22 + C32
STATEV(2)=A1
C STATEV(3)=B1 C STATEV(4)=C1
FMN = ZERO
IF (A1.GT.ZERO) THEN FMN =SQRT(A1)
C IF (B1.GT.ONE) THEN C FMN =FMN+SQRT(B1)
C IF(C1.GT.ONE) THEN C FMN =FMN+SQRT(C1) C END IF C END IF END IF STATEV(5)=FMN
C write(*,*) FMN DM = ZERO
DDMDFMN = ZERO DO I = 1, 6 DFMNDE(I) = ZERO DDMDE(I) = ZERO END DO
IF (FMN .GT. ONE) THEN C CALCULATE DM, DDMDFMN C WRITE(6,*)FMN
T1 = (C(1,1)-2*V12*C(1,2)) * XET**2 * CELENT / GX T2 = (ONE - FMN) * T1 DM = ONE - EXP(T2)/FMN
C WRITE(6,*)'T1 ',T1,' T2', T2, ' DM', DM C write(*,*) DM
C CALCULATE THE DERIVATIVE OF DAMAGE VARIABLE WITH RESPECT TO FAILURE C RITERION
DDMDFMN = (ONE / FMN + T1) * (ONE - DM)
C CALCULATE DFMNDE
IF (DM .GT. DMOLD) THEN
DFMNDE(1) = HALF/FMN*(TWO*STRANT(1)+XEC-XET)/XET/XEC DFMNDE(2) = HALF/FMN*(TWO*STRANT(2)+XEC-XET)/XET/XEC DFMNDE(3) = HALF/FMN*(TWO*STRANT(3)+XEC-XET)/XET/XEC
DFMNDE(4) = ONE/FMN*TWO*STRANT(4)/XES**TWO DFMNDE(5) = ONE/FMN*TWO*STRANT(5)/XES**TWO DFMNDE(6) = ONE/FMN*TWO*STRANT(6)/XES**TWO DO I = 1, 6
DDMDE(I) = DFMNDE(I) * DDMDFMN END DO END IF END IF
DM = MAX (DM, DMOLD) C write(6,*) DM
C SAVE THE OLD STRESS TO OLD_STRESS DO I = 1, NTENS
OLD_STRESS(I) = STRESS(I) END DO
C Effective stiffness DO I = 1, 6 DO J = 1, 6 CD(I,J)=C(I,J)
END DO END DO
IF (DM.NE.ZERO) THEN CD(1,1) = (ONE - DM)*C(1,1) CD(1,2) = (ONE - DM)*C(1,2) CD(2,1) = CD(1,2)
CD(2,2) = (ONE - DM)*C(2,2) CD(1,3) = (ONE - DM)*C(1,3) CD(3,1) = CD(1,3)
CD(2,3) = (ONE- DM)*C(2,3) CD(3,2) = CD(2,3)
CD(4,4) = (ONE - DM)*C(4,4) CD(5,5) = (ONE - DM)*C(5,5) CD(6,6) = (ONE - DM)*C(5,5) END IF
C Elastic derivative
DO I = 1, 6 DO J = 1, 6 DCDDM(I,J) = ZERO END DO END DO C
C CALCULATE DC/DDM C
DCDDM(1,1) = -C(1,1) DCDDM(1,2) = -C(1,2) DCDDM(2,1) = -C(2,1) DCDDM(2,2) = -C(2,2) DCDDM(2,3) = -C(2,3) DCDDM(3,2) = -C(3,2)
DCDDM(3,3) = -C(3,3) DCDDM(3,1) = -C(3,1) DCDDM(1,3) = -C(1,3) DCDDM(4,4) = -C(4,4) DCDDM(5,5) = -C(5,5) DCDDM(6,6) = -C(6,6)
C UPDATE THE JACOBIAN DO I = 1, NTENS ATEMP1(I) = ZERO
DO J = 1, NTENS
ATEMP1(I) = ATEMP1(I) + DCDDM(I,J) * STRANT(J) END DO END DO
DDSDDE=0 DO I = 1, NTENS DO J = 1, NTENS
DDSDDE(I,J)=CD(I,J)+(ATEMP1(I)*DDMDE(J))*DTIME/(DTIME+ETA) END DO END DO C Update stresses
DO I = 1, NTENS STRESS(I)=ZERO
DO J = 1, NTENS
C IF(DM.LT.0.5) THEN
STRESS(I)=STRESS(I)+ CD(I,J) * STRANT(J) C ELSE
C STRESS(I)=STRESS(I)+ CD(I,J) * STRANT(J) * (1-DM) C ENDIF END DO END DO
C Energy
DO I = 1, NDI
SSE = SSE + HALF * (STRESS(I) + OLD_STRESS(I)) * DSTRAN(I) END DO
DO I = NDI+1, NTENS
SSE = SSE + (STRESS(I) + OLD_STRESS(I)) * DSTRAN(I) END DO
STATEV(1)=DM
RETURN END