求VUMAT 塑性各向同性强化模型
求VUMAT塑性各向同性强化模型感谢了 急用！！！ 小弟对塑性力学有了一些基本认识愿意和大侠们讨论向大侠们学习here is the VUMAT for Isotropic Hardening:C
C Coding for Isotropic Hardening Plasticity VUMAT
C
SUBROUTINE VUMAT(
C Read only -
1 NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV, NPROPS, LANNEAL,
2 STEPTIME, TOTALTIME, DT, CMNAME, COORDMP, CHARLENGTH,
3 PROPS, DENSITY, STRAININC, RELSPININC,
4 TEMPOLD, STRETCHOLD, DEFGRADOLD, FIELDOLD,
5 STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,
6 TEMPNEW, STRETCHNEW, DEFGRADNEW, FIELDNEW,
C Write only -
7 STRESSNEW, STATENEW, ENERINTERNNEW, ENERINELASNEW)
C
INCLUDE 'VABA_PARAM.INC'
C
DIMENSION PROPS(NPROPS), DENSITY(NBLOCK), COORDMP(NBLOCK),
1 CHARLENGTH(NBLOCK), STRAININC(NBLOCK, NDIR+NSHR),
2 RELSPININC(NBLOCK, NSHR), TEMPOLD(NBLOCK),
3 STRETCHOLD(NBLOCK, NDIR+NSHR),DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR),
4 FIELDOLD(NBLOCK, NFIELDV), STRESSOLD(NBLOCK, NDIR+NSHR),
5 STATEOLD(NBLOCK, NSTATEV), ENERINTERNOLD(NBLOCK),
6 ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),
7 STRETCHNEW(NBLOCK, NDIR+NSHR),DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),
8 FIELDNEW(NBLOCK, NFIELDV), STRESSNEW(NBLOCK,NDIR+NSHR),
9 STATENEW(NBLOCK, NSTATEV), ENERINTERNNEW(NBLOCK),
1 ENERINELASNEW(NBLOCK)
CHARACTER*80 CMNAME
C LOCAL ARRAYS
C ----------------------------------------------------------------
C EELAS - ELASTIC STRAINS
C EPLAS - PLASTIC STRAINS
C FLOW - DIRECTION OF PLASTIC FLOW
C ----------------------------------------------------------------
C
DIMENSION EELAS(6),EPLAS(6),FLOW(6), HARD(3)
C
PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, SIX=6.D0,
1 ENUMAX=.4999D0, NEWTON=10, TOLER=1.0D-6,HALF=0.5D0)
C ----------------------------------------------------------------
C UMAT FOR ISOTROPIC ELASTICITY AND ISOTROPIC MISES PLASTICITY
C CANNOT BE USED FOR PLANE STRESS
C ----------------------------------------------------------------
C PROPS(1) - E
C PROPS(2) - NU
C PROPS(3..) - SYIELD AN HARDENING DATA
C CALLS UHARD FOR CURVE OF YIELD STRESS VS. PLASTIC STRAIN
C ----------------------------------------------------------------
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
NVALUE=NPROPS/2-1
C
C CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN
C
IF ( STEPTIME .EQ. ZERO ) THEN
DO K = 1, NBLOCK
TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
STRESSNEW(K,1) = STRESSOLD(K,1)
* + EG2 * STRAININC(K,1) + ELAM * TRACE
STRESSNEW(K,2) = STRESSOLD(K,2)
* + EG2 * STRAININC(K,2) + ELAM * TRACE
STRESSNEW(K,3) = STRESSOLD(K,3)
* + EG2 * STRAININC(K,3) + ELAM * TRACE
STRESSNEW(K,4)=STRESSOLD(K,4) + EG2 * STRAININC(K,4)
IF ( NSHR .GT. 1 ) THEN
STRESSNEW(K,5)=STRESSOLD(K,5) + EG2 * STRAININC(K,5)
STRESSNEW(K,6)=STRESSOLD(K,6) + EG2 * STRAININC(K,6)
END IF
END DO
ELSE
C
C111111111111
C
DOK = 1, NBLOCK
C
TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
STRESSNEW(K,1) = STRESSOLD(K,1)
* + EG2 * STRAININC(K,1) + ELAM * TRACE
STRESSNEW(K,2) = STRESSOLD(K,2)
* + EG2 * STRAININC(K,2) + ELAM * TRACE
STRESSNEW(K,3) = STRESSOLD(K,3)
* + EG2 * STRAININC(K,3) + ELAM * TRACE
STRESSNEW(K,4)=STRESSOLD(K,4) + EG2 * STRAININC(K,4)
IF ( NSHR .GT. 1 ) THEN
STRESSNEW(K,5)=STRESSOLD(K,5) + EG2 * STRAININC(K,5)
STRESSNEW(K,6)=STRESSOLD(K,6) + EG2 * STRAININC(K,6)
END IF
C
SMISES=(STRESSNEW(K,1)-STRESSNEW(K,2))**2 +
1 (STRESSNEW(K,2)-STRESSNEW(K,3))**2 +
2 (STRESSNEW(K,3)-STRESSNEW(K,1))**2
DO 90 K1=NDIR+1,NDIR+NSHR
SMISES=SMISES+SIX*STRESSNEW(K,K1)**2
90 CONTINUE
SMISES=SQRT(SMISES/TWO)
C
C RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE FORWARD
C ALSO RECOVER EQUIVALENT PLASTIC STRAIN
C
C CALL ROTSIG(STATEOLD(K,1), DROT, EELAS, 2, NDIR, NSHR)
C CALL ROTSIG(STATEOLD(K,1+NDIR+NSHR), DROT, EPLAS, 2, NDIR, NSHR)
DO K1=1, NDIR+NSHR
EELAS(K1)=STATEOLD(K,K1)
EPLAS(K1)=STATEOLD(K,K1+NDIR+NSHR)
END DO
EQPLAS=STATEOLD(K,1+2*(NDIR+NSHR))!!!
C
C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE
C
CALL VUHARD(SYIEL0, HARD, EQPLAS, EQPLASRT,TOTALTIME,DT,TEMP,
1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,
2 NUMFIELDV,NVALUE,PROPS(3))
C
C DETERMINE IF ACTIVELY YIELDING
C
IF (SMISES.GT.(ONE+TOLER)*SYIEL0) THEN
C
C ACTIVELY YIELDING
C SEPARATE THE HYDROSTATIC FROM THE DEVIATORIC STRESS
C CALCULATE THE FLOW DIRECTION
C
SHYDRO=(STRESSNEW(K,1)+STRESSNEW(K,2)+STRESSNEW(K,3))/THREE
DO K1=1,NDIR
FLOW(K1)=(STRESSNEW(K,K1)-SHYDRO)/SMISES
END DO
DO K1=NDIR+1,NDIR+NSHR
FLOW(K1)=STRESSNEW(K,K1)/SMISES
END DO
C
C SOLVE FOR EQUIVALENT VON MISES STRESS
C AND EQUIVALENT PLASTIC STRAIN INCREMENT USING NEWTON ITERATION
C
SYIELD=SYIEL0
DEQPL=ZERO
DO KEWTON=1, NEWTON
RHS=SMISES-EG3*DEQPL-SYIELD
DEQPL=DEQPL+RHS/(EG3+HARD(1))
CALL VUHARD(SYIELD,HARD,EQPLAS+DEQPL,EQPLASRT,TOTALTIME,DT,TEMP,
1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,NUMFIELDV,
2 NVALUE,PROPS(3))
IF(ABS(RHS).LT.TOLER*SYIEL0) GOTO 10
END DO
C
C WRITE WARNING MESSAGE TO THE .MSG FILE
C
WRITE(7,2) NEWTON
2FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT',
1 'CONVERGE AFTER ',I3,' ITERATIONS')
10CONTINUE
C
C UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND
C EQUIVALENT PLASTIC STRAIN
C
DO K1=1,NDIR
STRESSNEW(K,K1)=FLOW(K1)*SYIELD+SHYDRO
EPLAS(K1)=EPLAS(K1)+THREE/TWO*FLOW(K1)*DEQPL
EELAS(K1)=EELAS(K1)-THREE/TWO*FLOW(K1)*DEQPL
END DO
DO K1=NDIR+1,NDIR+NSHR
STRESSNEW(K,K1)=FLOW(K1)*SYIELD
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL
EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL
END DO
EQPLAS=EQPLAS+DEQPL
C
C CALCULATE PLASTIC DISSIPATION
C
SPD=DEQPL*(SYIEL0+SYIELD)/TWO
ENDIF
C
C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS
C IN STATE VARIABLE ARRAY
C
DO K1=1, NDIR+NSHR
STATENEW(K,K1)=EELAS(K1)
STATENEW(K,K1+NDIR+NSHR)=EPLAS(K1)
END DO
STATENEW(K,1+2*(NDIR+NSHR)) = STATEOLD(K,1+2*(NDIR+NSHR)) + DEQPL
C
C UPDATE THE SPECIFIC INTERNAL ENERGY -
C
IF ( NSHR .EQ. 1 ) THEN
STRESSPOWER = HALF * (
* ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
* ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
* ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3) ) +
* ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4)
ELSE
STRESSPOWER = HALF * (
* ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
* ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
* ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3) ) +
* ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4) +
* ( STRESSOLD(K,5) + STRESSNEW(K,5) ) * STRAININC(K,5) +
* ( STRESSOLD(K,6) + STRESSNEW(K,6) ) * STRAININC(K,6)
END IF
ENERINTERNNEW(K) = ENERINTERNOLD(K) + STRESSPOWER / DENSITY(K)
C
C UPDATE THE DISSIPATED INELASTIC SPECIFIC ENERGY -
C
PLASTICWORKINC = HALF * ( SYIEL0 + SYIELD ) * DEQPL
ENERINELASNEW(K) = ENERINELASOLD(K) + PLASTICWORKINC / DENSITY(K)
ENDDO
END IF
C
RETURN
END
SUBROUTINE VUHARD(SYIELD,HARD,EQPLAS,EQPLASRT,TIME,DTIME,TEMP,
1 DTEMP,NOEL,NPT,LAYER,KSPT,KSTEP,KINC,CMNAME,
2 NUMFIELDV, NVALUE,TABLE)
INCLUDE 'VABA_PARAM.INC'
CHARACTER*80 CMNAME
DIMENSION HARD(3)
C
DIMENSION TABLE(2, NVALUE)
C
PARAMETER(ZERO=0.D0)
C
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)
1FORMAT(//, 30X, '***ERROR - PLASTIC STRAIN MUST BE ',
1 'ENTERED IN ASCENDING ORDER')
CALL EXIT
ENDIF
C
C CURRENT YIELD STRESS AND HARDENING
C
DEQPL=EQPL1-EQPL0
SYIEL0=TABLE(1, K1)
SYIEL1=TABLE(1, K1+1)
DSYIEL=SYIEL1-SYIEL0
HARD(1)=DSYIEL/DEQPL
SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD(1)
GOTO 10
ENDIF
END DO
10CONTINUE
ENDIF
RETURN
END RathLusia 发表于 2013-11-20 14:25
here is the VUMAT for Isotropic Hardening:
这个你验证过吗？ 收了，走过路过 RathLusia 发表于 2013-11-20 14:25
here is the VUMAT for Isotropic Hardening:
这个你验证过吗 走过路过顺手牵过 感谢楼主 RathLusia 发表于 2013-11-20 14:25
here is the VUMAT for Isotropic Hardening:
同问，验证过吗？ RathLusia 发表于 2013-11-20 14:25
here is the VUMAT for Isotropic Hardening:
已验证，感谢大神无私的教诲 顺手拿走
FSDSFSFSSFSFSFSF 66666666666666666666
页:
[1]