text_image

F=0 dσe^r=Ddε^r (1-r)dσe^r rdσe^r σ2 σr-1 C B Ddλa= dλd_D σε^r Ddλa=Ddε_p^r σr σ3 σ1 σr=σr(σγ^0+H'ε_p^r)/σ̃^r D' σr

Fig. 7.10(b) Incremental stress changes at a point in an elasto-plastic continuum at initial yield. Step e. For yielded Gauss points only compute the portion of the total stress which satisfies the yield criterion as $\sigma^{r-1} + (1 - R)d\sigma_{e}^{r}$ . Step f. The remaining portion of stress, $Rd\sigma_{e}^{r}$ must be effectively eliminated in some way. The point A must be brought onto the yield surface by allowing plastic deformation to occur. Physically this can be described as follows. On loading from point C, the stress point moves elastically until the yield surface is met at B. Elastic behaviour beyond this point would result in a final stress state defined by point A. However in order to satisfy the yield criterion, the stress point cannot move outside the yield surface and consequently the stress point can only traverse the surface until both equilibrium conditions and the constitutive relation are satisfied. From (7.45), (7.46) and (7.47) we have $$ d \sigma^ {r} = D d \epsilon^ {r} - d \lambda d _ {D}, \tag {7.91} $$ or $$ \sigma^ {r} = \sigma^ {r - 1} + d \sigma_ {e} ^ {r} - d \lambda d _ {D}, \tag {7.92} $$ which gives the total stresses $\sigma^{r}$ satisfying elasto-plastic conditions when the stresses are incremented from $\sigma^{r-1}$ . Expression (7.92) is illustrated vectorially in Fig. 7.10 and the reader should note the similarity to Fig. 3.7(a). It is seen that if a finite sized stress increment is taken, the final stress point D, corresponding to $\sigma^{r}$ , may depart from the yield surface. This discrepancy can be practically eliminated by ensuring that the load increments considered in solution are sufficiently small. However the point D can be reduced to the yield surface by simply scaling the vector $\sigma^{r}$ . Denoting the effective stress, given by Col. 3, Table 7.2, due to stress $\sigma^{r}$ as $\bar{\sigma}^{r}$ and noting that this value should coincide with $\sigma_{Y} = \sigma_{Y}^{\circ} + H' \bar{\epsilon}_{p}^{r}$ if the point D lies on the yield surface, the appropriate scaling factor is readily seen to be $$ \sigma^ {r} = \sigma^ {r} \left(\frac {\sigma_ {Y} ^ {0} + H ^ {\prime} \bar {\epsilon} _ {p} ^ {r}}{\bar {\sigma} ^ {r}}\right). \tag {7.93} $$ This represents a scaling of the vector $\sigma^{r}$ which implies that the individual stress components are proportionally reduced. The normality condition for the plastic strain increment is evident from Fig. 7.10 since $Dd\lambda a = Dd\varepsilon_{p}$ . 

text_image

σ₂ σᵣ₋₁ C σᵣ = σᵣ \left( \frac{\sigma_{\gamma^{0}} + H' \epsilon_{p}}{\bar{\sigma}^r} \right) σ₁ σ₃ A B E' D' D σ₃

Fig. 7.11 Refined process for reducing a stress point to the yield surface. If relatively large load increment sizes are to be permitted the process described above can lead to an inaccurate prediction of the final point D on the yield surface if the stress point is in the vicinity of a region of large curvature of the yield surface. This is illustrated in Fig. 7.11 where the process of reducing the elastic stress to the yield surface is shown to end in the stress point D which is then scaled down to the yield surface to give point $D'$ . Greater accuracy can be achieved by relaxing the excess stress to the yield surface in several stages.\* Fig. 7.11 shows the case where the excess stress is divided into three equal parts and each increment reduced to the yield surface in turn. After the three reduction cycles to the stress point E the drift away from the yield surface can be corrected by simple scaling to give the final stress point $E'$ . It is seen that the final • Alternative procedures for this operation are presented in Refs. 18 and 19 whilst a completely different approach to stress projection is followed in Ref. 20. points $D'$ and $E'$ can be significantly different. An additional refinement which can be introduced is to scale the stress point to the yield surface after the reduction process for each cycle and not only after the final cycle as shown in Fig. 7.11. Obviously the greater the number of steps into which the excess stress AB is divided, the greater the accuracy. However the computation for each step is relatively expensive since the vectors a and $d_{D}$ have to be calculated at each stage. Clearly a balance must be sought and in this text the following criterion is adopted. The excess stress $Rd\sigma_{e}^{r}$ is divided into m parts where m is given by the nearest integer which is less than $$ \left(\frac {\bar {\sigma} _ {e} ^ {r} - \sigma_ {Y}}{\sigma_ {Y} \cdot^ {0}}\right) 8 \div 1, \tag {7.94} $$ where $\bar{\sigma}_{e}^{r}-\sigma_{Y}$ gives a measure of the excess stress AB and $\sigma_{Y}^{\circ}$ is the initial uniaxial yield stress in Col. 4, Table 7.2 before the onset of work hardening. This criterion can be readily amended by the user. Step g. For elastic Gauss points only calculate $\sigma^r$ as $\sigma^r = \sigma^{r-1} + d\sigma_e^r$ . Step h. Finally, calculate the equivalent nodal forces from the element stresses according to $$ (f ^ {(e)}) ^ {r} = \int_ {\Omega^ {(e)}} \boldsymbol {B} ^ {T} \boldsymbol {\sigma} ^ {r} d \Omega . \tag {7.95} $$ Subroutine RESIDU is now listed and described. SUBROUTINE RESIDU(ASDIS,COORD,EFFST,ELOAD,FACTO,IITER,LNODS, RSDU 1 . LPROP,MATNO,MELEM,MMATS,MPOIN,MTOTG,MTOTV,NDOFN, RSDU 2 . NELEM,NEVAB,NGAUS,NNODE,NSTR1,NTYPE,POSGP,PROPS, RSDU 3 . NSTRE,NCrit,STRSG,WEIGP,TDISP,EPSTN) RSDU 4 C********** RSDU 5 C RSDU 6 C*** THIS SUBROUTINE REDUCES THE STRESSES TO THE YIELD SURFACE AND RSDU 7 C EVALUATES THE EQUIVALENT NODAL FORCES RSDU 8 C RSDU 9 C********** RSDU 10 DIMENSION ASDIS(MTOTV),AVECT(4),CARTD(2,9),COORD(MPOIN,2), RSDU 11 . DEVIA(4),DVECT(4),EFFST(MTOTG),ELCOD(2,9),ELDIS(2,9), RSDU 12 . ELOAD(MELEM,18),LNODS(MELEM,9),POSGP(4),PROPS(MMATS,7), RSDU 13 . STRAN(4),STRES(4),STRSG(4,MTOTG), RSDU 14 . WEIGP(4),DLCOD(2,9),DESIG(4),SIGMA(4),SGTOT(4), RSDU 15 . DMATX(4,4),DERIV(2,9),SHAPE(9),GPCOD(2,9), RSDU 16 . EPSTN(MTOTG),TDISP(MTOTV),MATNO(MELEM),BMATX(4,18) RSDU 17 ROOT3=1.73205080757 RSDU 18 TWOPI=6.283185308 RSDU 19 DO 10 IELEM=1,NELEM RSDU 20 DO 10 IEVAB=1,NEVAB RSDU 21 10 ELOAD(IELEM,IEVAB)=0.0 RSDU 22 KGAUS=0 RSDU 23 DO 20 IELEM=1,NELEM RSDU 24 LPROP=MATNO(IELEM) RSDU 25 UNIAX=PROPS(LPROP,5) RSDU 26 HARDS=PROPS(LPROP,6) RSDU 27 ```csv FRICT=PROPS(LPROP,7) RSDU 28 IF(NCRIT.EQ.3)UNIAX=PROPS(LPROP,5)*COS(FRICT*0.017453292) RSDU 29 IF(NCRIT.EQ.4)UNIAX=6.0*PROPS(LPROP,5)*COS(FRICT*0.017453292)/ RSDU 30 .(ROOT3*(3.0-SIN(FRICT*0.017453292))) RSDU 31 C RSDU 32 C*** COMPUTE COORDINATE. AND INCREMENTAL DISPLACEMENTS OF THE RSDU 33 C ELEMENT NODAL POINTS RSDU 34 C RSDU 35 DO 30 INODE =1,NNODE RSDU 36 LNODE=IABS(LNODS(IELEM,INODE)) RSDU 37 NPOSN=(LNODE-1)*NDOFN RSDU 38 DO 30 IDOFN=1,NDOFN RSDU 39 NPOSN=NPOSN+1 RSDU 40 ELCOD(IDOFN,INODE)=COORD(LNODE,IDOFN) RSDU 41 30 ELDIS(IDOFN,INODE)=ASDIS(NPOSN) RSDU 42 CALL MODPS(DMATX,LPROP,MMATS,NTYPE,PROPS) RSDU 43 THICK=PROPS(LPROP,3) RSDU 44 KGASP=0 RSDU 45 DO 40 IGAUS=1,NGAUS RSDU 46 DO 40 JGAUS=1,NGAUS RSDU 47 EXISP=POSGP(IGAUS) RSDU 48 ETASP=POSGP(JGAUS) RSDU 49 KGAUS=KGAUS+1 RSDU 50 KGASP=KGASP+1 RSDU 51 CALL SFR2(DERIV,ETASP,EXISP,NNODE,SHAPE) RSDU 52 CALL JACOB2(CARTD,DERIV,DJACB,ELCOD,GPCOD,IELEM,KGASP, RSDU 53 NNODE,SHAPE) RSDU 54 DVOLU=DJACB*WEIGP(IGAUS)*WEIGP(JGAUS) RSDU 55 IF(NTYPE.EQ.3)DVOLU=DVOLU*TWOPI*GPCOD(1,KGASP) RSDU 56 IF(THICK.NE.0.0)DVOLU=DVOLU*THICK RSDU 57 CALL BMATPS(BMATX,CARTD,NNODE,SHAPE,GPCOD,NTYPE,KGASP) RSDU 58 CALL LINEAR(CARTD,DMATX,ELDIS,LPROP,MMATS,NDOFN,NNODE,NSTRE, RSDU 59 NTYPE,PROPS,STRAN,STRES,KGASP,GPCOD,SHAPE) RSDU 60 PREYS=UNIAX+EPSTN(KGAUS)*HARDS RSDU 61 DO 150 ISTR1=1,NSTR1 RSDU 62 DESIG(ISTR1)=STRES(ISTR1) RSDU 63 150 SIGMA(ISTR1)=STRSG(ISTR1,KGAUS)+STRES(ISTR1) RSDU 64 CALL INVAR(DEVIA,LPROP,MMATS,NCRIT,PROPS,SINT3,STEFF,SIGMA, RSDU 65 THETA,VARJ2,YIELD) RSDU 66 ESPRE=EFFST(KGAUS)-PREYS RSDU 67 IF(ESPRE.GE.0.0)GO TO 50 RSDU 68 ESCUR=YIELD-PREYS RSDU 69 IF(ESCUR.LE.0.0)GO TO 60 RSDU 70 RFACT=ESCUR/(YIELD-EFFST(KGAUS)) RSDU 71 GO TO 70 RSDU 72 50 ESCUR=YIELD-EFFST(KGAUS) RSDU 73 IF(ESCUR.LE.0.0)GO TO 60 RSDU 74 RFACT=1.0 RSDU 75 70 MSTEP=ESCUR*8.0/UNIAX+1.0 RSDU 76 ASTEP=MSTEP RSDU 77 REDUC=1.0-RFACT RSDU 78 DO 80 ISTR1=1,NSTR1 RSDU 79 SGTOT(ISTR1)=STRSG(ISTR1,KGAUS)+REDUC*STRES(ISTR1) RSDU 80 80 STRES(ISTR1)=RFACT*STRES(ISTR1)/ASTEP RSDU 81 DO 90 ISTEP=1,MSTEP RSDU 82 CALL INVAR(DEVIA,LPROP,MMATS,NCRIT,PROPS,SINT3,STEFF,SGTOT, RSDU 83 THETA,VARJ2,YIELD) RSDU 84 CALL YIELDF(AVECT,DEVIA,LPROP,MMATS,NCRIT,NSTR1, RSDU 85 PROPS,SINT3,STEFF,THETA,VARJ2) RSDU 86 CALL FLOWPL(AVECT,ABETA,DVECT,NTYPE,PROPS,LPROP,NSTR1,MMATS) RSDU 87 AGASH=0.0 RSDU 88 DO 100 ISTR1=1,NSTR1 RSDU 89 100 AGASH=AGASH+AVECT(ISTR1)*STRES(ISTR1) RSDU 90 DLAMD=AGASH*ABETA RSDU 91 ```

IF(DLAMD.LT.0.0) DLAMD=0.0	RSDU	92
BGASH=0.0	RSDU	93
DO 110 ISTR1=1,NSTR1	RSDU	94
BGASH=BGASH+AVECT(ISTR1)*SGTOT(ISTR1)	RSDU	95
110 SGTOT(ISTR1)=SGTOT(ISTR1)+STRES(ISTR1)-DLAMD*DVECT(ISTR1)	RSDU	96
EPSTN(KGAUS)=EPSTN(KGAUS)+DLAMD*BGASH/YIELD	RSDU	97
90 CONTINUE	RSDU	98
CALL INVAR(DEVIA,LPROP,MMATS,NCRIT,PROPS,SINT3,STEFF,SGTOT,THETA,VARJ2,YIELD)	RSDU	99
CURYS=UNIAX+EPSTN(KGAUS)*HARDS	RSDU	100
BRING=1.0	RSDU	101
IF(YIELD.GT.CURYS) BRING=CURYS/YIELD	RSDU	102
DO 130 ISTR1=1,NSTR1	RSDU	103
130 STRSG(ISTR1,KGAUS)=BRING*SGTOT(ISTR1)	RSDU	104
EFFST(KGAUS)=BRING*YIELD	RSDU	105
C*** ALTERNATIVE LOCATION OF STRESS REDUCTION LOOP TERMINATION CARD	RSDU	106
C 90 CONTINUE	RSDU	107
C***	RSDU	108
GO TO 190	RSDU	109
60 DO 180 ISTR1=1,NSTR1	RSDU	110
180 STRSG(ISTR1,KGAUS)=STRSG(ISTR1,KGAUS)+DESIG(ISTR1)	RSDU	111
EFFST(KGAUS)=YIELD	RSDU	112
C	RSDU	113
C*** CALCULATE THE EQUIVALENT NODAL FORCES AND ASSOCIATE WITH THE	RSDU	114
C ELEMENT NODES	RSDU	115
190 MGASH=0	RSDU	116
DO 140 INODE=1,NNODE	RSDU	117
DO 140 IDOFN=1,NDOFN	RSDU	118
MGASH=MGASH+1	RSDU	119
DO 140 ISTRE=1,NSTRE	RSDU	120
140 ELOAD(IELEM,MGASH)=ELOAD(IELEM,MGASH)+BMATX(ISTRE,MGASH)*	RSDU	121
.STRSG(ISTRE,KGAUS)*DVOLU	RSDU	122
40 CONTINUE	RSDU	123
20 CONTINUE	RSDU	124
RETURN	RSDU	125
END	RSDU	126
	RSDU	127

RSDU 18–19 Compute $\sqrt{(3)}$ and $2\pi$ . RSDU 20-22 Zero the array in which the equivalent nodal forces, calculated in Step h, will be stored. RSDU 23 Zero the Gauss point counter over all elements. RSDU 24 Loop over each element. RSDU 25 Identify the element material property number. RSDU 26-28 Identify the initial uniaxial yield stress, $\sigma_{Y}^{\circ}$ (or $c$ for Mohr-Coulomb or Drucker-Prager criteria), the linear strain hardening parameter $H'$ and the friction angle $\phi$ for Mohr-Coulomb and Drucker-Prager materials. RSDU 29 For a Mohr-Coulomb material evaluate the equivalent yield stress as $c \cos \phi$ . RSDU 30-31 For a Drucker–Prager material evaluate the equivalent yield stress as $k'$ according to (7.18). RSDU 36-42 Store the element nodal coordinates in array ELCOD and the nodal displacements due to the application of the residual forces in array ELDIS. RSDU 43 RSDU 44 RSDU 45 RSDU 46–49 RSDU 50-51 RSDU 52 RSDU 53-54 RSDU 55-57 RSDU 58 RSDU 59–60 RSDU 61 RSDU 62-64 RSDU 65-66 RSDU 67-68 RSDU 69-70 RSDU 71 Evaluate the elastic D matrix. Identify the element thickness. Zero the local Gauss point counter. Enter the loops for numerical integration and evaluate the local coordinates $(\xi, \eta)$ at the sampling point. Increment the local and global Gauss point counters. Evaluate the shape functions $N_{i}$ and their derivatives $\partial N_{i}/\partial\xi$ , $\partial N_{i}/\partial\eta$ . Evaluate the Gauss point coordinates GPCOD(IDIME, KGASP), the determinant of the Jacobian matrix $|J|$ and the Cartesian derivatives of the shape functions $\partial N_{i} / \partial_{x}$ , $\partial N_{i} / \partial y$ (or $\partial N_{i} / \partial r$ , $\partial N_{i} / \partial z$ for axisymmetric problems). Calculate the elemental volume for numerical integration as $|J|W_{\xi}W_{\eta}$ taking care to multiply by the appropriate thickness or by $2\pi r$ for axisymmetric problems. The default value of the thickness is 1.0. Compute the strain matrix $\pmb{B}$ for the Gauss point. Compute the stress increment STRES(ISTR1), assuming elastic behaviour as $d\sigma_{e}^{r} = Dd\epsilon^{r}$ . Compute the yield stress for the $(r - 1)^{\mathrm{th}}$ iteration as $\sigma_{Y}^{\circ} + H^{\prime}\bar{\epsilon}_{p}^{r - 1}$ . Store $d\sigma_{e}^{r}$ as DESIG(ISTR1) and $\sigma_{e}^{r}$ as SIGMA(ISTR1). Evaluate the effective stress in Col. 3, Table 7.2 and store as YIELD. Check if the Gauss point had yielded on the previous iteration, i.e. if $\bar{\sigma}^{r-1} > \sigma_Y^\circ + H' \bar{\epsilon}_p^{r-1}$ which is the first operation of Step $d$ . If the Gauss point was previously elastic, check to see if it has yielded during this iteration. For a Gauss point which yields during the iteration calculate $$ R = \frac {\bar {\sigma} _ {e} ^ {r} - \sigma_ {Y}}{\bar {\sigma} _ {e} ^ {r} - \bar {\sigma} ^ {r - 1}}. $$ RSDU 73-74 RSDU 75 RSDU 76–77 RSDU 78 RSDU 79-81 Check to see if a Gauss point which had previously yielded is unloading during this iteration. If yes, go to 60. Otherwise, set $R = 1$ . Evaluate the number of steps into which the excess stress, $Rd\sigma_{e}^{r}$ is to be divided according to (7.94). Compute $(1 - R)$ . Compute $\sigma^{r-1}+(1-R)d\sigma_{e}^{r}$ according to Step $e$ and store in SGTOT(ISTR1) and evaluate $Rd\sigma_{e}^{r}/m$ and store in STRES(ISTR1). RSDU 82 RSDU 83-87 Loop over each stress reduction step. Compute the vectors $a$ and $d_{D}$ . RSDU 88-92 Compute $d\lambda$ according to (7.45) and store as DLAMD. RSDU 93-96 Compute $\sigma^r = \sigma^{r-1} + (1 - R)d\sigma_e^r + Rd\sigma_e^r/m - d\lambda d_D/m$ . When the summation process from 1 to $m$ required in DO LOOP to index 90 is completed this will result in $\sigma^r = \sigma^{r-1} + d\sigma_e^r - d\lambda d_D$ to give the stress point $E$ in Fig. 7.11. RSDU 97 Compute the effective plastic strain as follows. From (7.51) we have $$ d \kappa = d \lambda a ^ {T} \sigma = \sigma^ {T} d \epsilon_ {p}, $$ or rewriting the right hand side in terms of the effective stress $\bar{\sigma}$ and effective plastic strain $\bar{\epsilon}_{p}$ we have $$ d \lambda a ^ {T} \sigma = \bar {\sigma} d \bar {\epsilon} _ {p}, $$ and therefore $$ \bar {\epsilon} _ {p} ^ {r} = \bar {\epsilon} _ {p} ^ {r - 1} + \frac {d \lambda \boldsymbol {a} ^ {T} \boldsymbol {\sigma}}{\bar {\sigma}}. \tag {7.96} $$ RSDU 98 Return to loop over the next stress reduction step. This statement is so placed that the final stresses $\sigma^{r}$ are scaled down to lie on the yield surface only after all the reduction steps have been completed. An additional refinement can be introduced where, with reference to Fig. 7.11, the stresses are scaled to the yield surface after each reduction step. Such a refinement is not normally required; however it can be introduced by moving statement RSDU 98 to the position indicated in RSDU 108. RSDU 99-100 Compute the effective stress $\bar{\sigma}^r$ . RSDU 101 Evaluate $\sigma_{Y}^{\circ} + H^{\prime}\bar{\epsilon}_{p}^{r}$ . RSDU 102–105 Factor the stresses $\sigma^{r}$ to ensure that they lie on the yield surface, according to $\sigma^{r} = \sigma^{r}(\sigma_{Y}^{\circ} + H' \bar{\epsilon}_{p}^{r}) / \bar{\sigma}^{r}$ as indicated in Fig. 7.11. RSDU 106 Store the effective stress $\bar{\sigma}^{r}$ in array EFFST. RSDU 108 Location of end of loop if the refinement indicated in RSDU 98 is to be included. RSDU 111-113 For elastic Gauss points compute $\sigma^r$ as $\sigma^{r-1} + d\sigma_e^r$ and store $\bar{\sigma}^r$ in EFFST. RSDU 117-123 Compute the equivalent nodal forces as $$ (f ^ {(e)}) ^ {r} = \int_ {\Omega} \boldsymbol {B} ^ {T} \sigma^ {r} d \Omega . $$ RSDU 124–125 Termination of loop for numerical integration and over each element respectively. # 7.8.8 Subroutine OUTPUT This subroutine outputs the results at a frequency determined by the output parameters NOUTP(1) and NOUTP(2) whose role is described in Section 6.5.3. The principal stresses and direction are also calculated in this subroutine and these are given by the following expressions $$ \sigma_ {\max} = \frac {\sigma_ {x} + \sigma_ {y}}{2} + \sqrt {\left(\frac {(\sigma_ {x} - \sigma_ {y}) ^ {2}}{4} + \tau_ {x y} ^ {2}\right)}, $$ $$ \sigma_ {\min} = \frac {\sigma_ {x} + \sigma_ {y}}{2} - \sqrt {\left(\frac {(\sigma_ {x} - \sigma_ {y}) ^ {2}}{4} + \tau_ {x y} ^ {2}\right)} $$ $$ \theta = \tan^ {- 1} \left(\frac {2 \tau_ {x y}}{\sigma_ {x} - \sigma_ {y}}\right). \tag {7.97} $$ with x and y being replaced by r and z for the axisymmetric case. The term $\theta$ defines the angle which the maximum principal stress makes with the y (or z) axis; a positive angle being measured anticlockwise. This subroutine is largely self-explanatory and is listed below. ```csv SUBROUTINE OUTPUT(IITER,MTOTG,MTOTV,MVFIX,NELEM,NGAUS,NOFIX, OTPT 1 . NOUTP,NPOIN,NVFIX,STRSG,TDISP,TREAC,EPSTN, OTPT 2 . NTYPE,NCHEK) OTPT 3 C******************************* C C**** THIS SUBROUTINE OUTPUTS DISPLACEMENTS.REACTIONS AND STRESSES OTPT 6 C C******************************* DIMENSION NOFIX(MVFIX),NOUTP(2),STRSG(4,MTOTG),STRSP(3), OTPT 9 . TDISP(MTOTV),TREAC(MVFIX,2),EPSTN(MTOTG) OTPT 10 KOUTP=NOUTP(1) OTPT 11 IF(IITER.GT.1) KOUTP=NOUTP(2) OTPT 12 IF(IITER.EQ.1.AND.NCHEK.EQ.0) KOUTP=NOUTP(2) OTPT 13 C C*** OUTPUT DISPLACEMENTS OTPT 15 C IF(KOUTP.LT.1) GO TO 10 OTPT 16 WRITE(6,900) OTPT 17 900 FORMAT(1HO,5X,13HDISPLACEMENTS) OTPT 18 IF(NTYPE.NE.3) WRITE(6,950) OTPT 19 950 FORMAT(1HO,6X,4HNODE,6X,7HX-DISP.,7X,7HY-DISP.) OTPT 20 IF(NTYPE.EQ.3) WRITE(6,955) OTPT 21 955 FORMAT(1HO,6X,4HNODE,6X,7HR-DISP.,7X,7HZ-DISP.) OTPT 22 DO 20 IPOIN=1,NPOIN OTPT 23 NGASH=IPOIN*2 OTPT 24 NGISH=NGASH-2+1 OTPT 25 20 WRITE(6,910) IPOIN,(TDISP(IGASH),IGASH=NGISH,NGASH) OTPT 26 910 FORMAT(I10,3E14.6) OTPT 27 10 CONTINUE OTPT 28 C C*** OUTPUT REACTIONS OTPT 30 C IF(KOUTP.LT.2) GO TO 30 OTPT 31 WRITE(6,920) OTPT 32 920 FORMAT(1HO,5X,9HREACTIONS) OTPT 33 IF(NTYPE.NE.3) WRITE(6,960) OTPT 34 OTPT 35 OTPT 36 ``` 960 FORMAT(1HO,6X,4HNODE,6X,7HX-REAC.,7X,7HY-REAC.) OTPT 37 IF(NTYPE.EQ.3) WRITE(6,965) OTPT 38 965 FORMAT(1HO,6X,4HNODE,6X,7HR-REAC.,7X,7HZ-REAC.) OTPT 39 DO 40 IVFIX=1,NVFIX OTPT 40 40 WRITE(6,910) NOFIX(IVFIX),(TREAC(IVFIX,IDOFN),IDOFN=1,2) OTPT 41 30 CONTINUE OTPT 42 C OTPT 43 C*** OUTPUT STRESSES OTPT 44 C OTPT 45 IF(KOUTP.LT.3) GO TO 50 OTPT 46 IF(NTYPE.NE.3) WRITE(6,970) OTPT 47 970 FORMAT(1HO,1X,4HG.P.,6X,9HXX-STRESS,5X,9HYY-STRESS,5X,9HXY-STRESS,OTPT 48 .5X,9HZZ-STRESS,6X,8HMAX P.S.,6X,8HMIN P.S.,3X,5HANGLE,3X, OTPT 49 .6HE.P.S.) OTPT 50 IF(NTYPE.EQ.3) WRITE(6,975) OTPT 51 975 FORMAT(1HO,1X,4HG.P.,6X,9HRR-STRESS,5X,9HZZ-STRESS,5X,9HRZ-STRESS,OTPT 52 .5X,9HTT-STRESS,6X,8HMAX P.S.,6X,8HMIN P.S.,3X,5HANGLE,3X, OTPT 53 .6HE.P.S.) OTPT 54 KGAUS=0 OTPT 55 DO 60 IELEM=1,NELEM OTPT 56 KELGS=0 OTPT 57 WRITE(6,930) IELEM OTPT 58 930 FORMAT(1HO,5X,13HELEMENT NO. =,I5) OTPT 59 DO 60 IGAUS=1,NGAUS OTPT 60 DO 60 JGAUS=1,NGAUS OTPT 61 KGAUS=KGAUS+1 OTPT 62 KELGS=KELGS+1 OTPT 63 XGASH=(STRSG(1,KGAUS)+STRSG(2,KGAUS))*0.5 OTPT 64 XGISH=(STRSG(1,KGAUS)-STRSG(2,KGAUS))*0.5 OTPT 65 XGESH=STRSG(3,KGAUS) OTPT 66 XGOSH=SQRT(XGISH*XGISH+XGESH*XGESH) OTPT 67 STRSP(1)=XGASH+XGOSH OTPT 68 STRSP(2)=XGASH-XGOSH OTPT 69 IF(XGISH.EQ.0.0) XGISH=0.1E-20 OTPT 70 STRSP(3)=ATAN(XGESH/XGISH)*28.647889757 OTPT 71 60 WRITE(6,940) KELGS,(STRSG(ISTR1,KGAUS),ISTR1=1,4), OTPT 72 . (STRSP(ISTRE),ISTRE=1,3),EPSTN(KGAUS) OTPT 73 940 FORMAT(I5,2X,6E14.6,F8.3,E14.6) OTPT 74 50 CONTINUE OTPT 75 RETURN OTPT 76 END OTPT 77 OTPT 11–13 Set the output indicator, KOUTP, according to whether or not this is the first iteration of a load increment or not. If it is the first iteration the results will be output according to NOUTP(1) but for a converged solution the results are output according to NOUTP(2). OTPT 17-29 For an output code value of 1 or greater, output the nodal displacements after printing the appropriate headings. OTPT 33-42 For an output code of 2 or greater, output appropriate headings and the reactions at restrained nodal points. OTPT 46 For an output code of 3 output the Gauss point stresses. OTPT 47-54 Write appropriate headings. OTPT 56–59 Loop over each element and write the element number. OTPT 60–61 Loop over each element Gauss point. OTPT 62–71 Evaluate the principal stresses and direction for each Gauss point according to (7.97). OTPT 72–74 Output the Cartesian stress components, the principal stresses and direction and the total effective plastic strain for each Gauss point. This latter quantity gives an immediate indication whether the Gauss point has yielded or not, since it will be zero for all elastic points. # 7.8.9 The main, master or controlling segment This segment controls the calling, in order, of the other subroutines and is similar in structure to the segment described in Section 3.8 for one-dimensional situations. Its other function is to control the iterative process and also the incrementing of the applied loads. The following channel numbers are employed by the program: 5 (card reader), 6 (line printer), 1, 2, 3, 4, 8 (scratch files). This routine is self-explanatory and is presented below without further comment. ```txt MASTER PLAST C********** PROGRAM FOR THE ELASTO-PLASTIC ANALYSIS OF PLANE STRESS, C PLANE STRAIN AND AXISYMMETRIC SOLIDS C********** DIMENSION ASDIS(300),COORD(150,2),ELOAD(40,18),ESTIF(18,18), EQRHS(10),EQUAT(80,10),FIXED(300),GLOAD(80),GSTIF(3240), IFFIX(300),LNODS(40,9),LOCEL(18),MATNO(40), NACVA(80),NAMEV(10),NDEST(18),NDFRO(40),NOFIX(30), NOUTP(2),NPIVO(10), POSGP(4),PRESC(30,2),PROPS(5,7),RLOAD(40,18), STFOR(300),TREAC(30,2),VECRV(80),WEIGP(4), STRSG(4,360),TDISP(300),TLOAD(40,18), TOFOR(300),EPSTN(360),EFFST(360) PLAS 6 PLAS 7 PLAS 8 PLAS 9 PLAS 10 PLAS 11 PLAS 12 PLAS 13 PLAS 14 PLAS 15 C C*** PRESET VARIABLES ASSOCIATED WITH DYNAMIC DIMENSIONING C CALL DIMEN(MBUFA,MELEM,MEVAB,MFRON,MMATS,MPOIN,MSTIF,MTOTG,MTOTV, MVFIX,NDOFN,NPROP,NSTRE) PLAS 16 PLAS 17 PLAS 18 PLAS 19 PLAS 20 C C*** CALL THE SUBROUTINE WHICH READS MOST OF THE PROBLEM DATA C CALL INPUT(COORD,IFFIX,LNODS,MATNO,MELEM,MEVAB,MFRON,MMATS, MPOIN,MTOTV,MVFIX,NALGO, NCRIT,NDFRO,NDOFN,NELEM,NEVAB,NGAUS,NGAU2, NINCS,NMATS,NNODE,NOFIX,NPOIN,NPROP,NSTRE, NSTR1,NTOTG,NTOTV, NTYPE,NVFIX,POSGP,PRES C,PROPS,WEIGP) PLAS 21 PLAS 22 PLAS 23 PLAS 24 PLAS 25 PLAS 26 PLAS 27 PLAS 28 PLAS 29 C C*** CALL THE SUBROUTINE WHICH COMPUTES THE CONSISTENT LOAD VECTORS C FOR EACH ELEMENT AFTER READING THE RELEVANT INPUT DATA PLAS 30 PLAS 31 PLAS 32 PLAS 33 PLAS 34 PLAS 35 PLAS 36 C C*** INITIALISE CERTAIN ARRAYS PLAS 37 PLAS 38 ```