[Fortran] 纯文本查看 复制代码
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),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
double precision E, v, lambda, mu
INTEGER i,j
double precision qTr, equPlasticStrain, dEquPlasticStrain,
1 yieldn, yieldn1, hydorStaticPressure, H, Sm,
2 Tm, yita, K, Ym
real*8 m
double precision, DIMENSION(NTENS):: flow
C
C user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD
C and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT
E=PROPS(1)
v=PROPS(2)
yield0=PROPS(3)
yield0strain=PROPS(4)
yieldx=PROPS(5)
yieldxstrain=PROPS(6)
R=PROPS(7)
yita0=PROPS(8)
X=PROPS(9)
Kb=PROPS(10)
omiga=PROPS(11)
a=PROPS(12)
v0=PROPS(13)
T=PROPS(14)
C=PROPS(15)
lambda=E*v/((1+v)*(1-2*v))
mu=E/(2*(1+v))
K=E/3*(1-2*v)
do i=1,NDI
do j=1,NDI
DDSDDE(i,j)=lambda
enddo
DDSDDE(i,i)=lambda+2*mu
DDSDDE(i+NDI,i+NDI)=mu
enddo
do i=1,NTENS
do j=1,NTENS
STRESS(i)=STRESS(i)+DDSDDE(i,j)*DSTRAN(j)
enddo
enddo
qTr=0
do i=1,NDI
if (i==NDI) then
j=1
else
j=i+1
End if
qTr=qTr+(STRESS(i)-STRESS(j))**2
enddo
do i=NDI+1, NTENS
qTr=qTr+6*STRESS(i)**2
enddo
qTr=sqrt(qTr/2)
equPlasticStrain=STATEV(1)
dEquPlasticStrain=0
call UMATHard(PROPS,NPROPS, equPlasticStrain+dEquPlasticStrain,
1 yieldn, H)
if (qTr>yieldn) then
c branch2
hydorStaticPressure=0
do i=1,NDI
hydorStaticPressure=hydorStaticPressure+STRESS(i)
enddo
hydorStaticPressure=hydorStaticPressure/3
do i=1,NDI
flow(i)=(STRESS(i)-hydorStaticPressure)/qTr
enddo
do i=NDI+1, NTENS
flow(i)=STRESS(i)/qTr
enddo
c
c yingbianzengliang
c
Sm=(2/3)*((1+v)/(1-v))*mu
Tm=sqrt((1/2)*Sm*Sm)
yita=COSH((Tm*omiga)/(2*Kb))-1
yita=yita*((2*a*Kb*T)/(Sm*v0*yita0))-1
yita=(1/X)*R*EXP(((-1)*a)/yita0)-(C*hydorStaticPressure)/K
Ym=SINH((Tm*omiga)/(2*Kb*T))
Ym=EXP(((-1)*a)/yita0)*Ym
Ym=2*R*Ym
m=1/2
dEquPlasticStrain=sqrt(m)*Ym*flow(i)+(1/3)*yita
call UMATHard(PROPS,NPROPS,
1 equPlasticStrain+dEquPlasticStrain, yieldn1, H)
do i=1, NDI
STRESS(i)=yieldn1*flow(i)+ hydorStaticPressure
enddo
do i=NDI+1, NTENS
STRESS(i)= yieldn1*flow(i)
enddo
else
c branch1
dEquPlasticStrain=0
end if
STATEV(1)=equPlasticStrain+dEquPlasticStrain
RETURN
END
SUBROUTINE UMATHard(PROPS,NPROPS, equPlasticStrain, yield, H)
INCLUDE 'ABA_PARAM.INC'
DIMENSION PROPS(NPROPS)
double precision equPlasticStrain, yield, yield0, H
INTEGER NPROPS
yield0=PROPS(3)
H=(PROPS(5)-PROPS(3))/(PROPS(6)-PROPS(4))
yield=yield0+H*equPlasticStrain
END