发布时间 : 星期三 文章一起学习更新完毕开始阅读
? ? ? ??2G ??? ???2G ????? ? ??2G ?J???
G ??? ? G?? G??? ?注意，在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) 读取等效塑性应变
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 c...
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