김경종
24 KiB
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
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.
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(XGISHXGISH+XGESHXGESH) 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.
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