$$ \boldsymbol {a} = C _ {1} \boldsymbol {a} _ {1} + C _ {2} \boldsymbol {a} _ {2} + C _ {3} \boldsymbol {a} _ {3}, \tag {7.69} $$ where $$ \boldsymbol {a} _ {1} ^ {T} = \frac {\partial J _ {1}}{\partial \sigma} = \{1, 1, 1, 0, 0, 0 \} $$ $$ \boldsymbol {a} _ {2} ^ {T} = \frac {\partial (J _ {2} ^ {\prime}) ^ {1 / 2}}{\partial \sigma} = \frac {1}{2 (J _ {2} ^ {\prime}) ^ {1 / 2}} \left\{\sigma_ {x} ^ {\prime}, \sigma_ {y} ^ {\prime}, \sigma_ {z} ^ {\prime}, 2 \tau_ {y z}, 2 \tau_ {z x}, 2 \tau_ {x y} \right\} $$ $$ \boldsymbol {a} _ {3} ^ {T} = \frac {\partial J _ {3}}{\partial \sigma} = \left\{\left(\sigma_ {y} ^ {\prime} \sigma_ {z} ^ {\prime} - \tau_ {y z} ^ {2} + \frac {J _ {2} ^ {\prime}}{3}\right), \quad \left(\sigma_ {x} ^ {\prime} \sigma_ {z} ^ {\prime} - \tau_ {x z} ^ {2} + \frac {J _ {2} ^ {\prime}}{3}\right), \right. $$ $$ \left(\sigma_ {x} ^ {\prime} \sigma_ {y} ^ {\prime} - \tau_ {x y} ^ {2} + \frac {J _ {2} ^ {\prime}}{3}\right), \quad 2 \left(\tau_ {x z} \tau_ {x y} - \sigma_ {x} ^ {\prime} \tau_ {y z}\right), $$ $$ \left. 2 \left(\tau_ {x y} \tau_ {y z} - \sigma_ {y} ^ {\prime} \tau_ {x z}\right), \quad 2 \left(\tau_ {y z} \tau_ {x z} - \sigma_ {z} ^ {\prime} \tau_ {x y}\right) \right\}, \tag {7.70} $$ and $$ C _ {1} = \frac {\partial F}{\partial J _ {1}}, \quad C _ {2} = \left(\frac {\partial F}{\partial \left(J _ {2} ^ {\prime}\right) ^ {1 / 2}} - \frac {\tan 3 \theta}{\left(J _ {2} ^ {\prime}\right) ^ {1 / 2}} \frac {\partial F}{\partial \theta}\right), $$ $$ C _ {3} = \frac {- \sqrt {3}}{2 \cos 3 \theta} \frac {1}{(J _ {2} ^ {\prime}) ^ {3 / 2}} \frac {\partial F}{\partial \theta}. \tag {7.71} $$ Only the constants $C_{1}$ , $C_{2}$ and $C_{3}$ are then necessary to define the yield surface. Thus we can achieve a simplicity of programming as only these three constants have to be varied between one yield surface and another. The constants $C_{i}$ are given in Table 7.1 for the four yield criteria considered in Section 7.2.1 and other yield functions can be expressed in the same form with equal ease. Table 7.1 Constants defining the yield surface in a form suitable for numerical analysis.

Yield Criterion	$C_1$	$C_2$	$C_3$
Tresca	0	$2\cos\theta(1+\tan\theta\tan3\theta)$	$\frac{\sqrt{3}}{J_2'}\frac{\sin\theta}{\cos3\theta}$
Von Mises	0	$\sqrt{3}$	0
Mohr–Coulomb	$\frac{1}{3}\sin\phi$	$\cos\theta[(1+\tan\theta\tan3\theta)+\sin\phi(\tan3\theta-\tan\theta)/\sqrt{3}]$	$\frac{(\sqrt{3}\sin\theta+\cos\theta\sin\phi)}{(2J_2'\cos3\theta)}$
Drucker–Prager	$\alpha$	1.0	0

# 7.5 Basic expressions for two dimensional problems For two dimensional problems, the general expressions derived so far in this chapter have to be modified. Primarily the main alteration required is the deletion of the stress (and strain) components which vanish under the conditions of plane stress, plane strain or axial symmetry. We have only four non-zero stress or strain components, namely $$ \sigma^ {T} = \{\sigma_ {x}, \sigma_ {y}, \tau_ {x y}, \sigma_ {z} \}, \quad \sigma_ {z} = 0 \quad \text { for Plane Stress } $$ $$ \{\sigma_ {x}, \sigma_ {y}, \tau_ {x y}, \sigma_ {z} \}, \quad \epsilon_ {z} = 0 \quad \text { Plane Strain } $$ $$ \left\{\sigma_ {r}, \sigma_ {z}, \tau_ {r z}, \sigma_ {\theta} \right\} \quad \text { Axial Symmetry. } \tag {7.72} $$ From Fig. 7.8 it is seen that the z direction is taken as the coordinate independent direction for plane stress and plane strain. It is also found convenient to order the stress components as indicated in (7.72) with the stress in the coordinate independent direction being last.  Fig. 7.8 Two-dimensional applications showing coordinate systems employed. The explicit form of the elasticity matrix D can be written $$ \boldsymbol {D} = \frac {E (1 - \nu)}{(1 + \nu) (1 - 2 \nu)} \left[ \begin{array}{c c c c} 1 & \frac {\nu}{1 - \nu} & 0 & \frac {\nu}{1 - \nu} \\ \frac {\nu}{1 - \nu} & 1 & 0 & \frac {\nu}{1 - \nu} \\ 0 & 0 & \frac {1 - 2 \nu}{2 (1 - \nu)} & 0 \\ \frac {\nu}{1 - \nu} & \frac {\nu}{1 - \nu} & 0 & 1 \end{array} \right] \text { for plane strain and axial symmetry, } $$ $$ \boldsymbol {D} = \frac {E}{1 - \nu^ {2}} \left[ \begin{array}{c c c c c} 1 & \nu & 0 & 0 \\ \nu & 1 & 0 & 0 \\ 0 & 0 & \frac {1 - \nu}{2} & 0 \\ - \frac {2}{0} & 0 & 0 & 1 \end{array} \right] \quad \text { for plane stress. } \tag {7.73} $$ Note that the components corresponding to the coordinate independent direction have been included for the plane stress and strain cases. These terms will be excluded for element stiffness formulation and only the first $3 \times 3$ portion indicated will be employed. By eliminating the appropriate stress terms the expressions developed to date can be readily modified. The flow vector a becomes $$ \boldsymbol {a} ^ {T} = \left\{\frac {\partial F}{\partial \sigma_ {x}}, \frac {\partial F}{\partial \sigma_ {y}}, \frac {\partial F}{\partial \tau_ {x y}}, \frac {\partial F}{\partial \sigma_ {z}} \right\}, \tag {7.74} $$ with x, y and z being replaced by r, z and $\theta$ respectively for the case of axial symmetry. The specific form of the vector, a is still given by (7.69) but in this case we have from (7.70) $$ \boldsymbol {a} _ {1} ^ {T} = \{1, 1, 0, 1 \} $$ $$ \boldsymbol {a} _ {2} ^ {T} = \frac {1}{2 (J _ {2} ^ {\prime}) ^ {1 / 2}} \left\{\sigma_ {x} ^ {\prime}, \sigma_ {y} ^ {\prime}, 2 \tau_ {x y}, \sigma_ {z} ^ {\prime} \right\} $$ $$ \boldsymbol {a} _ {3} ^ {T} = \left\{\left(\sigma_ {y} ^ {\prime} \sigma_ {z} ^ {\prime} + \frac {J _ {2} ^ {\prime}}{3}\right), \left(\sigma_ {x} ^ {\prime} \sigma_ {z} ^ {\prime} + \frac {J _ {2} ^ {\prime}}{3}\right), \right. $$ $$ \left. - 2 \sigma_ {z} ^ {\prime} \tau_ {x y}, \left(\sigma_ {x} ^ {\prime} \sigma_ {y} ^ {\prime} - \tau_ {x y} ^ {2} + \frac {J _ {2} ^ {\prime}}{3}\right) \right\}, \tag {7.75} $$ and the deviatoric stress invariants become, from (7.5) $$ J _ {2} ^ {\prime} = \frac {1}{2} \left(\sigma_ {x} ^ {\prime 2} + \sigma_ {y} ^ {\prime 2} + \sigma_ {z} ^ {\prime 2}\right) + \tau_ {x y} ^ {2} $$ $$ J _ {3} ^ {\prime} = \sigma_ {z} ^ {\prime} (\sigma_ {z} ^ {\prime 2} - J _ {2} ^ {\prime}). \tag {7.76} $$ To complete the prescription of the elasto-plastic matrix $D_{ep}$ given in (7.47) we require $d_{D}$ . Employing the relevant form of D from (7.73) in (7.47) results in, for plane strain and axial symmetry $$ \boldsymbol {d} _ {D} = \left\{ \begin{array}{l} d _ {1} \\ d _ {1} \\ d _ {3} \\ d _ {4} \end{array} \right\} = \left\{ \begin{array}{c} \frac {E}{1 + \nu} a _ {1} + M _ {1} \\ \frac {E}{1 + \nu} a _ {2} + M _ {1} \\ G a _ {3} \\ \frac {E}{1 + \nu} a _ {4} + M _ {1} \end{array} \right\}, \quad M _ {1} = \frac {E \nu \left(a _ {1} + a _ {2} + a _ {4}\right)}{(1 + \nu) (1 - 2 \nu)}, \tag {7.77} $$ where $G = E / 2(1 + \nu)$ is the shear modulus and $a_1 \ldots a_4$ are the components of $\pmb{a}$ . For plane stress we have $$ \boldsymbol {d} _ {D} = \left\{ \begin{array}{c} \frac {E}{1 + \nu} a _ {1} + M _ {2} \\ \frac {E}{1 + \nu} a _ {2} + M _ {2} \\ G a _ {3} \\ \frac {E}{1 + \nu} a _ {4} + M _ {2} \end{array} \right\}, \quad M _ {2} = \frac {E \nu \left(a _ {1} + a _ {2}\right)}{1 - \nu^ {2}}. \tag {7.78} $$ # 7.6 Singular points on the yield surface For many yield surfaces the flow vector a is not uniquely defined for certain stress combinations. For example this arises at the corners of the Tresca and Mohr–Coulomb criteria located by $\theta = \pm30^{\circ}$ and the direction of plastic straining there is indeterminate. Koiter $^{(12)}$ has provided limits within which the incremental plastic strain vector must lie. Numerical difficulties will be encountered as $\theta$ approaches $\pm30^{\circ}$ for the Tresca and Mohr–Coulomb laws since it is seen from Table 7.1 that for these values of $\theta$ both $C_{2}$ and $C_{3}$ become indeterminate. This difficulty can be overcome by returning to the original expressions (7.63) for the Tresca law and (7.65) for the Mohr–Coulomb criterion and rewriting these for the explicit values $\theta = \pm30^{\circ}$ . Thus we have for the Tresca law $$ \sqrt {(3)} \left(J _ {2} ^ {\prime}\right) ^ {\frac {1}{2}} = Y (\kappa) = \sqrt {(3)} k (\kappa), \tag {7.79} $$ and thus from (7.71) we have $$ C _ {1} = 0, \quad C _ {2} = \sqrt {(3)}, \quad C _ {3} = 0 \quad \text {for} \quad \theta = \pm 3 0 ^ {\circ}. \tag {7.80} $$ Physically, since (7.79) is the Von Mises criterion, this is equivalent to stating that the direction of plastic straining at the corners of the Tresca criterion is that given by the Von Mises circle which also passes through the corner (see Fig. 7.2). Similarly for the Mohr-Coulomb criterion we have from (7.65), $$ \frac {1}{3} J _ {1} \sin \phi + \left(J _ {2} ^ {\prime}\right) ^ {1 / 2} \frac {1}{2} \left(\sqrt {3} - \frac {\sin \phi}{\sqrt {3}}\right) - c \cos \phi = 0 \quad \text { for } \quad \theta = + 3 0 ^ {0} $$ $$ \frac {1}{3} J _ {1} \sin \phi + (J _ {2} ^ {\prime}) ^ {1 / 2} \frac {1}{2} \left(\sqrt {3 + \frac {\sin \phi}{\sqrt {3}}}\right) - c \cos \phi = 0 \quad \theta = - 3 0 ^ {0}, \tag {7.81} $$ or from (7.71) we have $$ C _ {1} = \frac {1}{3} \sin \phi , C _ {2} = \frac {1}{2} \left(\sqrt {3} - \frac {\sin \phi}{\sqrt {3}}\right), C _ {3} = 0 \quad \text { for } \quad \theta = + 3 0 ^ {0} $$ $$ C _ {1} = \frac {1}{3} \sin \phi , C _ {2} = \frac {1}{2} \left(\sqrt {3} + \frac {\sin \phi}{\sqrt {3}}\right), C _ {3} = 0 \quad \theta = - 3 0 ^ {0}. \tag {7.82} $$ The practical approach adopted in this text is to use the general expressions for $C_1, C_2, C_3$ given in Table 7.1 for all values of $|\theta| \leqslant 29^\circ$ and to then employ either (7.80) for Tresca or (7.82) for Mohr-Coulomb in the vicinity of the corners. This makes the direction of straining unique, and also satisfies the Koiter requirements. Physically this artifice corresponds to a 'rounding off' of the yield surface corners. # 7.7 Finite element expressions and program structure The basic expressions required for solution can be again obtained by use of the principle of virtual work. Consider the solid, in which the internal stresses $\sigma$ , the distributed loads/unit volume $b$ and external applied forces $f$ form an equilibrating field, to undergo an arbitrary virtual displacement pattern $\delta d^{*}$ which result in compatible strains $\delta \epsilon^{*}$ and internal displacements $\delta u^{*}$ . Then the principle of virtual work requires that $$ \int_ {\Omega} \left(\delta \epsilon^ {* T} \sigma - \delta u ^ {* T} b\right) d \Omega - \delta d ^ {* T} f = 0. \tag {7.83} $$ Then the normal finite element discretising procedure leads to the following expressions for the displacements and strains within any element $$ \delta \boldsymbol {u} ^ {*} = N \delta \boldsymbol {d} ^ {*}, \quad \delta \boldsymbol {\epsilon} ^ {*} = \boldsymbol {B} \delta \boldsymbol {d} ^ {*}, \tag {7.84} $$ where N and B are respectively the usual matrix of shape functions and the elastic strain matrix. Then the element assembly process gives $$ \int_ {\Omega} \delta \boldsymbol {d} ^ {* T} (\boldsymbol {B} ^ {T} \boldsymbol {\sigma} - \boldsymbol {N} ^ {T} \boldsymbol {b}) d \Omega - \delta \boldsymbol {d} ^ {* T} \boldsymbol {f} = 0, \tag {7.85} $$ where the volume integration over the solid is the sum of the individual element contributions. Since this expression must hold true for any arbitrary $\delta d^{*}$ value $$ \int_ {\Omega} \boldsymbol {B} ^ {T} \sigma d \Omega - \boldsymbol {f} - \int_ {\Omega} \boldsymbol {N} ^ {T} \boldsymbol {b} d \Omega = 0. \tag {7.86} $$ For the solution of nonlinear problems as described in Chapter 2, (7.86) will not generally be satisfied at any stage of the computation, and $$ \psi = \int_ {\Omega} \boldsymbol {B} ^ {T} \sigma d \Omega - \left(f + \int_ {\Omega} N ^ {T} \boldsymbol {b} d \Omega\right) \neq 0, \tag {7.87} $$ where $\psi$ is the residual force vector. For an elasto-plastic situation the material stiffness is continually varying, and instantaneously the incremental stress/strain relationship is given by (7.46). For the purpose of evaluating the material tangential stiffness matrix $K_{T}$ at any stage, the incremental form of (7.87) must be employed. Thus within an increment of load we have $$ \Delta \psi = \int_ {\Omega} \boldsymbol {B} ^ {T} \Delta \sigma d \Omega - \left(\Delta f + \int_ {\Omega} \boldsymbol {N} ^ {T} \Delta \boldsymbol {b} d \Omega\right). \tag {7.88} $$ Substituting for $\Delta\sigma$ from (7.46) results in $$ \Delta \psi = K _ {T} d - \left(\Delta f + \int_ {\Omega} N ^ {T} \Delta b d \Omega\right), \tag {7.89} $$ where $$ \boldsymbol {K} _ {T} = \int_ {\Omega} \boldsymbol {B} ^ {T} \boldsymbol {D} _ {e p} \boldsymbol {B} d \Omega . \tag {7.90} $$ Expression (7.89) is essentially identical to (2.4) and therefore the solution procedures developed in Chapter 2 can be again employed. The programming philosophy adopted for this application follows that employed in Chapter 3 for one-dimensional elasto-plastic problems. It is suggested that the reader reviews the appropriate sections of Chapter 3 before proceeding to the remainder of this chapter. The solution techniques discussed in Chapters 2 and 3 are utilised and in particular an initial stiffness algorithm, a tangential stiffness algorithm and two options of the combined initial/tangential stiffness approach are included. An outline of the program is provided in Fig. 7.9. Many of the subroutines required are common to the corresponding linear elastic solution program and their function and structure have already been described. In particular, subroutines BMATS, CHECK1, CHECK2, DBE, ECHO, FRONT, GAUSSQ, JACOB2, LOADPS, MODPS, NODEXY and SFR2 have been described in Section 6.4. Also the standard nonlinear subroutines ALGOR, CONVER, INCREM and INPUT have been presented in Section 6.5. We will now formulate the additional subroutines required and assemble them to form a working program. 

flowchart

```mermaid graph TD A["START"] --> B["DIMEN
Presents the variables associated with the dynamic dimensioning process."] B --> C["INPUT
Inputs data defining geometry, boundary conditions and material properties."] C --> D["LOADPS
Evaluates the equivalent nodal forces for pressure loading, gravity loading, etc."] D --> E["ZERO
Sets to zero arrays required for accumulation of data."] E --> F["INCREM
Increments the applied loads according to specified load factors."] F --> G["ALGOR
Sets indicator to identify the type of solution algorithm e.g. initial stiffness, tangential stiffness, etc."] G --> H["STIFFP
Calculates the element stiffnesses for elastic and elasto-plastic material behaviour."] H --> I["FRONT
Solves the simultaneous equation system by the frontal method."] I --> J["RESIDU
Calculates the residual force vector, ψ."] J --> K["CONVER
Checks to see if the solution process has converged."] K -->|NO| L["END"] K -->|YES| M["OUTPUT
Prints the results for this load increment."] L --> N["ITERATION LOOP"] M --> N ```

Fig. 7.9 Program organisation for two-dimensional elasto-plastic applications. # 7.8 Additional program subroutines A total of eight additional subroutines are required some of which will be common to other nonlinear applications considered in later chapters of this text. # 7.8.1 Subroutine DIMEN The function of this subroutine is to preset the values of variables employed in the program. In particular the variables associated with the dynamic dimensioning process described in Chapter 6 are defined. Thus if it is required to upgrade the magnitude of the maximum problem size which can be solved it is only necessary to modify the dimension statements in the main or master subroutine together with the variables set in subroutine DIMEN. All the variables preset in this subroutine have been previously defined and their specified values are indicated in the following listing. ```txt SUBROUTINE DIMEN(MBUFA,MELEM,MEVAB,MFRON,MMATS,MPOIN,MSTIF,MTOTG,MTOTV,MVFIX,NDOFN,NPROP,NSTRE) DIMN 1 C*************** DIMN 2 C C**** THIS SUBROUTINE PRESETS VARIABLES ASSOCIATED WITH DYNAMIC DIMN 3 C DIMENSIONING DIMN 4 C C*************** DIMN 5 MBUFA = 10 DIMN 6 MELEM=40 DIMN 7 MFRON=80 DIMN 8 MMATS = 5 DIMN 9 MPOIN=150 DIMN 10 MSTIF=(MFRON*MFRON-MFRON)/2.0+MFRON DIMN 11 MTOTG = MELEM*9 DIMN 12 NDOFN = 2 DIMN 13 MTOTV = MPOIN*NDOFN DIMN 14 MVFIX=30 DIMN 15 NPROP=7 DIMN 16 MEVAB = NDOFN*9 DIMN 17 RETURN DIMN 18 END DIMN 19 DIMN 20 DIMN 21 DIMN 22 ``` # 7.8.2 Subroutine ZERO This subroutine merely sets to zero the contents of several arrays employed in the program. These arrays will be employed to accumulate data as the incremental and iterative process continues and they therefore require to be initialised to zero. This subroutine is self-explanatory and is presented without further comment. ```txt SUBROUTINE ZERO(ELOAD,MELEM,MEVAB,MPOIN,MTOTG,MTOTV,NDOFN,NELEM, ZRO1 1 . NEVAB,NGAUS,NSTR1,NTOTG,EPSTN,EFFST, ZRO1 2 . NTOTV,NVFIX,STRSG,TDISP,TFACT, ZRO1 3 . TLOAD,TREAC,MVFIX) ZRO1 4 C******************************* C******************************* C******************************* C**** THIS SUBROUTINE INITIALISES VARIOUS ARRAYS TO ZERO ZRO1 7 C ZRO1 8 C******************************* DIMENSION ELOAD(MELEM,MEVAB),STRSG(4,MTOTG),TDISP(MTOTV), ZRO1 10 . TLOAD(MELEM,MEVAB),TREAC(MVFIX,2),EPSTN(MTOTG), ZRO1 11 . EFFST(MTOTG) ZRO1 12 TFACT=0.0 ZRO1 13 DO 30 IELEM=1,NELEM ZRO1 14 DO 30 IEVAB=1,NEVAB ZRO1 15 ELOAD(IELEM,IEVAB)=0.0 ZRO1 16 ```

30	TLOAD(IELEM,IEVAB)=0.0	ZRO1	17
	DO 40 ITOTV=1,NTOTV	ZRO1	18
40	TDISP(ITOTV)=0.0	ZRO1	19
	DO 50 IVFIX=1,NVFIX	ZRO1	20
	DO 50 IDOFN=1,NDOFN	ZRO1	21
50	TREAC(IVFIX,IDOFN)=0.0	ZRO1	22
	DO 60 ITOTG=1,NTOTG	ZRO1	23
	EPSTN(ITOTG)=0.0	ZRO1	24
	EFFST(ITOTG)=0.0	ZRO1	25
	DO 60 ISTR1=1,NSTR1	ZRO1	26
60	STRSG(ISTR1,ITOTG)=0.0	ZRO1	27
	RETURN	ZRO1	28
	END	ZRO1	29

# 7.8.3 Subroutine INVAR The role of this subroutine is to evaluate the various functions of stress used to indicate either initiation of or continuing plastic deformation for the four yield criteria considered in this text. More explicitly we need to calculate the items listed in Table 7.2. Table 7.2 Effective stress and uniaxial yield stress levels for the yield criteria included in the elasto-plastic computer code.

Equation No.	Yield criterion	Stress level (effective stress)	Uniaxial (or equivalent yield stress)
(7.63)	Tresca	$2(J_2')^{1/2} \cos \theta$	$\sigma_Y$
(7.64)	Von Mises	$\sqrt{3} (J_2')^{1/2}$	$\sigma_Y$
(7.65)	Mohr–Coulomb	$\frac{1}{3} J_1 \sin \phi + (J_2')^{1/2} \times (\cos \theta - \sin \theta \sin \phi/\sqrt{3})$	$c \cos \phi$
(7.66)	Drucker–Prager	$a J_1 + (J_2')^{1/2}$	$k'$

Whether or not plastic deformation takes place at any point is governed by its stress level as monitored by the functions in the third column of Table 7.2. For plastic flow to occur this stress level must achieve the values given in the final column of Table 7.2. For the Tresca and Von Mises criteria this value is precisely the uniaxial yield stress but for the Mohr–Coulomb and Drucker–Prager criteria it is an equivalent value defined by the stress-independent terms in (7.65) and (7.66) respectively. Note that all the values given in the final column of Table 7.2 can be functions of the hardening parameter, $\kappa$ . Subroutine INVAR merely computes the effective or deviatoric stress components and then evaluates the appropriate function in the third column of Table 7.2 depending on the yield criterion being employed. The choice of yield criterion is governed by the parameter NCRIT, input in subroutine INPUT, and the available options are provided below ```txt NCRIT = 1 Tresca yield criterion 2 Von Mises 3 Mohr-Coulomb 4 Drucker-Prager ``` Subroutine INVAR is now presented and descriptive notes provided. ```csv SUBROUTINE INVAR(DEVIA,LPROP,MMATS,NCRIT,PROPS,SINT3,STEFF,STEMP,INVR 1 THETA,VARJ2,YIELD) INVR 2 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 DEVIA(4),PROPS(MMATS,7),STEMP(4) INVR 9 ROOT3=1.73205080757 INVR 10 SMEAN=(STEMP(1)+STEMP(2)+STEMP(4))/3.0 INVR 11 DEVIA(1)=STEMP(1)-SMEAN INVR 12 DEVIA(2)=STEMP(2)-SMEAN INVR 13 DEVIA(3)=STEMP(3) INVR 14 DEVIA(4)=STEMP(4)-SMEAN INVR 15 VARJ2=DEVIA(3)*DEVIA(3)+0.5*(DEVIA(1)*DEVIA(1)+DEVIA(2)*DEVIA(2) INVR 16 .+DEVIA(4)*DEVIA(4)) INVR 17 VARJ3=DEVIA(4)*(DEVIA(4)*DEVIA(4)-VARJ2) INVR 18 STEFF=SQRT(VARJ2) INVR 19 IF(STEFF.EQ.0.0) GO TO 10 INVR 20 -SINT3=-3.0*ROOT3*VARJ3/(2.0*VARJ2*STEFF) INVR 21 IF(SINT3.GT.1.0) SINT3=1.0 INVR 22 GO TO 20 INVR 23 10 SINT3=0.0 INVR 24 20 CONTINUE INVR 25 IF(SINT3.LT.-1.0) SINT3=-1.0 INVR 26 IF(SINT3.GT.1.0) SINT3=1.0 INVR 27 THETA=ASIN(SINT3)/3.0 INVR 28 GO TO (1,2,3,4) NCRIT INVR 29 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 3 PHIRA=PROPS(LPROP,7)*0.017453292 INVR 37 SNPHI=SIN(PHIRA) INVR 38 YIELD=SMEAN*SNPHI+STEFF*(COS(THETA)-SIN(THETA)*SNPHI/ROOT3) INVR 39 RETURN INVR 40 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 1 YIELD=2.0*COS(THETA)*STEFF INVR 31 RETURN INVR 32 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 ``` INVR 11-15 Compute the deviatoric stresses according to (7.7) with the order of the components being as indicated in (7.72). INVR 16–17 Calculate the second deviatoric stress invariant, $J_{2}^{\prime}$ . INVR 18 Calculate the third deviatoric stress invariant, $J_{3}'$ .