验证子程序时当材料达到屈服时网格发生畸变

nicky_wang22 · 发表于 2020-6-17 11:01:16

我在验证子程序的时候，发现每当材料发生屈服时也就是应力应变曲线斜率为0时，网格发生畸变，我觉得会不会是此时应力增量变为0导致应力增量出了问题，这方面一直没搞懂，有没有大佬帮忙看看代码，多谢了

!C**********************************************************
!C VUMAT FOR ZWT MODEL WITH IMPACT ANALYSIS
!C WRITTEN BY ZZD
!C**********************************************************
SUBROUTINE VUMAT(
!C READ ONLY (UNMODIFIABLE)VARIABLES -
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 (MODIFIABLE) VARIABLES -
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),
4 DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR),
5 FIELDOLD(NBLOCK,NFIELDV), STRESSOLD(NBLOCK,NDIR+NSHR),
6 STATEOLD(NBLOCK,NSTATEV), ENERINTERNOLD(NBLOCK),
7 ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),
8 STRETCHNEW(NBLOCK,NDIR+NSHR),
8 DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),
9 FIELDNEW(NBLOCK,NFIELDV),
1 STRESSNEW(NBLOCK,NDIR+NSHR), STATENEW(NBLOCK,NSTATEV),
2 ENERINTERNNEW(NBLOCK), ENERINELASNEW(NBLOCK)
!C
CHARACTER*80 CMNAME
!C
!C STRAIN
!C STATE(*,1) = STRAIN11
!C STATE(*,2) = STRAIN22
!C STATE(*,3) = STRAIN33
!C STATE(*,4) = STRAIN12
!C STATE(*,5) = STRAIN23
!C STATE(*,6) = STRAIN31
!C
!C UNDAMAGED STRESS
!C STATE(*,7) = STRESS11 (UNDAMAGED)
!C STATE(*,8) = STRESS22 (UNDAMAGED)
!C STATE(*,9) = STRESS33 (UNDAMAGED)
!C STATE(*,10) = STRESS12 (UNDAMAGED)
!C STATE(*,11) = STRESS23 (UNDAMAGED)
!C STATE(*,12) = STRESS31 (UNDAMAGED)
!C
!C STATE(*,13) = DAMAGE
!C
PARAMETER( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0, THREE = 3.D0,
1 THIRD = ONE/THREE, HALF = ONE/TWO, TWOTHIRDS = TWO/THREE,
2 THREEHALFS = 1.5D0)
!C
!C
E0 = PROPS(1)
ALPHA = PROPS(2)
BETA = PROPS(3)
!c E1 = PROPS(4)
MU = PROPS(4)
!c THETA = PROPS(6)
!c E2 = PROPS(7)
!c THETA2 = PROPS(8)
!c E = PROPS(9)
TWOMU = TWO*MU
AF = ONE/((ONE+MU)*(ONE-TWOMU))
!C
!C
!C
DO I = 1, NBLOCK
!C CALCULATE DINCSIGMA=A*INCE(KL)
DINC_SIGMA11=AF*((ONE-MU)*STRAININC(I,1)+MU*STRAININC(I,2)
1 +MU*STRAININC(I,3))
DINC_SIGMA22=AF*(MU*STRAININC(I,1)+(ONE-MU)*STRAININC(I,2)
1 +MU*STRAININC(I,3))
DINC_SIGMA33=AF*(MU*STRAININC(I,1)+MU*STRAININC(I,2)
1 +(ONE-MU)*STRAININC(I,3))
DINC_SIGMA12 = AF*((ONE-TWOMU)/TWO)*STRAININC(I,4)
DINC_SIGMA23 = AF*((ONE-TWOMU)/TWO)*STRAININC(I,5)
DINC_SIGMA31 = AF*((ONE-TWOMU)/TWO)*STRAININC(I,6)
!C UPDATE STRAIN
STATENEW(I,1)=STATEOLD(I,1)+STRAININC(I,1)
STATENEW(I,2)=STATEOLD(I,2)+STRAININC(I,2)
STATENEW(I,3)=STATEOLD(I,3)+STRAININC(I,3)
STATENEW(I,4)=STATEOLD(I,4)+STRAININC(I,4)
STATENEW(I,5)=STATEOLD(I,5)+STRAININC(I,5)
STATENEW(I,6)=STATEOLD(I,6)+STRAININC(I,6)
!C CALCULATE INCREAMENT OF SIJ = INC_SIJ
!C
!C DINCF_S11=E0*DINC_SIGMA11+(TWO*ALPHA*STATENEW(I,1)+THREE
!C 1 *BETA*(STATENEW(I,1)**TWO))*DINC_SIGMA11
!C
!C DINCF_S22=E0*DINC_SIGMA22+(TWO*ALPHA*STATENEW(I,2)+THREE
!C 1 *BETA*(STATENEW(I,2)**TWO))*DINC_SIGMA22
!C
!C DINCF_S33=E0*DINC_SIGMA33+(TWO*ALPHA*STATENEW(I,3)+THREE
!C 1 *BETA*(STATENEW(I,3)**TWO))*DINC_SIGMA33
!C
!C DINCF_S12=E0*DINC_SIGMA12+(TWO*ALPHA*STATENEW(I,4)+THREE
!C 1 *BETA*(STATENEW(I,4)**TWO))*DINC_SIGMA12
!C
!C DINCF_S23=E0*DINC_SIGMA23+(TWO*ALPHA*STATENEW(I,5)+THREE
!C 1 *BETA*(STATENEW(I,5)**TWO))*DINC_SIGMA23
!C
!C DINCF_S31=E0*DINC_SIGMA31+(TWO*ALPHA*STATENEW(I,6)+THREE
!C 1 *BETA*(STATENEW(I,6)**TWO))*DINC_SIGMA31
!C CALCULATE LOW STRAIN RATE
!c DINCS_S11=STATENEW(I,7)*(EXP(-1.*DT/THETA)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA))/(DT+0.000000000000001))*E1*THETA)*DINC_SIGMA11
!c DINCS_S22=STATENEW(I,8)*(EXP(-1.*DT/THETA)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA))/(DT+0.000000000000001))*E1*THETA)*DINC_SIGMA22
!c DINCS_S33=STATENEW(I,9)*(EXP(-1.*DT/THETA)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA))/(DT+0.000000000000001))*E1*THETA)*DINC_SIGMA33
!c DINCS_S12=STATENEW(I,10)*(EXP(-1.*DT/THETA)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA))/(DT+0.000000000000001))*E1*THETA)*DINC_SIGMA12
!c DINCS_S23=STATENEW(I,11)*(EXP(-1.*DT/THETA)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA))/(DT+0.000000000000001))*E1*THETA)*DINC_SIGMA23
!c DINCS_S31=STATENEW(I,12)*(EXP(-1.*DT/THETA)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA))/(DT+0.000000000000001))*E1*THETA)*DINC_SIGMA31
!C CALCULATE HIGH STRAIN RATE PART
!c DINCT_S11=STATENEW(I,13)*(EXP(-1.*DT/THETA2)-ONE)+(((ONE-EXP(-1.
!c 1 *DT/THETA2))/(DT+0.000000000000001))*E2*THETA2)*DINC_SIGMA11
!c DINCT_S22=STATENEW(I,14)*(EXP(-1.*DT/THETA2)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA2))/(DT+0.000000000000001))*E2*THETA2)*DINC_SIGMA22
!c DINCT_S33=STATENEW(I,15)*(EXP(-1.*DT/THETA2)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA2))/(DT+0.000000000000001))*E2*THETA2)*DINC_SIGMA33
!c DINCT_S12=STATENEW(I,16)*(EXP(-1.*DT/THETA2)-ONE)+(((ONE-EXP(-1.*DT
!c 1 /THETA2))/(DT+0.000000000000001))*E2*THETA2)*DINC_SIGMA12
!c DINCT_S23=STATENEW(I,17)*(EXP(-1.*DT/THETA2)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA2))/(DT+0.000000000000001))*E2*THETA2)*DINC_SIGMA23
!c DINCT_S31=STATENEW(I,18)*(EXP(-1.*DT/THETA2)-ONE)+(((ONE-EXP(-1.*DT/
!c 1 THETA2))/(DT+0.000000000000001))*E2*THETA2)*DINC_SIGMA31
!C UPDATE THE HISTORY
!c STATENEW(I,7)=STATEOLD(I,7)+DINCS_S11
!c STATENEW(I,8)=STATEOLD(I,8)+DINCS_S22
!c STATENEW(I,9)=STATEOLD(I,9)+DINCS_S33
!c STATENEW(I,10)=STATEOLD(I,10)+DINCS_S12
!c STATENEW(I,11)=STATEOLD(I,11)+DINCS_S23
!c STATENEW(I,12)=STATEOLD(I,12)+DINCS_S31
!c STATENEW(I,13)=STATEOLD(I,13)+DINCT_S11
!c STATENEW(I,14)=STATEOLD(I,14)+DINCT_S22
!c STATENEW(I,15)=STATEOLD(I,15)+DINCT_S33
!c STATENEW(I,16)=STATEOLD(I,16)+DINCT_S12
!c STATENEW(I,17)=STATEOLD(I,17)+DINCT_S23
!c STATENEW(I,18)=STATEOLD(I,18)+DINCT_S31
STRESSNEW(I,1)=STRESSOLD(I,1)+E0*DINC_SIGMA11+
1 (TWO*ALPHA*STATENEW(I,1)+THREE
2 *BETA*(STATENEW(I,1)**TWO))*DINC_SIGMA11
!C
STRESSNEW(I,2)=STRESSOLD(I,2)+E0*DINC_SIGMA22+
1 (TWO*ALPHA*STATENEW(I,2)+THREE
2 *BETA*(STATENEW(I,2)**TWO))*DINC_SIGMA22
!C
STRESSNEW(I,3)=STRESSOLD(I,3)+E0*DINC_SIGMA33+
1 (TWO*ALPHA*STATENEW(I,3)+THREE
2 *BETA*(STATENEW(I,3)**TWO))*DINC_SIGMA33
!C
STRESSNEW(I,4)=STRESSOLD(I,4)+E0*DINC_SIGMA12+
1 (TWO*ALPHA*STATENEW(I,4)+THREE
2 *BETA*(STATENEW(I,4)**TWO))*DINC_SIGMA12
!C
STRESSNEW(I,5)=STRESSOLD(I,5)+E0*DINC_SIGMA23+
1 (TWO*ALPHA*STATENEW(I,5)+THREE
2 *BETA*(STATENEW(I,5)**TWO))*DINC_SIGMA23
!C
STRESSNEW(I,6)=STRESSOLD(I,6)+E0*DINC_SIGMA31+
1 (TWO*ALPHA*STATENEW(I,6)+THREE
2 *BETA*(STATENEW(I,6)**TWO))*DINC_SIGMA31
END DO
RETURN
END

复制代码

账号		自动登录	找回密码
密码			注册

账号		自动登录	找回密码
密码			立即注册

[二次开发] 验证子程序时当材料达到屈服时网格发生畸变

本帖子中包含更多资源