```txt C FACTS =FUNCTS (AZERO,BZERO,DTEND,DTIME,IFUNC,ISTEP,OMEGA) FACTH =FUNCTA (ACCEH,AFACT,DTEND,DTIME,IFUNC,ISTEP) FACTV =FUNCTA (ACCEV,AFACT,DTEND,DTIME,IFUNC,ISTEP) WRITE(6,910) FACTS,FACTH,FACTV 650 CONTINUE IF(ISTEP.EQ.1) GO TO 640 C C *** CALCULATES DAMPING AND K-STAR MATRICES C DO 530 ISIZE=1,NSIZE IMAXA=MAXAI(ISIZE) KMAXA=MAXAI(ISIZE+1)-1 JMAXA=MAXAJ(ISIZE) DO 530 LMAXA=IMAXA,KMAXA DAMPI(JMAXA)=AALFA*XMASS(LMAXA) 530 JMAXA=JMAXA+1 DO 580 IWKTL=1,NWKTL 580 DAMPI(IWKTL)=DAMPI(IWKTL)+BEETA*STIFF(IWKTL) CALL MULTPY (VELOL ,DAMPI ,VELOT ,MAXAJ ,NSIZE ,NWKTL ) KOUNT=(ISTEP/KSTEP)*KSTEP IF(KOUNT.NE.ISTEP) GO TO 660 640 DO 610 IWMTL=1,NWMTL 610 DAMPG(IWMTL)=CONSE*XMASS(IWMTL) DO 620 ISIZE=1,NSIZE IMAXA=MAXAI(ISIZE) 620 DAMPG(IMAXA)=DAMPG(IMAXA)-CONSH*YMASS(ISIZE) DO 630 IWMTL=1,NWMTL DAMPG(IWMTL)=DAMPG(IWMTL)+CONSG*STIFI(IWMTL) 630 STIFS(IWMTL)=STIFI(IWMTL)+DAMPG(IWMTL)*CONSF WRITE(6,900) (STIFS(I),I=1,NWMTL) CALL DECOMP (STIFS ,MAXAI ,NSIZE ,ISHOT ) C C *** CALCULATES PARTIAL EFFECTIVE LOAD VECTOR C 660 DO 520 ISIZE=1,NSIZE IF(IFUNC.NE.0) GO TO 570 IF(IFIXD.EQ.2) DISPL(ISIZE)=-VELOL(ISIZE)-FACTH*ACCEJ(ISIZE) +RLOAD(ISIZE) IF(IFIXD.EQ.1) DISPL(ISIZE)=-VELOL(ISIZE)-FACTV*ACCEK(ISIZE) +RLOAD(ISIZE) IF(IFIXD.EQ.0) DISPL(ISIZE)=-VELOL(ISIZE)-FACTH*ACCEJ(ISIZE) +RLOAD(ISIZE)-FACTV*ACCEK(ISIZE) IF(IFUNC.EQ.0) GO TO 520 570 DISPL(ISIZE)=-VELOL(ISIZE)+RLOAD(ISIZE)*FACTS 520 CONTINUE RETURN END IMEX 95 IMEX 96 IMEX 97 IMEX 98 IMEX 99 IMEX 100 IMEX 101 IMEX 102 IMEX 103 IMEX 104 IMEX 105 IMEX 106 IMEX 107 IMEX 108 IMEX 109 IMEX 110 IMEX 111 IMEX 112 IMEX 113 IMEX 114 IMEX 115 IMEX 116 IMEX 117 IMEX 118 IMEX 119 IMEX 120 IMEX 121 IMEX 122 IMEX 123 IMEX 124 IMEX 125 IMEX 126 IMEX 127 IMEX 128 IMEX 129 IMEX 130 IMEX 131 IMEX 132 IMEX 133 IMEX 134 IMEX 135 IMEX 136 IMEX 137 IMEX 138 IMEX 139 IMEX 140 IMEX 141 IMEX 142 ``` # 11.5.11 Subroutine ITRATE This routine generates the total effective load vector and solves for the incremental displacements. It then checks for convergence. ```txt SUBROUTINE ITRATE (ACCEI, ACCEL, CONSD, CONSF, XMASS, DISPI, ITER 1 . DISPL, DISPT, MAXAI, NCHEK, NSIZE, NWMTL, ITER 2 . RESID, STIFS, TOLER, VELOI, VELOL, VELOT, ITER 3 . IITER, MITER) ITER 4 C********** C C C *** CALCULATES INCREMENT IN DISPLACEMENT AND APPLIES CONVERGENCE C C********** C DIMENSION DISPI(1), VELOI(1), ACCEI(1), RESID(1), MAXAI(1), ITER 10 ITER 5 ITER 6 ITER 7 ITER 8 ITER 9 ITER 10 ``` ```fortran : DISPL(1), VELOL(1), ACCEL(1), STIFS(1), DISPT(1), ITER 11 : XMASS(1), VELOT(1) ITER 12 C C NCHEK=0 CALL MULTPY (ACCEL, XMASS, ACCEI, MAXAI, NSIZE, NWMTL) ITER 13 C C *** CALCULATES TOTAL EFFECTIVE LOAD VECTOR C DO 660 ISIZE=1, NSIZE 660 ACCEL(ISIZE)=DISPL(ISIZE)-ACCEL(ISIZE)-RESID(ISIZE) ITER 14 C C *** CALCULATES DELTA DISPLACEMENT C 210 CALL REDBAK (STIFS, ACCEL, MAXAI, NSIZE) ITER 15 C C *** APPLIES CONVERGENCE C SUMPP=0. SUMPQ=0. DO 670 ISIZE=1, NSIZE DISPP=ACCEL(ISIZE) DISPQ=DISPT(ISIZE)+DISPP DISPT(ISIZE)=DISPQ SUMPP=SUMPP+DISPP*DISPP SUMPQ=SUMPQ+DISPQ*DISPQ 670 CONTINUE DO 530 ISIZE=1, NSIZE ACCEI(ISIZE)=CONSF*(DISPT(ISIZE)-DISPI(ISIZE)) ITER 16 530 VELOT(ISIZE)=VELOI(ISIZE)+CONSD*ACCEI(ISIZE) ITER 17 220 SUMPP=SQRT(SUMPP/SUMPQ) IF(SUMPP.GT.TOLER) GO TO 550 NCHEK=1 GO TO 240 550 IF(IITER.LT.MITER) GO TO 230 240 DO 540 ISIZE=1, NSIZE VELOI(ISIZE)=VELOT(ISIZE) ITER 18 540 DISPI(ISIZE)=DISPT(ISIZE) ITER 19 230 CONTINUE RETURN END ITER 20 ITER 21 ITER 22 ITER 23 ITER 24 ITER 25 ITER 26 ITER 27 ITER 28 ITER 29 ITER 30 ITER 31 ITER 32 ITER 33 ITER 34 ITER 35 ITER 36 ITER 37 ITER 38 ITER 39 ITER 40 ITER 41 ITER 42 ITER 43 ITER 44 ITER 45 ITER 46 ITER 47 ITER 48 ITER 49 ITER 50 ITER 51 ``` ITER 20-21 Calculates total effective load vector. ITER 25 Solves for incremental displacements. ITER 28-37 Calculates norm of displacement increments. ITER 38-40 Calculates new and total displacement, velocities and accelerations. ITER 41-42 Applies convergence check. ITER 46–49 Stores the final velocities and displacements in vectors VELOI and DISPI respectively. # 11.5.12 Subroutine LINKIN This routine calculates the equation number from the array IFPRE which stores the information about the restrained degrees of freedom. ```fortran SUBROUTINE LINKIN (FORCE, IFPRE, INTGR, LEQNS, LNODS, MAXAI, LINK 1 MAXAJ, MHIGH, NDOFN, NELEM, NEQNS, NNODE, LINK 2 NPOIN, NWKTL, NWMTL, XMASS, YMASS) LINK 3 C C C *** LINKS WITH PROFILE SOLVER C C C DIMENSION LNODS(NELEM, 1), XMASS(1), MAXAI(1), INTGR(1), LINK 9 IFPRE(NDOFN, 1), YMASS(1), MAXAJ(1), MHIGH(1), LINK 10 LEQNS(18, 1), FORCE(1), EMASS(171) LINK 11 C IMASS=1 REWIND 3 NEVAB=NNODE*NDOFN C C C NUMBER OF UNKNOWNS C NEQNS=0 DO 100 IPOIN=1, NPOIN DO 150 IDOFN=1, NDOFN IF(IFPRE(IDOFN, IPOIN)) 110, 120, 110 NEQNS=NEQNS+1 IFPRE(IDOFN, IPOIN)=NEQNS GO TO 150 IFPRE(IDOFN, IPOIN)=0 CONTINUE WRITE(6, 7) IPOIN, (IFPRE(IDOFN, IPOIN), IDOFN=1, NDOFN) CONTINUE MEQNS=1+NEQNS C C C C C CONNECTIVITY ARRAY LEQNS C DO 70 IELEM=1, NELEM DO 70 IEVAB=1, NEVAB LEQNS(IEVAB, IELEM)=0 DO 50 IELEM=1, NELEM IEVAB=1 DO 80 INODE=1, NNODE IDENT=LNODS(IELEM, INODE) DO 80 IDOFN=1, NDOFN LEQNS(IEVAB, IELEM)=IFPRE(IDOFN, IDENT) IEVAB=IEVAB+1 WRITE(6, 6) IELEM, (LEQNS(IEVAB, IELEM), IEVAB=1, NEVAB) CONTINUE FORMAT(I10, 24I3) FORMAT(4I10) FORMAT(8E12.4) C C C C LOOP OVER ALL ELEMENTS C 250 DO 190 IELEM=1, NELEM IF(INTGR(IELEM).NE.IMASS) GO TO 190 CALL COLMHT (MHIGH, NEVAB, LEQNS(1, IELEM)) CONTINUE ADDRESES OF DIAGONAL ELEMENTS - MAXA ARRAY CALL ADDRES(MAXAJ, MHIGH, NEQNS, NWKTL, MKOUN) IF(IMASS.EQ.2) GO TO 205 DO 580 IEQNS=1, MEQNS MAXAI(IEQNS)=MAXAJ(IEQNS) IMASS=2 NWMTL=NWKTL LINK 1 LINK 5 LINK 6 LINK 7 LINK 8 LINK 9 LINK 10 LINK 11 LINK 12 LINK 13 LINK 14 LINK 15 LINK 16 LINK 17 LINK 18 LINK 19 LINK 20 LINK 21 LINK 22 LINK 23 LINK 24 LINK 25 LINK 26 LINK 27 LINK 28 LINK 29 LINK 30 LINK 31 LINK 32 LINK 33 LINK 34 LINK 35 LINK 36 LINK 37 LINK 38 LINK 39 LINK 40 LINK 41 LINK 42 LINK 43 LINK 44 LINK 45 LINK 46 LINK 47 LINK 48 LINK 49 LINK 50 LINK 51 LINK 52 LINK 53 LINK 54 LINK 55 LINK 56 LINK 57 LINK 58 LINK 59 LINK 60 LINK 61 LINK 62 LINK 63 LINK 64 ``` ```csv GO TO 250 205 CONTINUE WRITE(6,920) NEQNS,NWMTL,NWKTL WRITE(6,930) (MAXAI(I),I=1,MEQNS) WRITE(6,930) (MAXAJ(I),I=1,MEQNS) 930 FORMAT(5X,20I5) 920 FORMAT(/5X,'NEQNS=',I5,5X,'NWMTL=',I5,5X,'NWKTL=',I5/) IF(NWKTL.GT.6000) GO TO 210 GO TO 220 210 WRITE(6,910) STOP 220 CONTINUE 910 FORMAT (/'SET DIMENSION EXCEEDED - CHECK LINKIN '/) C C*** GLOBAL MASS MATRIX C DO 500 IELEM=1,NELEM IMASS=INTGR(IELEM) IF(IMASS.EQ.2) GO TO 500 READ (3) EMASS CALL ADDBAN (XMASS,MAXAI,EMASS,LEQNS(1,IELEM),NEVAB) 500 CONTINUE C C*** GLOBAL MASS VECTOR C NPOSM=0 DO 510 IPOIN =1,NPOIN DO 510 IDOFN =1,NDOFN NPOSM=NPOSM+1 NPOSN=IFPRE(IDOFN,IPOIN) IF(NPOSN.EQ.0) GO TO 510 YMASS(NPOSN)=YMASS(NPOSM) FORCE(NPOSN)=FORCE(NPOSM) 510 CONTINUE RETURN END LINK 65 LINK 66 LINK 67 LINK 68 LINK 69 LINK 70 LINK 71 LINK 72 LINK 73 LINK 74 LINK 75 LINK 76 LINK 77 LINK 78 LINK 79 LINK 80 LINK 81 LINK 82 LINK 83 LINK 84 LINK 85 LINK 86 LINK 87 LINK 88 LINK 89 LINK 90 LINK 91 LINK 92 LINK 93 LINK 94 LINK 95 LINK 96 LINK 97 LINK 98 LINK 99 LINK 100 ``` LINK 18-29 Reassigns IFPRE vector with equation numbers. If IFPRE is not zero than IFPRE is reassigned as zero. LINK 34-45 Evaluates the vector LEQNS on element level for assigning equation number corresponding to each node in an element. LINK 52-55 Calculates column height above the diagonal in global matrix. LINK 59–62 Assigns location for diagonal elements in global matrix. LINK 80-85 IMASS = 1 calculates stiffness matrix for only implicit elements. IMASS = 2 calculates stiffness matrix for complete mesh. # 11.5.13 Subroutine MULTPY This routine $^{(9)}$ evaluates the product of square matrix AMATX and an array START and stores the result in FINAL. ```txt SUBROUTINE MULTPY (FINAL, AMATX, START, MAXAI, NEQNS, NWMTL) MULT 1 C*************** MULT 2 C MULT 3 C *** TO EVALUATE PRODUCT OF B TIMES RR AND STORE RESULT IN TT C MULT 4 C*************** MULT 5 DIMENSION FINAL(1), AMATX(1), START(1), MAXAI(1) MULT 6 C MULT 7 MULT 8 ``` ```txt IF(NWMTL.GT.NEQNS) GO TO 20 DO 10 IEQNS=1,NEQNS 10 FINAL(IEQNS)=AMATX(IEQNS)*START(IEQNS) RETURN C 20 DO 40 IEQNS=1,NEQNS 40 FINAL(IEQNS)=0.0 DO 100 IEQNS=1,NEQNS LOWER=MAXAI(IEQNS) KUPER=MAXAI(IEQNS+1)-1 JEQNS=IEQNS+1 TERMI=START(IEQNS) DO 100 ICOLM=LOWER,KUPER JEQNS=JEQNS-1 100 FINAL(JEQNS)=FINAL(JEQNS)+AMATX(ICOLM)*TERMI IF(NEQNS.EQ.1) RETURN DO 200 IEQNS=2,NEQNS LOWER=MAXAI(IEQNS)+1 KUPER=MAXAI(IEQNS+1)-1 IF(KUPER-LOWER) 200,210,210 210 JEQNS=IEQNS SUMAA=0.0 DO 220 ICOLM=LOWER,KUPER JEQNS=JEQNS-1 220 SUMAA=SUMAA+AMATX(ICOLM)*START(JEQNS) FINAL(IEQNS)=FINAL(IEQNS)+SUMAA 200 CONTINUE RETURN END MULT 9 MULT 10 MULT 11 MULT 12 MULT 13 MULT 14 MULT 15 MULT 16 MULT 17 MULT 18 MULT 19 MULT 20 MULT 21 MULT 22 MULT 23 MULT 24 MULT 25 MULT 26 MULT 27 MULT 28 MULT 29 MULT 30 MULT 31 MULT 32 MULT 33 MULT 34 MULT 35 MULT 36 MULT 37 ``` # 11.5.14 Subroutine REDBAK This routine $^{(9)}$ solves the equations after the matrix is decomposed (into the form $LDL^{T}$ ) using forward and backward substitution. ```txt SUBROUTINE REDBAK (STIFF ,FORCE ,MAXAI ,NEQNS ) C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C DIMENSION STIFF(1) ,FORCE(1) ,MAXAI(1) C DO 400 IEQNS=1,NEQNS LOWER =MAXAI(IEQNS)+1 KUPER=MAXAI(IEQNS+1)-1 IF(KUPER-LOWER) 400,410,410 410 JEQNS=IEQNS SUMCC=0.0 DO 420 ICOLM=LOWER,KUPER JEQNS=JEQNS-1 420 SUMCC=SUMCC+STIFF(ICOLM)*FORCE(JEQNS) FORCE(IEQNS)=FORCE(IEQNS)-SUMCC 400 CONTINUE C DO 480 IEQNS=1,NEQNS KMAXA=MAXAI(IEQNS) 480 FORCE(IEQNS)=FORCE(IEQNS)/STIFF(KMAXA) IF(NEQNS.EQ.1) RETURN JEQNS=NEQNS DO 500 IEQNS=2,NEQNS LOWER=MAXAI(JEQNS)+1 KUPER=MAXAI(JEQNS+1)-1 RBAK 1 RBAK 2 RBAK 3 RBAK 4 RBAK 5 RBAK 6 RBAK 7 RBAK 8 RBAK 9 RBAK 10 RBAK 11 RBAK 12 RBAK 13 RBAK 14 RBAK 15 RBAK 16 RBAK 17 RBAK 18 RBAK 19 RBAK 20 RBAK 21 RBAK 22 RBAK 23 RBAK 24 RBAK 25 RBAK 26 RBAK 27 RBAK 28 ```
IF(KUPER-LOWER) 500,510,510RBAK29
510KEQNS=JEQNSRBAK30
DO 520 ICOLM=LOWER,KUPERRBAK31
KEQNS=KEQNS-1RBAK32
520FORCE(KEQNS)=FORCE(KEQNS)-STIFF(ICOLM)*FORCE(JEQNS)RBAK33
500JEQNS=JEQNS-1RBAK34
RETURNRBAK35
ENDRBAK36
# 11.5.15 Subroutine RESEPL This routine evaluates the internal force vector for elasto-plastic materials. (See Section 7.8.7.) ```fortran SUBROUTINE RESEPL (COORD ,DISPL ,EFFST ,ELOAD ,EPSTN ,IITER , RESD 1 INTGR ,LEQNS ,LNODS ,MATNO ,NCRIT ,NDIME , RESD 2 NDOFN ,NELEM ,NGAUS ,NLAPS ,NMATS ,NNODE , RESD 3 NPOIN ,NSTRE ,NTYPE ,POSGP ,PROPS ,RESID , RESD 4 STRAG ,STRIN ,STRSG ,WEIGP ,IPRED ,ISTEP ) RESD 5 C****************************************************************************************** RESD 6 C RESD 7 C *** EVALUATES RESIDUAL FORCES RESD 8 C RESD 9 C****************************************************************************************** RESD 10 DIMENSION COORD(NPOIN,1),DERIV(2,9),DMATX(4, 4),AVECT(4),MATNO(1), RESD 11 PROPS(NMATS,1),DLCOD(2,9),BMATX(4,18),DEVIA(4),DISPL(1), RESD 12 LNODS(NELEM,1),GPCOD(2,9),DJACM(2, 2),STRAN(4),POSGP(1), RESD 13 ELOAD(NELEM,1),CARTD(2,9),SHAPE( 9),STRES(4),WEIGP(1), RESD 14 STRIN( 4,1),ELCOD(2,9),SIGMA( 4),SGTOT(4),EFFST(1), RESD 15 STRSG( 4,1),ELDIS(2,9),DESIG( 4),DVECT(4),EPSTN(1), RESD 16 STRAG( 4,1),RESID( 1),LEQNS(18,1),INTGR(1) RESD 17 TWOPI=6.283185307179586 RESD 18 NEVAB=NNODE*NDOFN RESD 19 NTOTV=NPOIN*NDOFN RESD 20 NSTR1=4 RESD 21 DO 530 IELEM=1,NELEM RESD 22 IF(INTGR(IELEM).EQ.2.AND.IITER.GT.1.AND.IPRED.EQ.1) GO TO 530 RESD 23 DO 540 IEVAB=1,NEVAB RESD 24 540 ELOAD(IELEM,IEVAB)=0.0 RESD 25 530 CONTINUE RESD 26 DO 510 ITOTV=1,NTOTV RESD 27 510 RESID(ITOTV)=0.0 RESD 28 KGAUS=0 RESD 29 DO 20 IELEM=1,NELEM RESD 30 IF(INTGR(IELEM).EQ.2.AND.IITER.GT.1.AND.IPRED.EQ.1) GO TO 20 RESD 31 LPROP=MATNO(IELEM) RESD 32 YOUNG=PROPS(LPROP,1) RESD 33 POISS=PROPS(LPROP,2) RESD 34 THICK=PROPS(LPROP,3) RESD 35 UNIAX=PROPS(LPROP,6) RESD 36 HARDS=PROPS(LPROP,7) RESD 37 FRICT=PROPS(LPROP,8) RESD 38 FRICT=FRICT*0.017453292 RESD 39 IF(NCRIT.EQ.3) UNIAX=UNIAX*COS(FRICT) RESD 40 IF(NCRIT.EQ.4) UNIAX=6.0*UNIAX*COS(FRICT)/ RESD 41 (1.73205080757*(3.0-SIN(FRICT))) RESD 42 C RESD 43 C*** COMPUTE COORDINATE AND INCREMENTAL DISPLACEMENTS OF THE RESD 44 C ELEMENT NODAL POINTS RESD 45 C RESD 46 IPOSN=0 RESD 47 DO 30 INODE=1,NNODE RESD 48 LNODE=LNODS(IELEM,INODE) RESD 49 ``` DO 30 IDIME=1,NDIME RESD 50 IPOSN=IPOSN+1 RESD 51 NPOSN=LEQNS(IPOSN, IELEM) RESD 52 IF(NPOSN.EQ.0) DISPT=0. RESD 53 IF(NPOSN.NE.0) DISPT=DISPL(NPOSN) RESD 54 DLCOD(IDIME, INODE)=COORD(LNODE, IDIME)+DISPT RESD 55 ELCOD(IDIME, INODE)=COORD(LNODE, IDIME) RESD 56 30 ELDIS(IDIME, INODE)=DISPT RESD 57 CALL MODPS (DMATX, LPROP, NMATS, NSTRE, NTYPE, PROPS) RESD 58 KGASP=0 RESD 59 DO 40 IGAUS=1,NGAUS RESD 60 DO 40 JGAUS=1,NGAUS RESD 61 EXISP=POSGP(IGAUS) RESD 62 ETASP=POSGP(JGAUS) RESD 63 KGAUS=KGAUS+1 RESD 64 KGASP=KGASP+1 RESD 65 CALL SFR2 (DERIV, NNODE, SHAPE, EXISP, ETASP) RESD 66 CALL JACOB2 (CARTD, DERIV, DJACB, ELCOD, GPCOD, RESD 67 IELEM, KGASP, NNODE, SHAPE) RESD 68 CALL JACOBD (CARTD, DLCOD, DJACM, NDIME, NLAPS, NNODE) RESD 69 DVOLU=DJACB\*WEIGP(IGAUS)\*WEIGP(JGAUS) RESD 70 IF(NTYPE.EQ.3) DVOLU=DVOLU\*TWOPI\*GPCOD(1,KGASP) RESD 71 IF(NTYPE.EQ.1) DVOLU=DVOLU\*THICK RESD 72 CALL BLARGE (BMATX, CARTD, DJACM, DLCOD, GPCOD, RESD 73 KGASP, NLAPS, NNODE, NTYPE, SHAPE) RESD 74 CALL LINGNL (CARTD, DJACM, DMATX, ELDIS, GPCOD, KGASP, RESD 77) KGAUS, NDOFN, NLAPS, NNODE, NSTRE, NTYPE, RESD 78 POISS, SHAPE, STRAN, STRES, STRAG) RESD 79 DO 560 ISTR1=1,NSTR1 RESD 80 560 STRAG(ISTR1,KGAUS)=STRAG(ISTR1,KGAUS)+STRAN(ISTR1) RESD 81 IF(ISTEP.GT.1.AND.IITER.GT.1) GO TO 160 RESD 82 DO 170 ISTR1=1,NSTR1 RESD 83 170 STRES(ISTR1)=STRES(ISTR1)+STRIN(ISTR1,KGAUS) RESD 84 160 CONTINUE RESD 85 PREYS=UNIAX+EPSTN(KGAUS)\*HARDS RESD 86 DO 150 ISTR1=1,NSTR1 RESD 87 DESIG(ISTR1)=STRES(ISTR1) RESD 88 150 SIGMA(ISTR1)=STRSG(ISTR1,KGAUS)+STRES(ISTR1) RESD 89 IF(NLAPS.EQ.2.OR.NLAPS.EQ.0) GO TO 60 RESD 90 CALL INVAR (DEVIA, LPROP, NCRIT, NMATS, PROPS, SINT3, STEFF, RESD 91 SIGMA, THETA, VARJ2, YIELD) RESD 92 ESPRE=EFFST(KGAUS)-PREYS RESD 93 IF(ESPRE.GE.0.0) GO TO 50 RESD 94 ESCUR=YIELD-PREYS RESD 95 IF(ESCUR.LE.0.0) GO TO 60 RESD 96 RFACT=ESCUR/(YIELD-EFFST(KGAUS)) RESD 97 GO TO 70 RESD 98 50 ESCUR=YIELD-EFFST(KGAUS) RESD 99 IF(ESCUR.LE.0.0) GO TO 60 RESD 100 RFACT=1.0 RESD 101 70 MSTEP=ESCUR\*8.0/UNIAX+1.0 RESD 102 IF(MSTEP.GT.10) MSTEP=10 RESD 103 ASTEP=MSTEP RESD 104 REDUC=1.0-RFACT RESD 105 DO 80 ISTR1=1,NSTR1 RESD 106 SGTOT(ISTR1)=STRSG(ISTR1,KGAUS)+REDUC\*STRES(ISTR1) RESD 107 80 STRES(ISTR1)=RFACT\*STRES(ISTR1)/ASTEP RESD 108 DO 90 JSTEP=1,MSTEP RESD 109 CALL INVAR (DEVIA, LPROP, NCRIT, NMATS, PROPS, SINT3, STEFF, RESD 110) SGTOT, THETA, VARJ2, YIELD) RESD 111 CALL YIELDF (AVECT,DEVIA,FRICT,NCRIT,SINT3,STEFF, RESD 112 . THETA,VARJ2) RESD 113 CALL FLOWPL (AVECT, ABETA, DVECT, HARDS, NTYPE, POISS, YOUNG) RESD 114 ```csv AGASH=0.0 RESD 115 DO 100 ISTR1=1,NSTR1 RESD 116 100 AGASH=AGASH+AVECT(ISTR1)*STRES(ISTR1) RESD 117 DLAMD=AGASH*ABETA RESD 118 IF(DLAMD.LT.0.0) DLAMD=0.0 RESD 119 BGASH=0.0 RESD 120 DO 110 ISTR1=1,NSTR1 RESD 121 BGASH=BGASH+AVECT(ISTR1)*SGTOT(ISTR1) RESD 122 110 SGTOT(ISTR1)=SGTOT(ISTR1)+STRES(ISTR1)-DLAMD*DVECT(ISTR1) RESD 123 EPSTN(KGAUS)=EPSTN(KGAUS)+DLAMD*BGASH/YIELD RESD 124 90 CONTINUE RESD 125 CALL INVAR (DEVIA,LPROP,NCRIT,NMATS,PROPS,SINT3,STEFF, RESD 126 .SGTOT,THETA,VARJ2,YIELD) RESD 127 CURYS=UNIAX+EPSTN(KGAUS)*HARDS RESD 128 BRING=1.0 RESD 129 IF(YIELD.GT.CURYS) BRING=CURYS/YIELD RESD 130 DO 130 ISTR1=1,NSTR1 RESD 131 130 STRSG(ISTR1,KGAUS)=BRING*SGTOT(ISTR1) RESD 132 EFFST(KGAUS)=BRING*YIELD RESD 133 C*** ALTERNATIVE LOCATION OF STRESS REDUCTION LOOP TERMINATION CARD RESD 134 C 90 CONTINUE RESD 135 C*** GO TO 190 RESD 136 60 DO 180 ISTR1=1,NSTR1 RESD 137 180 STRSG(ISTR1,KGAUS)=STRSG(ISTR1,KGAUS)+DESIG(ISTR1) RESD 138 EFFST(KGAUS)=YIELD RESD 139 C RESD 140 C*** CALCULATE THE EQUIVALENT NODAL FORCES AND ASSOCIATE WITH THE RESD 141 C ELEMENT NODES RESD 142 190 MGASH=0 RESD 143 DO 140 INODE=1,NNODE RESD 144 DO 140 IDQFN=1,NDOFN RESD 145 MGASH=MGASH+1 RESD 146 DO 140 ISTRE=1,NSTRE RESD 147 140 ELOAD(IELEM,MGASH)=ELOAD(IELEM,MGASH)+BMATX(ISTRE,MGASH)* RESD 148 .STRSG(ISTRE,KGAUS)*DVOLU RESD 149 40 CONTINUE RESD 150 20 CONTINUE RESD 151 DO 500 IELEM=1,NELEM RESD 152 DO 500 IEVAB=1,NEVAB RESD 153 LMVEB=LEQNS(IEVAB,IELEM) RESD 154 IF(LMVEB.EQ.0) GO TO 550 RESD 155 RESID(LMVEB)=RESID(LMVEB)+ELOAD(IELEM,IEVAB) RESD 156 550 CONTINUE RESD 157 500 CONTINUE RESD 158 RETURN RESD 159 END RESD 160 RESD 161 ``` # 11.6 Examples # 11.6.1 Spherical shell example Some of the capabilities $^{(10)}$ of the program MIXDYN are explained by analysing some simple problems. The spherical shell problem described $^{(11,12)}$ in Chapter 10 is again solved for the following cases: (i) Elastic small deformation (all implicit elements) (ii) Elastic geometrically nonlinear (all implicit elements) (iii) Elasto-plastic small deformation (all implicit elements) (iv) Elastic small deformation (all explicit elements) (v) Elastic geometrically nonlinear (all explicit elements) (vi) Elasto-plastic small deformation (all explicit elements) ![](images/page-469_bc51ca27ea45df3b15bf84651b88a05ecaa35ae89a43f3eff4a464b7a4fae73e.jpg)
text_image p = 600 lb/in² Two stiff elements R = 22.27 in α
Fig. 11.4 Modified spherical shell example with stiff elements. To demonstrate the capabilities of program MIXDYN we also solve a slightly modified version of the spherical shell example. Two stiff and dense elements are added to the finite element mesh at the crown as shown in Fig. 11.4. The stiff elements have the following properties: $$ \text { Elastic modulus } \quad E = 0. 1 0 5 \times 1 0 ^ {9} \text { lb / in } ^ {2} $$ $$ \text { poisson's ratio } \quad \nu = 0. 3 $$ $$ \text { mass density } \quad \rho = 0. 7 8 0 \times 1 0 ^ {- 3} \text { lb.sec } ^ {2} / \text { in } ^ {4} $$ $$ \text { yield stress } \quad \sigma_ {0} = 0. 5 \times 1 0 ^ {5} \mathrm{lb/in} ^ {2} $$ The following modified shell examples are also analysed: (vii) Elasto-plastic small deformations (all implicit elements) (viii) Elasto-plastic small deformations (all explicit elements) (ix) Elasto-plastic small deformations (stiff elements are implicit elements, the remaining elements are explicit). The highest and lowest eigenvalues are evaluated for both the original and the modified spherical shells. For the original spherical shell the fundamental period is $0.547 \times 10^{-3}$ sec and the smallest time period is $1.380 \times 10^{-6}$ sec. For the modified spherical shell the fundamental period $T_{f}$ is $0.592 \times 10^{-3}$ sec and the smallest time period $T_{h}$ is $0.776 \times 10^{-6}$ sec. Thus the addition of the stiff elements does not significantly change the largest period but it does change the smallest period quite dramatically. For an accurate solution based on implicit time integration the time step length $\Delta t$ is taken as $T_{f}/100 \simeq 0.6 \times 10^{-5}$ sec for both the original and the modified spherical shell. For a stable and accurate solution based on explicit time integration the time step length $\Delta t \leqslant T_{h}/\pi$ which is $0.25 \times 10^{-6}$ sec for the modified spherical shell or $0.40 \times 10^{-6}$ sec for the original spherical shell. Thus the addition of two stiff elements reduces the critical time step length to 1/1.6 of the original critical time step length. Hence the explicit analysis becomes more expensive. However, if the stiff elements are taken as implicit elements in case (ix) for implicit-explicit analysis, then the critical time step is governed by the remaining explicit elements so that the time step must be less than or equal to $0.40 \times 10^{-6}$ sec. ![](images/page-470_1223531b36a8f36460cb4cf02a18883696b6a76784ae154818e3217a08179ebc.jpg)
line | t (x10⁻³ secs) | Linear small displacement formulation | Linear large displacement formulation | | -------------- | ------------------------------------- | ------------------------------------- | | 0 | 0.00 | 0.00 | | 2 | 0.03 | 0.035 | | 4 | 0.015 | 0.018 | | 6 | 0.035 | 0.04 | | 8 | 0.01 | 0.015 | | 10 | 0.00 | 0.005 |
Fig. 11.5(a) Spherical shell results. Cases (i), (ii), (iv) and (v). Figure 11.5(a) compares the response of the elastic analyses with small and large deformations.\* The results are similar to the results obtained using DYNPAK. The response with the large deformation gives a time period which is elongated. \* Note that the implicit and explicit results overlap.