通用有限元软件的二次开发.pptVIP

C MISES STRESSC SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2)) + (STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3)) + (STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1)) DO 90 K1=NDI+1,NTENS SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1) CONTINUE SMISES=SQRT(SMISES/TWO) IF (SMISES.GT.(1.0+TOLER1)*SYIEL0) THEN C FLOW DIRECTIONC SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE ONESY=ONE/SMISES DO 110 K1=1,NDI FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO) CONTINUE DO 120 K1=NDI+1,NTENS FLOW(K1)=STRESS(K1)*ONESY 120 CONTINUE C SOLVE FOR EQUIV STRESS, NEWTON ITERATION C SYIELD=SYIEL0 DEQPL=0.0 DO 130 KEWTON=1,NEWTON RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD) EQPLAS= EQPLAS+DEQPL SYIELD=SYIEL0+HARD*DEQPL IF(ABS(RHS).LT.TOLER2*SYIEL0) GOTO 140 CONTINUE WRITE(6,2) NEWTON 2 FORMAT(//,30X,***WARNING - PLASTICITY ALGORITHM DID NOT , CONVERGE AFTER ,I3, ITERATIONS) 140 CONTINUE C CALC STRESS AND UPDATE STRAINSC DO 150 K1=1,NDI STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO CONTINUE DO 160 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 CONTINUE EQPLAS=EQPLAS+DEQPL C JACOBIAN EFFHRD=EG3*HARD/(EG3+HARD) EFFG=EG*SYIELD/SMISES EFFG2=TWO*EFFG EFFG3=THREE*EFFG2/TWO EFFLAM=(EBULK3-EFFG2)/THREE DO 220 K1=1,NDI DO 210 K2=1,NDI