```txt SUBROUTINE LUMASS (COORD, INTGR, LNODS, MATNO, NCONM, NDIME, NDOFN, NELEM, NGAUM, NMATS, NNODE, NPOIN, NTYPE, PROPS, YMASS) ``` ```txt C ``` ```txt C ``` ```txt C *** CALCULATES LUMPED MASS FOR 4, 8 AND 9 NODED ELEMENT ``` ```txt C ``` ```txt C****************************************************************************************** ``` ```txt DIMENSION COORD(NPOIN,1),ELCOD(2,9),DIAGM(9),POSGP(4), . LNODS(NELEM,1),CARTD(2,9),SHAPE(9),WEIGP(4), . PROPS(NMATS,1),GPCOD(2,9),MATNO(1),YMASS(1), . EMASS(171),DERIV(2,9),INTGR(1) ``` ```txt C ``` ```csv REWIND 3 TWOPI=6.283185307179586 NEVAB=NNODE*NDOFN NTOTV=NPOIN*NDOFN DO 500 ITOTV =1,NTOTV ``` ```txt 500 YMASS(ITOTV)=0.0 CALL GAUSSQ (NGAUM, POSGP, WEIGP) DO 100 IELEM=1, NELEM DO 5 IEVAB=1, 171 ``` ```txt 5 EMASS(IEVAB)=0.0 IMASS=INTGR(IELEM) KGASP=0 TAREA=0.0 LPROP=MATNO(IELEM) THICK=PROPS(LPROP,3) RHOEL=PROPS(LPROP,4) DO 10 INODE=1,NNODE DIAGM(INODE)=0.0 LNODE=LNODS(IELEM,INODE) DO 10 IDIME=1,NDIME ELCOD(IDIME,INODE)=COORD(LNODE,IDIME) ``` ```csv 10 CONTINUE DO 70 IGAUS=1,NGAUM EXISP=POSGP(IGAUS) DO 70 JGAUS=1,NGAUM KGASP=KGASP+1 ETASP=POSGP(JGAUS) CALL SFR2 (DERIV,NNODE,SHAPE,EXISP,ETASP) CALL JACOB2 (CARTD,DERIV,DJACB,ELCOD,GPCOD,IELEM,KGASP,NNODE,SHAPE) ``` ```txt DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS) IF(NTYPE.EQ.1) DVOLU=DVOLU*THICK IF(NTYPE.EQ.3) DVOLU=DVOLU*TWOPI*GPCOD(1,KGASP) IF(IMASS.EQ.1) GO TO 210 DO 20 INODE=1,NNODE SHAPI=SHAPE(INODE) ``` ```txt 20 DIAGM(INODE)=DIAGM(INODE)+SHAPI*SHAPI*DVOLU TAREA=TAREA+DVOLU ``` ```matlab 210 IF(IMASS.EQ.2) GO TO 70 DVOLU=DVOLU*RHOEL IEVAB=1 KOUNT=NEVAB DO 30 INODE=1,NNODE SHAPI=SHAPE(INODE) DO 60 JNODE=INODE,NNODE DMASS=DVOLU*SHAPI*SHAPE(JNODE) EMASS(IEVAB)=EMASS(IEVAB)+DMASS JEVAB=IEVAB+KOUNT EMASS(JEVAB)=EMASS(JEVAB)+DMASS 60 IEVAB=IEVAB+2 KOUNT=KOUNT-2 ``` ```fortran IEVAB=JEVAB+1 MASS 65 30 CONTINUE MASS 66 70 CONTINUE MASS 67 C MASS 68 C*** WRITES CONSISTENT MASS MATRIX ON TAPE 3 MASS 69 C MASS 70 IF(IMASS.EQ.2) GO TO 200 MASS 71 WRITE(3) EMASS MASS 72 WRITE(6,90) (EMASS(I),I=1,171) MASS 73 200 IF(IMASS.EQ.1) GO TO 100 MASS 74 C MASS 75 C *** GENERATES LUMPED MASS MATRIX PROPORTIONAL TO DIAGONAL MASS 76 C MASS 77 SUMAS=0. MASS 78 DO 40 INODE=1,NNODE MASS 79 40 SUMAS=SUMAS+DIAGM(INODE) MASS 80 TAREA=TAREA*RHOEL MASS 81 SUMAS=TAREA/SUMAS MASS 82 DO 50 INODE=1,NNODE MASS 83 LNODE=LNODS(IELEM,INODE) MASS 84 IPOSN=(LNODE-1)*NDOFN MASS 85 DO 50 IDOFN=1,NDOFN MASS 86 IPOSN=IPOSN+1 MASS 87 YMASS(IPOSN)=YMASS(IPOSN)+DIAGM(INODE)*SUMAS MASS 88 50 CONTINUE MASS 89 90 FORMAT(2X,9E12.3) MASS 90 100 CONTINUE MASS 91 C MASS 92 C CONCENTRATED MASSES MASS 93 C MASS 94 IF(NCONM.EQ.0) RETURN MASS 95 WRITE(6,900) MASS 96 DO 520 ICONM=1,NNOM MASS 97 READ(5,910) IPOIN,XCMAS,YCMAS MASS 98 900 FORMAT(5X,19HCONCENTRATED MASSES) MASS 99 WRITE(6,910) IPOIN,XCMAS,YCMAS MASS 100 NPOSN=(IPOIN-1)*NDOFN+1 MASS 101 YMASS(NPOSN)=YMASS(NPOSN)+XCMAS MASS 102 NPOSN=NPOSN+1 MASS 103 YMASS(NPOSN)=YMASS(NPOSN)+YCMAS MASS 104 520 CONTINUE MASS 105 C WRITE(6,90) (YMASS(I),I=1,NTOTV) MASS 106 910 FORMAT(I5,2F10.3) MASS 107 RETURN MASS 108 END MASS 109 ``` # MASS 24 Sets indicator for mass matrix evaluation. INTGR(I) = 1 for the consistent mass matrix and INTGR(I) = 2 for the special lumped mass vector. # MASS 35-52 Evaluate the diagonal element of the consistent mass matrix DIAGM. # MASS 53-63 Evaluates the element consistent mass matrix. # MASS 72 Writes element consistent mass matrix on tape 3. # MASS 78-80 Evaluates ELMAS, the sum of the diagonal elements. # MASS 81 Determines the total element mass from the element volume TAREA and mass density RHOEL. MASS 83–89 Scales the diagonal terms using the factor TAREA/ELMAS to preserve element mass and assembles the result into diagonal mass vector YMASS. MASS 95–107 Reads the concentrated masses and assembles them into YMASS. # 10.6.17 Subroutine MODPS This subroutine evaluates the elasticity matrix and has been described earlier in Chapter 6. The only changes involved are given below. ```csv SUBROUTINE MODPS (DMATX ,LPROP ,NMATS ,NSTRE ,NTYPE ,PROPS ) MODP 1 C************************** C C ** ELASTICITY D MATRIX C C************************** DIMENSION DMATX(4,4),PROPS(NMATS,1) ``` # 10.6.18 Subroutine NODXYR It calculates $(r, z)$ coordinates from $(R, \Theta)$ coordinates for axisymmetric problems. If coordinates of midside nodes are not read, it evaluates them by linear interpolation. An almost identical subroutine was described in Chapter 6. ```fortran SUBROUTINE NODXYR (COORD, LNODS, NELEM, NNODE, NPOIN, NRADS, NTYPE) NODX 1 C*************** NODX 2 C NODX 3 C*** INTERPOLATION OF MIDSIDE AND CENTER NODES NODX 4 C NODX 5 C*************** NODX 6 DIMENSION COORD(NPOIN, 1), LNODS(NELEM, 1) NODX 7 C NODX 8 IF(NTYPE.NE.3.OR.NRADS.EQ.0) GO TO 40 NODX 9 C NODX 10 C*** CHANGE POLAR COORDINATES TO CARTISIAN NODX 11 DO 50 IPOIN=1, NPOIN NODX 12 RADDI=COORD(IPOIN, 1) NODX 13 THETA=COORD(IPOIN, 2) NODX 14 THETA=0.017453292*THETA NODX 15 COORD(IPOIN, 1)=RADDI*SIN(THETA) NODX 16 50 COORD(IPOIN, 2)=RADDI*COS(THETA) NODX 17 C NODX 18 40 IF(NNODE.EQ.4) RETURN NODX 19 C NODX 20 LNODE = NNODE - 1 NODX 21 DO 30 Ielem=1, NELEM NODX 22 C*** LOOP OVER EACH ELEMENT EDGE NODX 23 DO 20 INODE=1, NNODE, 2 NODX 24 IF(INODE.EQ.9) GO TO 20 NODX 25 C*** COMPUTE THE NODE NUMBER OF THE FIRST NODE NODX 26 NODST=LNODS(IELEM, INODE) NODX 27 IGASH=INODE+2 NODX 28 IF(IGASH.GT.LNODE) IGASH=1 NODX 29 C*** COMPUTE THE NODE NUMBER OF THE LAST NODE NODX 30 NODFN=LNODS(IELEM, IGASH) NODX 31 MIDPT=INODE+1 NODX 32 ``` ```fortran C*** COMPUTE THE NODE NUMBER OF THE INTERMEDIATE NODE NODX 33 NODMD=LNODS(IELEM,MIDPT) NODX 34 TOTAL=ABS(COORD(NODMD,1))+ABS(COORD(NODMD,2)) NODX 35 C*** IF THE COORDINATES OF THE INTERMEDIATE NODE ARE BOTH ZERO NODX 36 C INTERPOLATE BY A STRAIGHT LINE NODX 37 IF(TOTAL.GT.0.0) GO TO 20 NODX 38 KOUNT=1 NODX 39 10 COORD(NODMD,KOUNT)=(COORD(NODST,KOUNT)+COORD(NODFN,KOUNT))/2.0 NODX 40 KOUNT=KOUNT+1 NODX 41 IF(KOUNT.EQ.2) GO TO 10 NODX 42 20 CONTINUE NODX 43 30 CONTINUE NODX 44 RETURN NODX 45 END NODX 46 ``` # 10.6.19 Subroutine OUTDYN This routine writes out most of the output on the line printer and on various tapes for plotting purposes. It outputs the displacements and stresses every NOUTP steps. It also writes the displacement and stress histories of specified nodal and integration points at every NOUTP steps. The complete state of displacements is also written on tape 13 for a deformation plot. The complete state of the stresses is written on tape 4. The principal stresses and their directions are also calculated and output. ```txt SUBROUTINE OUTDYN (DISPL, DTIME, ISTEP, NDOFN, NELEM, NGAUS, OUTP 1 NGRQS, NOUTD, NOUTP, NPOIN, NPRQD, NREQD, OUTP 2 NREQS, NTYPE, STRSG, TDISP, VIVEL) OUTP 3 C******************************* OUTP 4 C C** OUTPUT ROUTINE OUTP 5 C C******************************* OUTP 6 OUTP 7 OUTP 8 DIMENSION STRSG(4,1), DISPL(1), NPRQD(1), STRSP(3), OUTP 9 VIVEL(5,1), TDISP(1), NGRQS(1) OUTP 10 NSTR1=4 OUTP 11 KSTEP=ISTEP OUTP 12 MGAUS=NELEM*NGAUS*NGAUS OUTP 13 IF(ISTEP.EQ.1) WRITE(10,925) OUTP 14 TTIME=TTIME+DTIME OUTP 15 C OUTP 16 C *** WRITES DISPLACEMENT HISTORY AT REQUESTED NODAL POINTS ON TAPE 10 OUTP 17 C *** AND STRESS HISTORY AT REQUESTED GAUSS POINTS AT EVERY NOUTD STEPSOUTP 18 C OUTP 19 KOUNT=0 OUTP 20 KOUTD=(ISTEP/NOUTD)*NOUTD OUTP 21 IF(KOUTD.NE.ISTEP) GO TO 510 OUTP 22 DO 500 IPOIN=1,NPOIN OUTP 23 DO 500 IREQD=1,NREQD OUTP 24 IF(IPOIN.NE.NPRQD(IREQD)) GO TO 500 OUTP 25 NPOSN=(IPOIN-1)*NDOFN+1 OUTP 26 NPOSM=NPOSN+1 OUTP 27 KOUNT=KOUNT+1 OUTP 28 DISPL(KOUNT)=TDISP(NPOSN) OUTP 29 KOUNT=KOUNT+1 OUTP 30 DISPL(KOUNT)=TDISP(NPOSM) OUTP 31 500 CONTINUE OUTP 32 WRITE(10,960) (DISPL(IKOUN), IKOUN=1,KOUNT), TTIME OUTP 33 ``` ```csv DO 520 IGAUS=1,MGAUS DO 520 IREQS=1,NREQS IF(IGAUS.NE.NGRQS(IREQS)) GO TO 520 WRITE(11,950) (STRSG(ISTR1,IGAUS),ISTR1=1,NSTR1) 520 CONTINUE 510 KOUTD=(KSTEP/NOUTP)*NOUTP IF(KOUTD.NE.KSTEP) RETURN XTIME=FLOAT(KSTEP)*DTIME WRITE(6,604) KSTEP,XTIME 604 FORMAT(//5X,28H DISPLACEMENTS AT TIME STEP ,I10,5X,5HTIME ,E20.11) C C *** REARRANGE DISPLACEMENT VECTOR C NODEI=0 DO 550 IPOIN=1,NPOIN DO 550 IDOFN=1,NDOFN NODEI=NODEI+1 DISPL(NODEI)=TDISP(NODEI) 550 CONTINUE C C*** OUTPUT DISPLACEMENTS C 925 FORMAT(5X,' DISPLACEMENTS ') WRITE(6,990) 990 FORMAT(/3(1X,'NNODE',3X,'X-DISP',6X,'Y-DISP',3X)/) DO 560 IPOIN=1,NPOIN,3 NGASI=NDOFN*IPOIN-1 NGASJ=NGASI+NDOFN NGASK=NGASJ+NDOFN MGASI=NGASI+1 MGASJ=NGASJ+1 MGASK=NGASK+1 JPOIN=IPOIN+1 KPOIN=JPOIN+1 C C *** WRITES DISPLACEMENTS ON TAPE 13 FOR DEFORMATION PLOT C WRITE(13,910) IPOIN ,(DISPL(IGASI),IGASI=NGASI,MGASI) IF(JPOIN.GT.NPOIN) GO TO 200 WRITE(13,910) JPOIN ,(DISPL(IGASJ),IGASJ=NGASJ,MGASJ) IF(KPOIN.GT.NPOIN) GO TO 200 WRITE(13,910) KPOIN ,(DISPL(IGASK),IGASK=NGASK,MGASK) 200 CONTINUE C C *** WRITES DISPLACEMENTS ON OUTPUT FILE C 560 WRITE(6,920) IPOIN,DISPL(NGASI),DISPL(MGASI), JPOIN,DISPL(NGASJ),DISPL(MGASJ), KPOIN,DISPL(NGASK),DISPL(MGASK) C C *** WRITES STRESSES ON OUTPUT FILE C WRITE(6,900) IF(NTYPE.NE.3) WRITE(6,970) 970 FORMAT(1HO,1X,4HG.P.,6X,9HXX-STRESS,5X,9HYY-STRESS,5X,9HXY-STRESS,OUTP .5X,9HZZ-STRESS,6X,8HMAX P.S.,6X,8HMIN P.S.,3X,5HANGLE,3X,6H P.S.)OUTP IF(NTYPE.EQ.3) WRITE(6,975) 975 FORMAT(1HO,1X,4HG.P.,6X,9HRR-STRESS,5X,9HZZ-STRESS,5X,9HRZ-STRESS,OUTP .5X,9HTT-STRESS,6X,8HMAX P.S.,6X,8HMIN P.S.,3X,5HANGLE,3X,6H P.S.)OUTP KGAUS=0 DO 570 IELEM=1,NELEM KELGS=0 WRITE(6,930) IELEM 930 FORMAT(1HO,5X,13HELEMENT NO. =,I5) ``` ```csv DO 570 IGAUS=1,NGAUS DO 570 JGAUS=1,NGAUS KGAUS=KGAUS+1 KELGS=KELGS+1 XGASH=(STRSG(1,KGAUS)+STRSG(2,KGAUS))*0.5 XGISH=(STRSG(1,KGAUS)-STRSG(2,KGAUS))*0.5 XGESH=STRSG(3,KGAUS) XGOSH=SQRT(XGISH*XGISH+XGESH*XGESH) STRSP(1)=XGASH+XGOSH STRSP(2)=XGASH-XGOSH IF(XGISH.EQ.0.0) XGISH=0.1E-20 STRSP(3)=ATAN(XGESH/XGISH)*28.647889757 C C *** WRITES COMPLETE STRESS STATE ON TAPE 4 C WRITE(4,950) (STRSG(ISTR1,KGAUS),ISTR1=1,NSTR1), .(STRSP(ISTRE),ISTRE=1,3) 570 WRITE(6,940) KELGS,(STRSG(ISTR1,KGAUS),ISTR1=1,NSTR1), .(STRSP(ISTRE),ISTRE=1,3),VIVEL(5,KGAUS) 980 FORMAT(1X,60I2) 960 FORMAT(1X,10E11.4) 950 FORMAT(7E10.4) 940 FORMAT(I5,2X,6E14.6,F8.3,E14.6) 900 FORMAT(/,10X,8HSTRESSES,/) 920 FORMAT(3(1X,I5,2E12.5)) 910 FORMAT(I5,2E15.6) RETURN END OUTP 98 OUTP 99 OUTP 100 OUTP 101 OUTP 102 OUTP 103 OUTP 104 OUTP 105 OUTP 106 OUTP 107 OUTP 108 OUTP 109 OUTP 110 OUTP 111 OUTP 112 OUTP 113 OUTP 114 OUTP 115 OUTP 116 OUTP 117 OUTP 118 OUTP 119 OUTP 120 OUTP 121 OUTP 122 OUTP 123 OUTP 124 OUTP 125 ``` # 10.6.20 Subroutine PREVOS This routine reads and write the initial forces and stresses. ```csv SUBROUTINE PREVOS (FORCE ,NDOFN ,NELEM ,NGAUS ,NPOIN ,NPREV , STRIN ) C*************** PREV 2 C C*** GRAVITY LOADS AND STRESSES PREV 3 C C*************** PREV 4 C DIMENSION FORCE(1) ,STRIN(4,1) PREV 5 C IF(NPREV.EQ.0) RETURN PREV 6 C NSTR1=4 PREV 7 NGAU2=NGAUS*NGAUS PREV 8 C C*** READ GRAVITY LOADS PREV 9 C WRITE(6,920) PREV 10 920 FORMAT(//4X,6H NODE ,17H GRAVITY X-LOAD: ,17H GRAVITY Y-LOAD: /) PREV 11 200 READ (5,900) NGASH,XGASH,YGASH PREV 12 900 FORMAT(I5,4F10.3) PREV 13 910 FORMAT(I10,4E18.5) PREV 14 NPOSN=(NGASH-1)*NDOFN+1 PREV 15 FORCE(NPOSN)=XGASH PREV 16 NPOSN=NPOSN+1 PREV 17 FORCE(NPOSN)=YGASH PREV 18 WRITE(6,910) NGASH,XGASH,YGASH PREV 19 IF (NGASH.NE.NPOIN) GO TO 200 PREV 20 C C*** READ GRAVITY STRESS PREV 21 C PREV 22 PREV 23 PREV 24 PREV 25 PREV 26 PREV 27 PREV 28 PREV 29 PREV 30 ``` WRITE(6,930) PREV 31 930 FORMAT(//2X,9HGAUSS PT.,17H GRAVITY X-STRESS,17H GRAVITY Y-STRESS,PREV 32 .18H GRAVITY XY-STRESS,17H GRAVITY Z-STRESS/) PREV 33 DO 500 IELEM=1,NELEM PREV 34 DO 500 IGAUS=1,NGAU2 PREV 35 READ(5,900) KGAUS,(STRIN(ISTRI,KGAUS),ISTRI=1,NSTR1) PREV 36 500 WRITE(6,910)KGAUS,(STRIN(ISTRI,KGAUS),ISTRI=1,NSTR1) PREV 37 RETURN PREV 38 END PREV 39 # 10.6.21 Subroutine RESVPL This routine evaluates the internal resisting force vector $$ \boldsymbol {p} _ {n} = \int_ {\Omega} [ \boldsymbol {B} ] _ {n} ^ {T} \boldsymbol {\sigma} _ {n} d \Omega . $$ It is very similar to the routine described in Section 8.8. SUBROUTINE RESVPL (COORD, DTIME, LNODS, MATNO, NCRIT, NDIME, RESD 1 NDOFN, NELEM, NGAUS, NLAPS, NNODE, NMATS, RESD 2 NPOIN, NSTRE, NTYPE, POSGP, PROPS, RESID, RESD 3 RLOAD, STRIN, STRSG, TDISP, VISTN, VIVEL, RESD 4 WEIGP) C****************************************************************************************** C C*** EVALUATION OF INTEGRAL (B)**T*(SIGMA) C C****************************************************************************************** DIMENSION COORD(NPOIN,1), DERIV(2,9), DJACM(2,2), AVECT(4), MATNO(1), RESD 11 PROPS(NMATS,1), DLCOD(2,9), STRIN(4,1), DEVIA(4), TDISP(1), RESD 12 LNODS(NELEM,1), GPCOD(2,9), STRSG(4,1), STRAN(4), POSGP(1), RESD 13 RLOAD(NELEM,1), CARTD(2,9), VISTN(4,1), STRES(4), WEIGP(1), RESD 14 DMATX(4,4), ELCOD(2,9), VIVEL(5,1), SHAPE(9), RESID(1), RESD 15 BMATX(4,18), ELDIS(2,9), DESTN(4) RESD 16 KGAUS=0 NSTR1=4 NEVAB=NNODE*NDOFN NTOTV=NPOIN*NDOFN TWOPI=6.283185307179586 DO 530 IELEM=1, NELEM DO 540 IEVAB=1, NEVAB 540 RLOAD(IELEM, IEVAB)=0.0 530 CONTINUE DO 510 ITOTV=1, NTOTV 510 RESID(ITOTV)=0.0 C C*** LOOP OVER ALL THE ELEMENTS C DO 20 IELEM=1, NELEM LPROP=MATNO(IELEM) THICK=PROPS(LPROP,3) POISS=PROPS(LPROP,2) FRICT=PROPS(LPROP,8) C C*** COMPUTE NEW COORDINATES AND DISPLACEMENTS OF THE C ELEMENT NODAL POINTS C DO 30 INODE =1, NNODE LNODE=IABS(LNODS(IELEM, INODE)) NPOSN=(LNODE-1)*NDOFN RESD 1 RESD 2 RESD 3 RESD 4 RESD 5 RESD 6 RESD 7 RESD 8 RESD 9 RESD 10 RESD 11 RESD 12 RESD 13 RESD 14 RESD 15 RESD 16 RESD 17 RESD 18 RESD 19 RESD 20 RESD 21 RESD 22 RESD 23 RESD 24 RESD 25 RESD 26 RESD 27 RESD 28 RESD 29 RESD 30 RESD 31 RESD 32 RESD 33 RESD 34 RESD 35 RESD 36 RESD 37 RESD 38 RESD 39 RESD 40 RESD 41 RESD 42 ```csv DO 30 IDOFN=1,NDOFN RESD 43 NPOSN=NPOSN+1 RESD 44 ELCOD(IDOFN,INODE)=COORD(LNODE,IDOFN) RESD 45 DLCOD(IDOFN,INODE)=COORD(LNODE,IDOFN)+TDISP(NPOSN) RESD 46 30 ELDIS(IDOFN,INODE)=TDISP(NPOSN) RESD 47 CALL MODPS (DMATX,LPROP,NMATS,NSTRE,NTYPE,PROPS) RESD 48 KGASP=0 RESD 49 DO 40 IGAUS=1,NGAUS RESD 50 DO 40 JGAUS=1,NGAUS RESD 51 KGAUS=KGAUS+1 RESD 52 KGASP=KGASP+1 RESD 53 EXISP=POSGP(IGAUS) RESD 54 ETASP=POSGP(JGAUS) RESD 55 C CALL SFR2 (DERIV,NNODE,SHAPE,EXISP,ETASP) RESD 57 CALL JACOB2 (CARTD,DERIV,DJACB,ELCOD,GPCOD, RESD 58 IELEM,KGASP,NNODE,SHAPE) RESD 59 CALL JACOBD (CARTD,DLCOD,DJACM,NDIME,NLAPS,NNODE) RESD 60 DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS) RESD 61 IF(NTYPE.EQ.1) DVOLU=DVOLU*THICK RESD 62 IF(NTYPE.EQ.3) DVOLU=DVOLU*TWOPI*GPCOD(1,KGASP) RESD 63 CALL BLARGE (BMATX,CARTD,DJACM,DLCOD,GPCOD, RESD 64 KGASP,NLAPS,NNODE,NTYPE,SHAPE) RESD 65 CALL LINGNL (CARTD,DJACM,DMATX,ELDIS,GPCOD,KGASP, RESD 66 KGAUS,NDOFN,NLAPS,NNODE,NSTRE,NTYPE, RESD 67 POISS,SHAPE,STRAN,STRES,VISTN) RESD 68 C DO 580 ISTR1=1,NSTR1 RESD 69 580 STRES(ISTR1)=STRES(ISTR1)+STRIN(ISTR1,KGAUS) RESD 70 DO 570 ISTR1=1,NSTR1 RESD 71 570 STRSG(ISTR1,KGAUS)=STRES(ISTR1) RESD 72 C IF(NLAPS.EQ.2.OR.NLAPS.EQ.0) GO TO 200 RESD 73 C CALL INVAR (DEVIA,LPROP,NCrit,NMATS,PROPS,SINT3,STEFF, RESD 74 STRES,THETA,VARJ2,YIELD) RESD 75 CALL YIELDF (AVECT,DEVIA,FRICT,NCrit,SINT3,STEFF,THETA,VARJ2) RESD 76 CALL FLOWVP (AVECT,KGAUS,LPROP,NCrit,NMATS,PROPS, RESD 77 STEFF,VIVEL,YIELD) RESD 78 C C*** VISCOPLASTIC STRAIN INCREMENT AND A MEASURE FOR HARDENING RESD 79 C DO 60 ISTR1=1,NSTR1 RESD 80 DESTN(ISTR1)=VIVEL(ISTR1,KGAUS)*DTIME RESD 81 60 VISTN(ISTR1,KGAUS)=VISTN(ISTR1,KGAUS)+DESTN(ISTR1) RESD 82 DEBAR=SQRT((2.0*(DESTN(1)*DESTN(1)+DESTN(2)*DESTN(2)+ RESD 83 DESTN(4)*DESTN(4))+DESTN(3)*DESTN(3))/3.0) RESD 84 VIVEL(5,KGAUS)=DEBAR RESD 85 C C*** COMPUT INT(B**T*SICMA) ON ELEMENT LEVEL RESD 86 C 200 CONTINUE RESD 94 KEVAB=0 RESD 95 DO 502 INODE=1,NNODE RESD 96 DO 502 IDOFN=1,NDOFN RESD 97 KEVAB=KEVAB+1 RESD 98 DO 501 ISTRE=1,NSTRE RESD 99 501 RLOAD(IELEM,KEVAB)=RLOAD(IELEM,KEVAB)+ .BMATX(ISTRE,KEVAB)*STRSG(ISTRE,KGAUS)*DVOLU RESD 100 502 CONTINUE RESD 101 40 CONTINUE RESD 102 20 CONTINUE RESD 103 C C*** ASSEMBLY OF RESID VECTOR RESD 104 C ```
CRESD 107
DO 500 IELEM=1,NELEMRESD 108
KEVAB=0RESD 109
DO 500 INODE=1,NNODERESD 110
LNODE=LNODS(IELEM,INODE)RESD 111
NPOSN=(LNODE-1)*NDOFNRESD 112
DO 500 IDOFN=1,NDOFNRESD 113
KEVAB=KEVAB+1RESD 114
NPOSN=NPOSN+1RESD 115
RESID(NPOSN)=RESID(NPOSN)+RLOAD(IELEM,KEVAB)RESD 116
500CONTINUERESD 117
RETURNRESD 118
ENDRESD 119
RESD 66-68Call LINGNL to determine the state of stress at the current Gauss point.
RESD 77-78Call INVAR to evaluate stress invariants at the current Gauss point.
RESD 79Call YIELDF to select the yield function and calculate the a vector.
RESD 80-81Call FLOWVP to define the rate of viscoplastic straining VIVEL if the stress point is outside the current yield surface.
RESD 86Evaluate the increments of viscoplastic strains DESTN.
RESD 87Evaluate the viscoplastic strains $(\epsilon_{vp})_{n+1}$ for the next time station $t_n + \Delta t$ , VISTN.
RESD 88-90Determine a measure of hardening for the current yield surface.
RESD 95-101Evaluate $p_n^{(e)}$ at the element level, RLOAD.
RESD 108-117Assemble $p_n$ , RESID.
# 10.6.22 Subroutine YIELDF This subroutine selects the yield function and calculates the vector a (AVECT) and is almost identical to the version described in Section 7.8.4.1. ```fortran SUBROUTINE YIELDF (AVECT, DEVIA, FRICT, NCRIT, SINT3, STEFF, YELD 1 THETA, VARJ2) YELD 2 C********** YELD 3 C YELD 4 C *** SELECTS YIELD FUNCTION AND CALCULATES VECTOR 'AVECT' YELD 5 C YELD 6 C********** YELD 7 DIMENSION AVECT(4), DEVIA(4), VECA1(4), VECA2(4), VECA3(4) YELD 8 IF(STEFF.EQ.0.0) RETURN YELD 9 NSTR1=4 YELD 10 TANTH=TAN(THETA) YELD 11 SINTH=SIN(THETA) YELD 12 COSTH=COS(THETA) YELD 13 COST3=COS(3.0*THETA) YELD 14 ROOT3=1.73205080757 YELD 15 ``` ```txt C*** CALCULATE VECTOR A1 VECA1(1)=1.0 VECA1(2)=1.0 VECA1(3)=0.0 VECA1(4)=1.0 C*** CALCULATE VECTOR A2 DO 10 ISTR1=1,NSTR1 10 VECA2(ISTR1)=DEVIA(ISTR1)/(2.0*STEFF) VECA2(3)=DEVIA(3)/STEFF C*** CALCULATE VECTOR A3 VECA3(1)=DEVIA(2)*DEVIA(4)+VARJ2/3.0 VECA3(2)=DEVIA(1)*DEVIA(4)+VARJ2/3.0 VECA3(3)=-2.0*DEVIA(3)*DEVIA(4) VECA3(4)=DEVIA(1)*DEVIA(2)-DEVIA(3)*DEVIA(3)+VARJ2/3.0 GO TO (1,2,3,4) NCRIT C*** TRESCA 1 CONS1=0.0 ABTHE=ABS(THETA*57.29577951308) IF(ABTHE.LT.29.0) GO TO 20 CONS2=ROOT3 CONS3=0.0 GO TO 40 20 CONS2=2.0*(COSTH+SINTH*TAN(3.0*THETA)) CONS3=ROOT3*SINTH/(VARJ2*COST3) GO TO 40 C*** VON MISES 2 CONS1=0.0 CONS2=ROOT3 CONS3=0.0 GO TO 40 C*** MOHR-COULOMB 3 CONS1=SIN(FRICT*0.017453292)/3.0 ABTHE=ABS(THETA*57.29577951308) IF(ABTHE.LT.29.0) GO TO 30 CONS3=0.0 PLUMI=1.0 IF(THETA.GT.0.0) PLUMI=-1.0 CONS2=0.5*(ROOT3+PLUMI*CONS1/ROOT3) GO TO 40 30 TANT3=TAN(3.0*THETA) CONS2=COSTH*(1.0+TANTH*TANT3)+CONS1*(TANT3-TANTH)/ROOT3) CONS3=(ROOT3*SINTH+CONS1*COSTH)/(2.0*VARJ2*COST3) GO TO 40 C*** DRUCKER-PRAGER 4 SNPHI=SIN(FRICT*0.017453292) CONS1=2.0*SNPHI/(ROOT3*(3.0-SNPHI)) CONS2=1.0 CONS3=0.0 40 CONTINUE DO 50 ISTR1=1,NSTR1 50 AVECT(ISTR1)=CONS1*VECA1(ISTR1)+CONS2* .VECA2(ISTR1)+CONS3*VECA3(ISTR1) RETURN END YELD 16 YELD 17 YELD 18 YELD 19 YELD 20 YELD 21 YELD 22 YELD 23 YELD 24 YELD 25 YELD 26 YELD 27 YELD 28 YELD 29 YELD 30 YELD 31 YELD 32 YELD 33 YELD 34 YELD 35 YELD 36 YELD 37 YELD 38 YELD 39 YELD 40 YELD 41 YELD 42 YELD 43 YELD 44 YELD 45 YELD 46 YELD 47 YELD 48 YELD 49 YELD 50 YELD 51 YELD 52 YELD 53 YELD 54 YELD 55 YELD 56 YELD 57 YELD 58 YELD 59 YELD 60 YELD 61 YELD 62 YELD 63 YELD 64 YELD 65 YELD 66 YELD 67 YELD 68 YELD 69 ``` # 10.7 Examples # 10.7.1 Introduction To illustrate the use of DYNPAK we now describe the nonlinear transient dynamic analysis of (i) a spherical shell and (ii) a concrete gravity dam.