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
END DO c...
99 continue c...
ENDIF
C将弹性应变,塑性应变分量保存到状态变量中,并传到下一个增量步 C STORE ELASTIC AND (EQUIVALENT) PLASTIC STRAINS C IN STATE VARIABLE ARRAY C
DO K1=1,NTENS
STATEV(K1)=EELAS(K1)
STATEV(K1+NTENS)=EPLAS(K1) END DO
STATEV(1+2*NTENS)=EQPLAS C
RETURN END c...
c...子程序,根据等效塑性应变,利用插值的方法得到对应的屈服应力 SUBROUTINE HARDSUB(SYIELD,HARD,EQPLAS,TABLE,NVALUE) C
INCLUDE 'ABA_PARAM.INC' 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=ZERO
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 `, 1 `ENTERED IN ASCENDING ORDER') CALL XIT 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=DSYIEL/DEQPL
SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD GOTO 10 ENDIF END DO 10 CONTINUE ENDIF RETURN END
相关推荐: