为什么写的子程序只弹性部分起了作用，定义的塑性行为没用？

chenyaorock

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)
     
      INCLUDE'ABA_PARAM.INC'

      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)
     
DIMENSION PS(NDI)
DIMENSION EELAS(NTENS),EPLAS(NTENS),FLOW(NTENS),FS(NTENS)
      PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
      DATA NEWTON,TOLER/100,1.D-6/
      E=PROPS(1)
ENU=PROPS(2)
      KA=PROPS(3)
KB=PROPS(4)

      EBULK3=E/(ONE-TWO*ENU)
      EG2=E/(ONE+ENU)
      EG=EG2/TWO
      EG3=THREE*EG
      ELAM=(EBULK3-EG2)/THREE
C
C     ELASTIC STIFFNESS
C
      DO 20 K1=1,NTENS
        DO 10 K2=1,NTENS
           DDSDDE(K2,K1)=0.0
10     CONTINUE
20   CONTINUE
C
      DO 40 K1=1,NDI
        DO 30 K2=1,NDI
           DDSDDE(K2,K1)=ELAM
30     CONTINUE
        DDSDDE(K1,K1)=EG2+ELAM
40   CONTINUE
      DO 50 K1=NDI+1,NTENS
        DDSDDE(K1,K1)=EG
50   CONTINUE
C
C    CALCULATE STRESS FROM ELASTIC STRAINS
C    计算试应力
      DO 70 K1=1,NTENS
        DO 60 K2=1,NTENS
           STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
60     CONTINUE
70   CONTINUE
      
C
C    RECOVER ELASTIC AND PLASTIC STRAINS
C    从状态变量中调用弹性应变、塑性应变、和等效塑性应变
DO 80 K1=1,NTENS
         EELAS(K1)=STATEV(K1)+DSTRAN(K1)
         EPLAS(K1)=STATEV(K1+NTENS)
80   CONTINUE
         EQPLAS=STATEV(1+2*NTENS)
      WRITE(6,1000) KINC,NOEL,NPT,EQPLAS
1000  FORMAT(//,30X,"INCREMENT",2X,I2,2X,'ELEMENT NUMBER',
     & 2X,I3,2X,'NPT',2X,I1,2X,'CURRENT EQ PLASTIC STRAIN',2X,F12.6)         
     
      !塑性材料行为                                                          不起作用？
        CALL SPRINC(S,PS,1,NDI,NSHR)

SMISES=(PS(1)+PS(2)+PS(3))**THREE/(PS(1)*PS(2)*PS(3))
IF (SMISES.GT.(1.0+TOLER)*KA) THEN
C     屈服后的流动方向
      I1=PS(1)+PS(2)+PS(3)
FLOW(1)=KB*(THREE*I1*I1/KB-STRESS(2)*STRESS(3)+STRESS(5)**TWO)
      FLOW(2)=KB*(THREE*I1*I1/KB-STRESS(3)*STRESS(1)+STRESS(6)**TWO)
      FLOW(3)=KB*(THREE*I1*I1/KB-STRESS(1)*STRESS(2)+STRESS(4)**TWO)
      FLOW(4)=KB*(TWO*STRESS(1)*STRESS(5)-TWO*STRESS(4)*STRESS(6))
      FLOW(5)=KB*(TWO*STRESS(2)*STRESS(6)-TWO*STRESS(5)*STRESS(4))
      FLOW(6)=KB*(TWO*STRESS(3)*STRESS(4)-TWO*STRESS(6)*STRESS(5))
      FS(1)=KA*(THREE*I1*I1/KA-STRESS(2)*STRESS(3)+STRESS(5)**TWO)
      FS(2)=KA*(THREE*I1*I1/KA-STRESS(3)*STRESS(1)+STRESS(6)**TWO)
      FS(3)=KA*(THREE*I1*I1/KA-STRESS(1)*STRESS(2)+STRESS(4)**TWO)
      FS(4)=KA*(TWO*STRESS(1)*STRESS(5)-TWO*STRESS(4)*STRESS(6))
      FS(5)=KA*(TWO*STRESS(2)*STRESS(6)-TWO*STRESS(5)*STRESS(4))
      FS(6)=KA*(TWO*STRESS(3)*STRESS(4)-TWO*STRESS(6)*STRESS(5))
      WRITE(6,4000) FS(1),PS(1),FLOW(1)
4000  FORMAT(//,3X,'FS',2X,I2,2X,'PS',2X,I2,2X,'FLOW',2X,I2,2X)
C     半隐式算法求等效塑性应变
      SYIELD=KA
DEQPL=0.0
DO 130 KEWTON=1,NEWTON
RHS=SMISES-SYIELD
RHD=EELAM*(FS(1)*FLOW(1)+FS(2)*FLOW(2)+FS(3)*FLOW(3))+
     &     ELAM*(FS(2)*FLOW(1)+FS(3)*FLOW(1)+FS(1)*FLOW(2)+FS(3)*FLOW(2)
     &    +FS(1)*FLOW(3)+FS(2)*FLOW(3))+EG*(FS(4)*FLOW(4)+FS(5)*FLOW(5)+
     &     FS(6)*FLOW(6))
DEQPL=DEQPL+RHS/RHD
      IF(ABS(RHS).LT.TOLER*KA)THEN
WRITE(6,2000) KINC,NOEL,NPT,DEQPL
2000  FORMAT(//,30X,"INCREMENT",2X,I2,2X,'ELEMENT NUMBER',
     & 2X,I3,2X,'NPT',2X,I1,2X,'EQ PLASTIC STRAIN INC',2X,F12.6)
      GOTO 140
ENDIF
130   CONTINUE
140   CONTINUE


C     计算应力更新应变
STRESS(1)=STRESS(1)-DEQPL*(EELAM*FLOW(1)+ELAM*(FLOW(2)
     & +FLOW(3)))
      STRESS(2)=STRESS(2)-DEQPL*(EELAM*FLOW(2)+ELAM*(FLOW(1)
     & +FLOW(3)))
      STRESS(3)=STRESS(3)-DEQPL*(EELAM*FLOW(3)+ELAM*(FLOW(1)
     &    +FLOW(2)))
  DO 150 K1=NDI+1,NTENS
STRESS(K1)=STRESS(K1)-DEQPL*EG*FLOW(K1)
150   CONTINUE
  DO 160 K1=1,NDI
EPLAS(K1)=EPLAS(K1)+DEQPL*FLOW(K1)
EELAS(k1)=EELAS(K1)-DEQPL*FLOW(K1)
160   CONTINUE
      EQPLAS=EQPLAS+DEQPL
      
C     更新雅克比矩阵
            EELAM=ELAM+EG2
EELAM2=EELAM*2
ELAMG=(ELAM+EG2)*EG


DO 170 K1=1,3     
BEITA=EELAM*FLOW(K1)*FLOW(K1)
170 CONTINUE
      DO 180 K2=4,6
BEITA=BEITA+EG*FLOW(K2)*FLOW(K2)
180 CONTINUE  
      DO 200 K1=1,3
DO 200 K2=1,3
DDSDDE(K1,K2)=EELAM2*FLOW(K1)*FLOW(K2)/BEITA
190   CONTINUE
200   CONTINUE
      DO 220 K1=1,3
DO 210 K2=4,6
DDSDDE(K1,K2)=ELAMG*FLOW(K1)*FLOW(K2)/BEITA
210   CONTINUE
220   CONTINUE
      DO 240 K1=4,6
DO 230 K2=3,1
DDSDDE(K1,K2)=DDSDDE(K2,K1)
230   CONTINUE
240   CONTINUE
      DO 260 K1=4,6
DO 250 K2=4,6
DDSDDE(K1,K2)=EG**2*FLOW(K1)*FLOW(K2)
250   CONTINUE
260   CONTINUE
       ENDIF


DO 310 K1=1,NTENS
STATEV(K1)=EELAS(K1)
STATEV(K1+NTENS)=EPLAS(K1)
310    CONTINUE
      STATEV(1+2*NTENS)=EQPLAS
WRITE(6,3000) KINC,NOEL,NPT,STATEV(1+2*NTENS)
3000  FORMAT(//,30X,"INCREMENT",2X,I2,2X,'ELEMENT NUMBER',
     & 2X,I3,2X,'NPT',2X,I1,2X,'TOTAL EQ PLASTIC STRAIN',2X,F12.6)
RETURN
END
      
      这是子程序，请大家帮忙看下为什么定义的塑性行为没有发生，

OptimusPrime

1# chenyaorock
结构有进入塑性阶段吗？

chenyaorock

你是说屈服了没？施加的力比屈服应力大许多了，应该是进入了吧 2# OptimusPrime