SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) C
INCLUDE 'ABA_PARAM.INC' C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),
2 DDSDDT(NTENS),DRPLDE(NTENS),
3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) C
DIMENSION STRANT(6),CFULL(6,6),CDFULL(6,6),DCDDFV(6,6), 1 DCDDMV(6,6),DCDDDV(6,6),DFFDE(6),DFMDE(6),DFDDE(6), 2 TEMP1(6),TEMP2(6),TEMP3(6),TEMP4(6),TEMP5(6),TEMP6(6), 3 TEMP7(6,6)
PARAMETER (ZERO = 0.D0,ONE = 1.D0,TWO = 2.D0,THREE = 3.D0) C 输入材料弹性参数
e1=props(1) e2=props(2) e3=props(3) u12=props(4) u13=props(5) u23=props(6) g12=props(7) g13=props(8) g23=props(9)
C 输入材料的强度参数，用于失效判断。 xt=props(10) xc=props(11) yt=props(12) yc=props(13) s12=props(14) s13=props(15) s23=props(16) zt=props(17) zc=props(18)
eta=props(19) 粘性系数
gc1=props(20) 复合材料的1方向的临界能量释放率 gc2=props(21) gc3=props(22)
!CFULL 无损伤下的刚度阵 u21=(u12/e1)*e2 u31=(u13/e1)*e3 u32=(u23/e2)*e3
fact=1-u12*u21-u23*u32-u13*u31- & 2*u21*u32*u13
factor=fact/(e1*e2*e3) CFULL(1,1)=(1-u23*u32)/(e2*e3*factor) CFULL(2,2)=(1-u13*u31)/(e1*e3*factor) CFULL(3,3)=(1-u12*u21)/(e1*e2*factor)
CFULL(1,2)=(u12+u32*u13)/(e1*e3*factor) CFULL(2,3)=(u32+u12*u31)/(e1*e3*factor) CFULL(1,3)=(u31+u21*u32)/(e3*e2*factor) CFULL(3,1)=CFULL(1,3) CFULL(3,2)=CFULL(2,3) CFULL(2,1)=CFULL(1,2) CFULL(4,4)=g12 CFULL(5,5)=g13 CFULL(6,6)=g23 C
C RECOVER ELASTIC STRAINS 更新应变 C
DO 10 K1=1,NTENS
STRANT(K1)=STATEV(K1)+DSTRAN(K1) 10 CONTINUE C
C STORE STRAINS IN STATE VARIABLE ARRAY 存储应变值，用于下一步的应变更新 C
DO 20 K1=1,NTENS
STATEV(K1)=STRANT(K1) 20 CONTINUE
!damage strain 计算各种失效情况下发生失效时的应变极限值 epsilonxt=xt/CFULL(1,1) epsilonxc=xc/CFULL(1,1) epsilonyt=yt/CFULL(2,2) epsilonyc=yc/CFULL(2,2) epsilonzt=zt/CFULL(3,3) epsilonzc=zc/CFULL(3,3)
epsilons12=s12/CFULL(4,4) epsilons13=s13/CFULL(5,5) epsilons23=s23/CFULL(6,6)
!fiber tense damage 失效准则
if(STRANT(1).GE.0.0)then
eft=(STRANT(1)/epsilonxt)**2+(STRANT(4)/epsilons12)**2+
& (STRANT(5)/epsilons13)**2 efc=0.0 else eft=0.0
!fiber compress buckling
efc=(STRANT(1)/xc)**2 end if !matrix crack
epm=STRANT(2)+STRANT(3) if(epm.GE.0.0)then
emt=(STRANT(2)+STRANT(3))**2/(epsilonyt*epsilonzt)- &STRANT(2)*STRANT(3)/epsilons23**2+ &(STRANT(4)/epsilons12)**2+ &(STRANT(5)/epsilons13)**2+ &(STRANT(6)/epsilons23)**2 emc=0.0 else emt=0.0 !matrix cruck
emc=(STRANT(2)+STRANT(3))**2/(epsilonyc*epsilonzc)+ &((STRANT(2)+STRANT(3))/epsilonyc)* &(epsilonyc/2.0/epsilons12-1.0)-
&STRANT(2)*STRANT(3)/epsilons23**2+ &(STRANT(4)/epsilons12)**2+ &(STRANT(5)/epsilons13)**2+ &(STRANT(6)/epsilons23)**2 end if !delemination
if(STRANT(3).GE.0.0)then
edt=(STRANT(3)/epsilonzt)**2+(STRANT(6)/epsilons23)**2+ & (STRANT(5)/epsilons13)**2 edc=0.0 else edt=0.0
edc=(STRANT(3)/epsilonzc)**2+(STRANT(6)/epsilons23)**2+ & (STRANT(5)/epsilons13)**2 end if
! Open(16,File='D:/1.txt',status='unknown') ! write(16,*)ef,em,ed ef=eft+efc em=emt+emc ed=edt+edc
!damage or not 判断是否发生失效 !fiber damage
if(ef.GT.1.0)then 纤维发生损伤 ff=SQRT(ef) 计算ff
if(STRANT(1).GE.0.0)then 若为纤维拉断
df=1.0-exp((-1.0)*xt*epsilonxt*CELENT*(ff-1.0)/gc1)/ff 计算纤维的损伤变量df
dfdff=xt*epsilonxt*CELENT/gc1*
计算，
更新雅克比矩阵用
& exp((-1.0)*xt*epsilonxt*CELENT*(ff-1.0)/gc1)/ff+ & exp((-1.0)*xt*epsilonxt*CELENT*(ff-1.0)/gc1)/ff**2 DO i=1,6
DFFDE(i)=0.0
计算
，更新雅克比矩阵用
END DO
DFFDE(1)=1.0/ff*STRANT(1)/epsilonxt**2 DFFDE(4)=1.0/ff*STRANT(4)/epsilons12**2 DFFDE(5)=1.0/ff*STRANT(5)/epsilons13**2
if( df.GT.STATEV(7))then 因为损伤不可逆，保存df的最大值，用到粘性规律中 STATEV(7)=df else
df=STATEV(7) endif
else 若为纤维压缩屈曲
df=1.0-exp((-1.0)*xc*epsilonxc*CELENT*(ff-1.0)/gc1)/ff 计算纤维的损伤变量df
dfdff=xc*epsilonxc*CELENT/gc1* & exp((-1.0)*xc*epsilonxc*CELENT*(ff-1.0)/gc1)/ff+ & exp((-1.0)*xc*epsilonxc*CELENT*(ff-1.0)/gc1)/ff**2
计算计算
DO i=1,6
DFFDE(i)=0.0 END DO
DFFDE(1)=1.0/ff*STRANT(1)/epsilonxc**2
if( df.GT.STATEV(7))then 因为损伤不可逆，保存df的最大值，用到粘性规律中
STATEV(7)=df else
df=STATEV(7) endif endif
else 纤维没有发生损伤，纤维的损伤变量等都置零 df=0.0 dfdff=0.0 DO i=1,6 DFFDE(i)=0.0 END DO endif
!matrix damage 基体损伤，思路与纤维损伤相同 if(em.GT.1.0)then fm=SQRT(em)
if(epm.GE.0.0)then
dm=1.0-exp((-1.0)*yt*epsilonyt*CELENT*(fm-1.0)/gc2)/fm dmdfm=yt*epsilonyt*CELENT/gc2*
& exp((-1.0)*yt*epsilonyt*CELENT*(fm-1.0)/gc2)/fm+ & exp((-1.0)*yt*epsilonyt*CELENT*(fm-1.0)/gc2)/fm**2 DO i=1,6 DFMDE(i)=0.0 END DO
DFMDE(2)=1.0/2.0/fm*
&(2.0*(STRANT(2)+STRANT(3))/epsilonyt/epsilonzt- &STRANT(3)/epsilons23**2) DFMDE(3)=1.0/2.0/fm*
& (2.0*(STRANT(2)+STRANT(3))/epsilonyt/epsilonzt- & STRANT(2)/epsilons23**2)
DFMDE(4)=1.0/fm*STRANT(4)/epsilons12**2 DFMDE(5)=1.0/fm*STRANT(5)/epsilons13**2 DFMDE(6)=1.0/fm*STRANT(6)/epsilons23**2 if(dm.GT.STATEV(8))then STATEV(8)=dm else
dm=STATEV(8) endif else
dm=1.0-exp(-1.0*yc*epsilonyc*CELENT*(fm-1.0)/gc2)/fm dmdfm=yc*epsilonyc*CELENT/gc2*
& exp((-1.0)*yc*epsilonyc*CELENT*(fm-1.0)/gc2)/fm+ & exp((-1.0)*yc*epsilonyc*CELENT*(fm-1.0)/gc2)/fm**2 DO i=1,6
DFMDE(i)=0.0 END DO
DFMDE(2)=1.0/2.0/fm*