内容发布更新时间 : 2024/5/4 22:08:59星期一 下面是文章的全部内容请认真阅读。

各向同性材料损伤本构模型

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