C
局部变量，用来暂时保存弹性应变、塑性应变分量以及流动方向 DIMENSION EELAS(6),EPLAS(6),FLOW(6) 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) C
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 HARDSUB FOR CURVE OF YIELD STRESS VS. PLASTIC STRAIN
C ---------------------------------------------------------------- C
C ELASTIC PROPERTIES C
获取杨氏模量，泊松比，作为已知量由PROPS向量传入 EMOD=PROPS(1) E ENU=PROPS(2) ν
EBULK3=EMOD/(ONE-TWO*ENU) 3K EG2=EMOD/(ONE+ENU) 2G EG=EG2/TWO G EG3=THREE*EG 3G
ELAM=(EBULK3-EG2)/THREE λ DO K1=1,NTENS DO K2=1,NTENS
DDSDDE(K1,K2)=ZERO END DO END DO
弹性部分，Jacobian矩阵很容易计算
注意，在ABAQUS中，剪切应变采用工程剪切应变的定义?ij?ui,j?uj,i，所以剪切部分模量是G而不是2G! 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
C RECOVER ELASTIC AND PLASTIC STRAINS AND ROTATE FORWARD
C ALSO RECOVER EQUIVALENT PLASTIC STRAIN C
读取弹性应变分量，塑性应变分量，并旋转（调用了ROTSIG），分别保存在EELAS和EPLAS中；
CALL ROTSIG(STATEV( 1),DROT,EELAS,2,NDI,NSHR)
CALL ROTSIG(STATEV(NTENS+1),DROT,EPLAS,2,NDI,NSHR) 读取等效塑性应变 STATEV中的第13个量 EQPLAS=STATEV(1+2*NTENS)
先假设没有发生塑性流动，按完全弹性变形计算试算应力?σ?J.?ε σn?1?σn??σ
C
C CALCULATE PREDICTOR STRESS AND ELASTIC STRAIN C
DO K1=1,NTENS DO K2=1,NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) END DO
EELAS(K1)=EELAS(K1)+DSTRAN(K1) END DO
C计算Mises应力
C CALCULATE EQUIVALENT VON MISES STRESS C
SMISES=(STRESS(1)-STRESS(2))**2+(STRESS(2)-STRESS(3))**2 1 +(STRESS(3)-STRESS(1))**2 DO K1=NDI+1,NTENS
SMISES=SMISES+SIX*STRESS(K1)**2 END DO
SMISES=SQRT(SMISES/TWO)
C 根据当前等效塑性应变，调用HARDSUB得到对应的屈服应力
C GET YIELD STRESS FROM THE SPECIFIED HARDENING CURVE C
NVALUE=NPROPS/2-1
CALL HARDSUB(SYIEL0,HARD,EQPLAS,PROPS(3),NVALUE)
C
C DETERMINE IF ACTIVELY YIELDING
C 如果Mises应力大余屈服应力，屈服发生，计算流动方向 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=(STRESS(1)+STRESS(2)+STRESS(3))/THREE 静水压力 DO K1=1,NDI
FLOW(K1)=(STRESS(K1)-SHYDRO)/SMISES 流动方向 END DO
DO K1=NDI+1,NTENS
FLOW(K1)=STRESS(K1)/SMISES END DO
C根据J2理论并应用Newton-Rampson方法求得等效塑性应变增量 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) CALL
HARDSUB(SYIELD,HARD,EQPLAS+DEQPL,PROPS(3),NVALUE) IF(ABS(RHS).LT.TOLER*SYIEL0) 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更新应力?n?1，应变分量
C UPDATE STRESS, ELASTIC AND PLASTIC STRAINS AND C EQUIVALENT PLASTIC STRAIN C
DO K1=1,NDI
STRESS(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=NDI+1,NTENS
STRESS(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 C
C 计算塑性变形下的Jacobian矩阵
FORMULATE THE JACOBIAN (MATERIAL TANGENT) C FIRST CALCULATE EFFECTIVE MODULI C
EFFG=EG*SYIELD/SMISES EFFG2=TWO*EFFG
EFFG3=THREE/TWO*EFFG2
EFFLAM=(EBULK3-EFFG2)/THREE
EFFHRD=EG3*HARD/(EG3+HARD)-EFFG3
经典弹塑性理论，应力表达式式如下：pr?ij??ij?s??ij?kk133其中：??ijpl??ij??pl,?ij?Sijpr/?pr2对时间取微分：pr?ij?s??ij??ij???s??ij??kk?13?ij?其中：??s???pr?pr?Spr??prSijij(?pr)2??prSij?pr?pr?pr3SijprSklSkl2(?pr)3??ij?1/3?ij??mm)2G(??pr??kl?1/3?kl??mm)3G?ij?kl(??prd?s?p3Gpr?kld??hSkl/?pr?pd?h?3G1pr?ij?s??ij??ij???s??ij??kk那么?3pr?ij?1/3?ij??mm)?kl?1/3?kl??mm)2G(?3G?ij?kl(??ijSkl/?pr?kl?K?ij??mm??s??s?h??pr?prh/3G?1?mm2G?s?mm?s2G?s?ij?3G?ij?kl?s?ij?klG?ij?kl?ij??????[K?ij?mm?]?pr?ij??kl?h?kl?3?pr??prh/3G?1?prG?ij?kl?kl?s2G?s2G?s3G?sh????mm?(K?)?????[?]?????ijmm3?pr?prijh/3G?1?prijklkl?prG?ij?kl?kl?sG?ij?sSkprk?mm??mm?0（Skpr注意：??）k?0prpr2?(?)G?d?s2如果令：G*=prs，?*?K-G*，h??3d?plh?mm?2G*??ij?[?kl?ij??*?ij?则：??3G*]?ij?kl?h/3G?1
=======================================
if (props(7).lt..001) go to 99 c...
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