tangvp.for

本程序是D.R.J.OWEN主编的FINITE ELEMENTS IN PLASTICITYG一书中第八章例题的有限元程序

C$DEBUG
C$LARGE
	SUBROUTINE TANGVP(LPROP,STRES,PROPS,TIMEX,DTIME,
     .                  NSTRE,NTYPE,MMATS,NCRIT,DMATX) !(Page 286)
C*********************************************************************
C
C*** THE SUBROUTINE EVALUATES THE PSEUDO D-MATRIX
C
C*********************************************************************
	DIMENSION STRES(4),CMATX(4,4),TMATX(4,4),TRIX1(4,4),TRIX2(4,4),
     .          PROPS(MMATS,10),DEVIA(4),DMATX(4,4)
	ROOT3=1.73205080757
	FDATM=PROPS(LPROP,5)
	GAMMA=PROPS(LPROP,8)
	DELTA=PROPS(LPROP,9)
	NFLOW=PROPS(LPROP,10)
	CALL INVAR(DEVIA,LPROP,MMATS,NCRIT,PROPS,SINT3,STEFF,STRES,THETA,
     .           VARJ2,YIELD)
	FCURR=YIELD-FDATM
	FNORM=FCURR/FDATM
	IF(FNORM.LE.0.0) RETURN
	IF(NFLOW.EQ.1) GO TO 10
	CMULT=EXP(DELTA*FNORM)-1.0
	GRADP=DELTA*(EXP(DELTA*FNORM))/FDATM
	GO TO 20
   10 CMULT=FNORM*DELTA
      GRADP=DELTA*(FNORM**(DELTA-1.0))/FDATM
   20 FACT1=GAMMA*ROOT3*CMULT/(2.0*STEFF)
      FACT2=GAMMA*(0.75*GRADP/VARJ2-3.0*CMULT/(4.0*ROOT3*STEFF*VARJ2))
C
C*** MATRICES M1 AND M2 FOR A VON MISES MATERIAL
C
      TRIX1(1,1)=0.666666667
      TRIX1(1,2)=-0.333333333
      TRIX1(1,3)=0.0
      TRIX1(2,2)=0.666666667
      TRIX1(2,3)=0.0
      TRIX1(3,3)=2.0
	IF(NTYPE.NE.3) GO TO 30
      TRIX1(1,4)=-0.333333333
      TRIX1(2,4)=-0.333333333
      TRIX1(3,4)=0.0
      TRIX1(4,4)=0.666666667
   30 TRIX2(1,1)=DEVIA(1)*DEVIA(1)
      TRIX2(1,2)=DEVIA(1)*DEVIA(2)
      TRIX2(1,3)=2.0*DEVIA(1)*DEVIA(3)
      TRIX2(2,2)=DEVIA(2)*DEVIA(2)
      TRIX2(2,3)=2.0*DEVIA(2)*DEVIA(3)
      TRIX2(3,3)=4.0*DEVIA(3)*DEVIA(3)
	IF(NTYPE.NE.3) GO TO 40
      TRIX2(1,4)=DEVIA(1)*DEVIA(4)
      TRIX2(2,4)=DEVIA(2)*DEVIA(4)
      TRIX2(3,4)=2.0*DEVIA(3)*DEVIA(4)
      TRIX2(4,4)=DEVIA(4)*DEVIA(4)
   40 DO 50 ISTRE=1,NSTRE
      DO 50 JSTRE=1,NSTRE
	TRIX1(JSTRE,ISTRE)=TRIX1(ISTRE,JSTRE)
   50	TRIX2(JSTRE,ISTRE)=TRIX2(ISTRE,JSTRE)
      DO 60 ISTRE=1,NSTRE
      DO 60 JSTRE=1,NSTRE
   60 CMATX(ISTRE,JSTRE)=TIMEX*DTIME*(FACT1*TRIX1(ISTRE,JSTRE)
     . +FACT2*TRIX2(ISTRE,JSTRE))
	CALL INVERT(DMATX,TMATX,NSTRE)
      DO 70 ISTRE=1,NSTRE
      DO 70 JSTRE=1,NSTRE
   70 TMATX(ISTRE,JSTRE)=TMATX(ISTRE,JSTRE)+CMATX(ISTRE,JSTRE)
      CALL INVERT(TMATX,DMATX,NSTRE)
	RETURN
	END