SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,NDI,NSHR,NTENS,
3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,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),
4 DFGRD0(3,3),DFGRD1(3,3)
C
C LOCAL ARRAYS
C ---------------------------------------------------------
C EELAS - ELASTIC STRAINS
C EPLAS - PLASTIC STRAINS
c AKPHA - SHIFT TENSOR
C FLOW - DIRECTION OF PLASTIC FLOW
C OLDS - STRESS AT START OF INCREMENT
C OLDPL - PLASTIC STRAINS AT START OF INCREMENT
C ---------------------------------------------------------
C
DIMENSION EELAS(6),EPLAS(6),ALPHA(6),FLOW(6),OLDS(6),OLDPL(6),
1 HARD(3)
C
DOUBLE PRECISION I1
C
PARAMETER (ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,SIX=6.D0,
1 ENUMAX=.4999D0,NEWTON=100,TOLER=1.0D-6)
C
C -----------------------------------------------------------
C UMAT FOR ISOTROPIC ELASTICITY AND DRUCKER-PRAGER PLASTICITY
C WITH MULTILINEAR ISOTROPIC AND KINEMATIC HARDENING
C CANNOT BE USED FOR PLANE STRESS
C -----------------------------------------------------------
C PROPS(1) - E
C PROPS(2) - NU
C PROPS(3) - A, DRUCKER-PRAGER CONSTANT, LET A=ZERO FOR
C CLASSICAL VON MISES PLASTICITY
C PROPS(4) - BETA, SCALAR PARAMETER TO DETERMINE
C PERCENTAGE OF EACH TYPE OF HARDENING
C 0.0 <= BETA <= 1.0
C PROPS(5..) - SYIELD VS PLASTIC STRAIN DATA
C -----------------------------------------------------------
C
C WRITE(7,100)
C 100 FORMAT(//,30X,'*** TEST MESSAGE')
IF (NDI.NE.3) THEN
WRITE(6,1)
1 FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',
$ 'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')
ENDIF竹报平安
C
C ELASTIC PROPERTIES
C
EMOD=PROPS(1)
ENU=MIN(PROPS(2),ENUMAX)
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG中国历史皇帝
ELAM=(EBULK3-EG2)/THREE
EBULK=EBULK3/THREE
A=PROPS(3)
BETA=PROPS(4)
C READ INITIAL YIELD STRESS
SYINIT=PROPS(5)
C
C ELASTIC STIFFNESS
C
DO K1=1,NDI
DO K2=1,NDI
DDSDDE(K2,K1)=ELAM
END DO
DDSDDE(K1,K1)=EG2+ELAM
END DO
DO K1=NDI+1,NTENS
DDSDDE(K1,K1)=EG
END DO
C RECOVER ELASTIC AND PLASTIC STRAINS AND AND SHIFT TENSOR AND
C ROTATE NOTE: USE CODE 1 FOR (TENSOR) STRESS, CODE 2 FOR
C (ENGINEERING) STRAIN
泛泛之交
C
CALL ROTSIG(STATEV( 1),DROT,EELAS,2,NDI,NSHR)
CALL ROTSIG(STATEV( NTENS+1),DROT,EPLAS,2,NDI,NSHR)
CALL ROTSIG(STATEV(2*
NTENS+1),DROT,ALPHA,1,NDI,NSHR)
EQPLAS=STATEV(3*NTENS+1)
C
C SAVE STRESS AND PLASTIC STRAINS AND
C CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN
C
DO K1=1,NTENS
OLDS(K1)=STRESS(K1)
OLDPL(K1)=EPLAS(K1)
EELAS(K1)=EELAS(K1)+DSTRAN(K1)
DO K2=1,NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
END DO
END DO
C
C CALCULATE EQUIVALENT VON MISES STRESS
C
SMISES=(STRESS(1)-ALPHA(1)-STRESS(2)+ALPHA(2))**2
1 +(STRESS(2)-ALPHA(2)-STRESS(3)+ALPHA(3))**2
1 +(STRESS(3)-ALPHA(3)-STRESS(1)+ALPHA(1))**2
DO K1=NDI+1,NTENS
SMISES=SMISES+SIX*(STRESS(K1)-ALPHA(K1))**2
END DO
SMISES=SQRT(SMISES/TWO)
SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
I1=SHYDRO*THREE
C CALCULATE EQUIVALENT DRUCKER-PRAGER STRESS
C EQUATION 3
DPYIEL=SMISES+A*I1
C
C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
C
NVALUE=(NPROPS-4)/2
CALL UHARD(SYIELD,HARD,EQPLAS,EQPLASRT,TIME,DTIME,TEMP,
1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NSTATV,
2 STATEV,NUMFIELDV,PREDEF,DPRED,NVALUE,PROPS)
C
C CALCULATE PREVIOUS YIELD VALUE AND DETERMINE IF
C ACTIVELY YIELDING
C
C STORE PREVIOUS TABLE VALUE FOR SYIELD
SYT1=SYIELD
SYPREV=SYINIT+(SYT1-SYINIT)*BETA
IF (DPYIEL.GT.(ONE+TOLER)*SYPREV) THEN
C
C ACTIVELY YIELDING
C SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS
C CALCULATE THE FLOW DIRECTION
C EQUATIONS 13 AND 30
C
DO K1=1,NDI
FLOW(K1)=(STRESS(K1)-ALPHA(K1)-SHYDRO)/SMISES
END DO
DO K1=NDI+1,NTENS
FLOW(K1)=(STRESS(K1)-ALPHA(K1))/SMISES
END DO
C
C SOLVE FOR NEW YIELD STRESS
C AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATION
C
SYIELD=SYPREV
DEQPL=ZERO
DO KEWTON=1,NEWTON
RHS=SMISES-SYPREV-(THREE*EBULK3*A*A+EG3)*DEQPL
$ -SYIELD+A*I1+SYT1
C EQUATION 38
DEQPL=DEQPL+RHS/(THREE*EBULK3*A*A+EG3+HARD(1))
EQPLA2=EQPLAS+DEQPL
CALL UHARD(SYIELD,HARD,EQPLA2,EQPLASRT,TIME,DTIME
$ ,TEMP,DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME
$ ,NSTATV,STATEV,NUMFIELDV,PREDEF,DPRED,NVALUE,PROPS)
C
IF(ABS(RHS).LT.TOLER*SYPREV) GOTO 10
END DO
C
C WRITE WARNING MESSAGE TO THE .MSG FILE
C
杨绛先生语录WRITE(7,2) NEWTON
2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
1 'CONVERGE AFTER ',I3,' ITERATIONS')
10 CONTINUE
C
C STORE THE CURRENT TABLE YIELD VALUE
SYT2=SYIELD
C UPDATE THE EQUIVALENT PLASTIC STRAIN
EQPLAS=EQPLAS+DEQPL
C
C SOLVE FOR UPDATED YIELD STRESS
C EQUATION 41
SYIELD=SYINIT+(SYT2-SYINIT)*BETA
C
C CALCULATE NEEDED CONSTANTS
C EQUATIONS 21 - 24
P=(-EBULK*A*SYIELD-SHYDRO*EG+EBULK*A*SMIS
ES)/(EBULK3*A*A+EG)
Q=(EBULK3*A*A*SMISES+SYIELD*EG-A*I1*EG)/(EBULK3*A*A+EG)
DEP=(A*(SMISES-SYIELD+A*I1))/(EBULK3*A*A+EG)
DEQ=(SMISES-SYIELD+A*I1)/(THREE*EBULK3*A*A+EG3)
C
C UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND SHIFT TENSOR
C
DO K1=1,NDI
C EQUATION 37
ALPHA(K1)=ALPHA(K1)+(SYT2-SYT1)*FLOW(K1)*(ONE-BETA)
C CHANGE IN PLASTIC STRAIN FROM EQUATION 26
EPLAS(K1)=EPLAS(K1)+DEP/THREE+THREE*DEQ*FLOW(K1)/TWO
EELAS(K1)=EELAS(K1)-(DEP/THREE+THREE*DEQ*FLOW(K1)/TWO)
C EQUATION 27
STRESS(K1)=ALPHA(K1)-P+Q*FLOW(K1)
END DO
DO K1=NDI+1,NTENS
C EQUATION 37
ALPHA(K1)=ALPHA(K1)+(SYT2-SYT1)*FLOW(K1)*(ONE-BETA)
C CHANGE IN PLASTIC STRAIN FROM EQUATION 26
EPLAS(K1)=EPLAS(K1)+THREE*DEQ*FLOW(K1)
EELAS(K1)=EELAS(K1)-THREE*DEQ*FLOW(K1)
C EQUATION 27
STRESS(K1)=ALPHA(K1)+Q*FLOW(K1)基冈
END DO
C
C CALCULATE THE PLASTIC DISSIPATION
C
SPD=ZERO
DO K1=1,NTENS
SPD=SPD+(STRESS(K1)+OLDS(K1))*(EPLAS(K1)-OLDPL(K1))/TWO
END DO
C
C FORMULATE THE JACOBIAN (MATERIAL TANGENT)
C FIRST CALCULATE THE EFFECTIVE MODULI
C
EFFG=EG*SYIELD/SMISES
吃槟榔
EFFG2=TWO*EFFG
EFFG3=THREE*EFFG
EFFLAM=(EBULK3-EFFG2)/THREE
EFFHRD=EG3*HARD(1)/(EG3+HARD(1))-EFFG3
DO K1=1,NDI
DO K2=1,NDI
DDSDDE(K2,K1)=EFFLAM
END DO
DDSDDE(K1,K1)=EFFG2+EFFLAM
END DO
DO K1=NDI+1,NTENS
DDSDDE(K1,K1)=EFFG
END DO
DO K1=1,NTENS
DO K2=1,NTENS
DDSDDE(K2,K1)=DDSDDE(K2,K1)+EFFHRD*FLOW(K2)*FLOW(K1)
END DO
END DO
ENDIF
C
C STORE ELASTIC STRAINS, PLASTIC STRAINS AND SHIFT TENSOR
C IN STATE VARIABLE ARRAY
C
DO K1=1,NTENS
STATEV(K1)=EELAS(K1)
STATEV(K1+NTENS)=EPLAS(K1)
STATEV(K1+2*NTENS)=ALPHA(K1)
END DO
STATEV(3*NTENS+1)=EQPLAS
C
RETURN
END
C
C
SUBROUTINE UHARD(SYIELD,HARD,EQPLAS,EQPLASRT,TIME,DTIME,TEMP,
$ DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NSTATV,STATEV
$ ,NUMFIELDV,PREDEF,DPRED,NVALUE,PROPS)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION HARD(3), STATEV(NSTATV),TIME(*),
1 PREDEF(NUMFIELDV),DPRED(*),PROPS(*)
C
DIMENSION TABLE(2,NVALUE)
C
电脑声音怎么设置
PARAMETER(ZERO=0.D0)
C
C FILL IN SY VS EQPLAS TABLE FROM PROPS ARRAY
I=5
DO K2=1,NVALUE
如何增强宝宝抵抗力
TABLE(1,K2)=PROPS(I)
TABLE(2,K2)=PROPS(I+1)
I=I+2
END DO
C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
C
SYIELD=TABLE(1,NVALUE)
HARD(1)=ZERO
C
C IF MORE THAN ONE ENTRY, SEARCH TABLE
C
IF(NVALUE.GT.1) THEN
DO K1=1,NVALUE-1
EQPL1=TABLE(2,K1+1)
IF(EQPLAS.LT.EQPL1) THEN
EQPL0=TABLE(2,K1)
IF(EQPL1.LE.EQPL0) THEN
WRITE(7,1)
1 FORMAT(//,30X,'***ERROR - PLASTIC STRAIN MUST BE ',
$ 'ENTERED IN ASCENDING ORDER')
CALL XIT
ENDIF
C
C CURRENT YIELD STRESS AND HARDENING
C
DEQPLS=EQPL1-EQPL0
SYIEL0=TABLE(1,K1)
SYIEL1=TABLE(1,K1+1)
DSYIEL=SYIEL1-SYIEL0
HARD(1)=DSYIEL/DEQPLS
SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD(1)
GOTO 10
ENDIF
END DO
10 CONTINUE
ENDIF
RETURN
END
