avoids locking and contains no spurious mechanisms. The Heterosis element is implemented here using a hierarchical formulation described later. We have already considered elastic Mindlin plate finite element analysis in Chapter 6. Nonlinear Mindlin plate finite element analysis is now considered. # 9.2 Equilibrium equations # 9.2.1 Three-dimensional equilibrium equations Let us begin with the equilibrium equations of three-dimensional stress analysis. We will assume that, for convenience, no tractions are present on the boundary $\Gamma_{t}$ of the three-dimensional domain $\Omega$ . The virtual work equation may be expressed as $$ \int_ {\Omega} \{[ \delta \epsilon ] ^ {T} \sigma - [ \delta u ] ^ {T} b \} d \Omega = 0 \tag {9.1} $$ where the vector of virtual displacements in the x, y and z directions is $\delta u = [\delta u, \delta v, \delta w]^{T}$ , the vector of associated virtual strains is $\delta \epsilon = [\delta \epsilon_{x}, \delta \epsilon_{y}, \delta \epsilon_{z}, \delta \gamma_{xy}, \delta \gamma_{xz}, \delta \gamma_{yz}]^{T}$ , the vector of stress is $\sigma = [\sigma_{x}, \sigma_{y}, \sigma_{z}, \tau_{xy}, \tau_{xz}, \tau_{yz}]^{T}$ and the vector of applied body forces is $b = [b_{x}, b_{y}, b_{z}]^{T}$ . Displacements u are prescribed on boundary $\Gamma_{u}$ of domain $\Omega$ . The stress-strain relationships for an isotropic material are given as $$ \boldsymbol {D} = a _ {1} \left[ \begin{array}{c c c c c c} a _ {2} & a _ {3} & a _ {3} & 0 & 0 & 0 \\ a _ {3} & a _ {2} & a _ {3} & 0 & 0 & 0 \\ a _ {3} & a _ {3} & a _ {2} & 0 & 0 & 0 \\ 0 & 0 & 0 & a _ {4} & 0 & 0 \\ 0 & 0 & 0 & 0 & a _ {4} & 0 \\ 0 & 0 & 0 & 0 & 0 & a _ {4} \end{array} \right] \tag {9.2} $$ where $a_1 = E / (1 + \nu)(1 - 2\nu)$ , $a_2 = 1 - \nu$ , $a_3 = \nu$ and $a_4 = (1 - 2\nu) / 2$ . Note that $E$ is the elastic modulus and $\nu$ is Poisson's ratio. # 9.2.2 Mindlin plate equilibrium equations In Mindlin plate theory, the domain of interest $\Omega$ is of the special form $$ \Omega = \{(x, y, z) \in R ^ {3} \mid z \in [ - t / 2, t / 2 ], (x, y) \in A \in R ^ {2} \} \tag {9.3} $$ where t is the plate thickness which may be a function of x and y and A is the plate area. The boundary of A is denoted by $\Gamma$ . We also make the following set of assumptions: (i) Normals to the midsurface (i.e., $z = 0$ ) before deformation remain straight but not necessarily normal to the midsurface after deformation. If $\theta_x$ and $\theta_y$ are the rotations of the midsurface normal in the $xz$ - and $yz$ -plane respectively, then $$ \mathbf {u} = \left[ \begin{array}{l} u (x, y, z) \\ v (x, y, z) \\ w (x, y, z) \end{array} \right] = \left[ \begin{array}{c} - z \theta_ {x} (x, y) \\ - z \theta_ {y} (x, y) \\ w (x, y) \end{array} \right] \tag {9.4} $$ The sign convention is illustrated in Fig. (9.1). Right hand rotations $\bar{\theta}_{x}$ and $\bar{\theta}_{y}$ are defined by the expression $$ \left[ \begin{array}{l} \theta_ {x} \\ \theta_ {y} \end{array} \right] = \left[ \begin{array}{l l} 0 & 1 \\ - 1 & 0 \end{array} \right] \left[ \begin{array}{l} \bar {\theta} _ {x} \\ \bar {\theta} _ {y} \end{array} \right]. \tag {9.5} $$ It is usually more convenient to develop the theory in terms of $\theta_{x}$ and $\theta_{y}$ rather than $\bar{\theta}_{x}$ and $\bar{\theta}_{y}$ since the resulting algebra is greatly simplified. (ii) The normal stress $\sigma_z$ is assumed equal to zero. The virtual work statement may be expressed as $$ \int_ {\Omega} \left[ \delta \boldsymbol {\epsilon} ^ {\prime} \right] ^ {T} \boldsymbol {\sigma} ^ {\prime} d \Omega - \int_ {\Omega} \left[ \delta \boldsymbol {u} \right] ^ {T} \boldsymbol {b} d \Omega = 0 \tag {9.6} $$ in which $$ [ \delta \epsilon^ {\prime} ] = [ \delta \epsilon_ {x}, \delta \epsilon_ {y}, \delta \gamma_ {x y} | \delta \gamma_ {x z}, \delta \gamma_ {y z} ] ^ {T} = [ (\delta \epsilon_ {f}) ^ {T}, (\delta \epsilon_ {s}) ^ {T} ] ^ {T} $$ and $$ \sigma^ {\prime} = \left[ \sigma_ {x}, \sigma_ {y}, \sigma_ {z} \mid \tau_ {x z}, \tau_ {y z} \right] ^ {T} = \left[ \left(\sigma_ {f}\right) ^ {T}, \left(\sigma_ {s}\right) ^ {T} \right] ^ {T}. $$ Note that $$ \delta \epsilon_ {f} = z \left[ - \frac {\partial (\delta \theta_ {x})}{\partial x}, - \frac {\partial (\delta \theta_ {y})}{\partial y}, - \left(\frac {\partial (\delta \theta_ {x})}{\partial y} + \frac {\partial (\delta \theta_ {y})}{\partial x}\right) \right] ^ {T} \doteq z \delta \hat {\epsilon} _ {f} \tag {9.7} $$ • In Mindlin plate theory a reduced form of the constitutive relations is obtained by making $\sigma_{z}=0$ and subsequently eliminating $\epsilon_{z}$ . Thus $$ \sigma^ {\prime} = D ^ {\prime} \epsilon^ {\prime} $$ where for elastic isotropic situations $$ \boldsymbol {D} ^ {\prime} = \frac {E}{(1 - \nu^ {2})} \left[ \begin{array}{c c c c c c} 1 & \nu & 0 & & 0 & 0 \\ \nu & 1 & 0 & & 0 & 0 \\ 0 & 0 & \frac {(1 - \nu)}{2} & & 0 & 0 \\ \hline 0 & 0 & 0 & & \frac {(1 - \nu)}{2} & 0 \\ 0 & 0 & 0 & & 0 & \frac {(1 - \nu)}{2} \end{array} \right] = \left[ \begin{array}{c c} \boldsymbol {D} _ {f ^ {\prime}} & \boldsymbol {0} \\ \boldsymbol {0} & \boldsymbol {D} _ {s ^ {\prime}} \end{array} \right] $$ † Terms symbolised thus ( $^{ }$ ) denote quantities integrated over the thickness. and $$ \delta \epsilon_ {s} = \left[ \frac {\partial (\delta w)}{\partial x} - \delta \theta_ {x}, \frac {\partial (\delta w)}{\partial y} - \delta \theta_ {y} \right] ^ {T} = \delta \hat {\epsilon} _ {s}. \tag {9.8} $$ Using (9.7) and (9.8) we find that (9.6) can be rewritten as $$ \int_ {A} \int_ {- t / 2} ^ {t / 2} [ z (\delta \hat {\epsilon} _ {f}) ^ {T} \sigma_ {f} + (\delta \hat {\epsilon} _ {s}) ^ {T} \sigma_ {s} - (\delta \boldsymbol {u}) ^ {T} \boldsymbol {b} ] d z d A = 0 \tag {9.9} $$ This equation is adopted in the layered approach. After integration over the thickness of the plate (9.9) can be rewritten in the form $$ \int_ {A} \left[ \left(\delta \hat {\epsilon} _ {f}\right) ^ {T} \hat {\sigma} _ {f} - \left(\delta \epsilon_ {s}\right) ^ {T} \hat {\sigma} _ {s} - (\delta \boldsymbol {u}) ^ {T} \hat {\boldsymbol {b}} \right] d A = 0 \tag {9.10} $$ where $$ \hat {\sigma} _ {f} = \int_ {- t / 2} ^ {t / 2} z \sigma_ {f} d z $$ $$ \hat {\sigma} _ {s} = \int_ {- t / 2} ^ {t / 2} \sigma_ {s} d z $$ and $$ \hat {\boldsymbol {b}} = \int_ {- t / 2} ^ {t / 2} \boldsymbol {b} d z. $$ We interpret $\hat{\sigma}_{f} = [M_{x}, M_{y}, M_{xy}]^{T}$ as the bending moments and $\hat{\sigma}_{s} = [Q_{x}, Q_{y}]^{T}$ as the shear force. Usually we take $\hat{b} = [q, 0, 0]^{T}$ in which q is the lateral distributed loading acting on the plate. We use (9.10) in the non-layered plate formulation. 

text_image

z, w Qx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myy Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx Myx

Fig. 9.1 Sign convention for Mindlin plate theory. # 9.3 Discretisation # 9.3.1 Standard representation If we adopt a standard $C(0)$ finite element representation then the displacements can be written as $$ u = \sum_ {i = 1} ^ {n} N _ {i} d _ {i} \tag {9.11} $$ in which the shape function matrix is $N_{i} = N_{i}I_{3}$ and the vector of nodal values $d_{i} = [w_{i}, \theta_{xi}, \theta_{yi}]^{T}$ . The flexural strain displacement equations are given as $$ \delta \hat {\epsilon} _ {f} = \sum_ {i = 1} ^ {n} B _ {f i} \delta d _ {i} \tag {9.12} $$ in which $$ \boldsymbol {B} _ {f i} = \left[ \begin{array}{c c c} 0 & - \frac {\partial N _ {i}}{\partial x} & 0 \\ 0 & 0 & - \frac {\partial N _ {i}}{\partial y} \\ 0 & - \frac {\partial N _ {i}}{\partial y} & - \frac {\partial N _ {i}}{\partial x} \end{array} \right] $$ The shear strain displacement equations have the form $$ \delta \hat {\epsilon} _ {s} = \sum_ {i = 1} ^ {n} B _ {s i} \delta d _ {i} \tag {9.13} $$ in which $$ \boldsymbol {B} _ {s i} = \left[ \begin{array}{c c c} \frac {\partial N _ {i}}{\partial x} & - N _ {i} & 0 \\ \frac {\partial N _ {i}}{\partial y} & 0 & - N _ {i} \end{array} \right] $$ If we substitute (9.11)-(9.13) in (9.9) we obtain the expression $$ \sum_ {i = 1} ^ {n} \left[ \delta \boldsymbol {d} _ {i} \right] ^ {T} \left\{\int_ {A} \int_ {- t / 2} ^ {t / 2} \left[ \boldsymbol {B} _ {f i} \right] ^ {T} \sigma_ {f} z + \left[ \boldsymbol {B} _ {s i} \right] ^ {T} \sigma_ {s} - \left[ N _ {i} \right] ^ {T} \boldsymbol {b} \right] d z d A \Bigg \} = 0. \tag {9.14} $$ Since (9.14) must be true for any set of virtual displacements we obtain the expression $$ \int_ {A} \int_ {- t / 2} ^ {t / 2} \left[ \left[ \boldsymbol {B} _ {f i} \right] ^ {T} \sigma_ {f} z + \left[ \boldsymbol {B} _ {s i} \right] ^ {T} \sigma_ {s} - \left[ N _ {i} \right] ^ {T} \boldsymbol {b} \right] d z d A = 0 \tag {9.15} $$ or $\psi_{i}(d) = 0.$ We use (9.15) in the layered approach. If we integrate the terms in square brackets over the thickness of the plate then we obtain the following equation $$ \int_ {A} \left[ \left[ \boldsymbol {B} _ {f t} \right] ^ {T} \hat {\boldsymbol {\sigma}} _ {f} + \left[ \boldsymbol {B} _ {s i} \right] ^ {T} \hat {\boldsymbol {\sigma}} _ {s} - \left[ \boldsymbol {N} _ {i} \right] ^ {T} \hat {\boldsymbol {b}} \right] d A = 0 \tag {9.16} $$ or $\psi_{i}(\pmb {d}) = 0.$ Equation (9.16) is used in the nonlayered approach. Note that we obtain equations for the residual force vector $\psi_{i}(d)$ for every node in the finite element discretisation. When the stresses are nonlinear then both (9.15) and (9.16) are sets of nonlinear simultaneous equations. Contributions to the residual force vector $\psi = [\psi_1^T, \ldots, \psi_n^T]^T$ may be evaluated at the element level and then assembled to form $\psi$ . We may use any standard $C(0)$ two-dimensional isoparametric element. Several elements have been presented in the literature and it emerges that the most convenient one is the 8/9 node 'Heterosis' element of Hughes. $^{(1)}$ In the programs described later we use 4, 8 and 9-noded isoparametric quadrilateral elements (see Chapter 6), as well as the Heterosis element. Selective integration is adopted and this will be described later. # 9.3.2 Hierarchical formulation of the Heterosis element In the implementation of the Heterosis and the 9-node element a hierarchical formulation is adopted. The first 8 shape functions are borrowed from the 8-noded Serendipity element and the shape function for the central $9^{th}$ node is the bubble function $$ N _ {9} ^ {(e)} (\xi , \eta) = (1 - \xi^ {2}) (1 - \eta^ {2}) \tag {9.17} $$ which is already available from the quadratic Lagrangian element. This means that all variables associated with the central node are hierarchical in nature. In other words, they are departures from the interpolated Serendipity values. The hierarchical representation can be used for geometrical representation as well as for interpolating displacements. In order to implement the heterosis element we adopt a hierarchical formulation either by adding a stiff spring (large number) to the leading diagonal term of the stiffness matrix associated with the lateral displacement parameter for node 9, or by prescribing displacement at this centre node to zero. This has the effect of forcing w to behave as though it was represented by Serendipity quadratic shape functions. Thus the desired effect is achieved. It is worth noting that if no spring is added the element obtained is identical to the 9-noded Lagrangian element provided that care is taken in evaluating the consistent nodal forces. Furthermore if stiff springs are added to all the terms of the leading diagonal associated with node 9, then the element reverts to a Serendipity 8-noded element. For convenience, in the present case, when representing the geometry of the heterosis element, the x and y coordinate departures from the interpolated Serendipity values are taken as equal to zero. In other words, as Serendipity geometrical representation is adopted this distinction is only of importance when elements with curved boundaries are present. (N.B. This is automatically taken care of by a modified version of Subroutine, NODEXY described in Section 6.4.1). # 9.4 Solution of nonlinear equations # 9.4.1 Plasticity in layered plates For Mindlin plates we may assume that the yield function $F$ is a function of $\sigma_f$ , the direct stresses associated with flexure, but not of the transverse shear stresses $\sigma_s$ . The yield function $F$ is also a function of the hardening parameter, $H$ . When yielding occurs at some point, it is assumed that, unless unloading takes place, the stresses always remain on the yield surface so that $$ F (\sigma_ {f}, H) = 0 \tag {9.18} $$ Thus the incremental stress-strain relationship is given as $$ \left[ \begin{array}{c} d \sigma_ {f} \\ d \sigma_ {s} \end{array} \right] = \left[ \begin{array}{c c} \left(\boldsymbol {D} _ {e p ^ {\prime}}\right) _ {f} & \boldsymbol {0} \\ \boldsymbol {0} & \boldsymbol {D} _ {s ^ {\prime}} \end{array} \right] \left[ \begin{array}{c} d \boldsymbol {\epsilon} _ {f} \\ d \boldsymbol {\epsilon} _ {s} \end{array} \right] \tag {9.19} $$ or $d\pmb{\sigma}^{\prime} = \pmb{D}_{ep}{}^{\prime}d\pmb{\epsilon}^{\prime}$ in which $(D_{ep}^{\prime})_{f}$ is identical to $D_{ep}$ given in Chapter 7 for the elasto-plastic plane stress problem. Note that $D_{s}^{\prime}$ always remains elastic. Recall from equation (7.47) that $$ \left(\boldsymbol {D} _ {e p} ^ {\prime}\right) _ {f} = \boldsymbol {D} _ {f} ^ {\prime} - \frac {\boldsymbol {d} _ {D} \boldsymbol {d} _ {D} ^ {T}}{A + \boldsymbol {d} _ {D} ^ {T} \boldsymbol {a} ^ {\prime}} \tag {9.20} $$ where $$ \boldsymbol {a} ^ {\prime} = \left[ \frac {\partial F}{\partial \sigma_ {x}}, \frac {\partial F}{\partial \sigma_ {y}}, \frac {\partial F}{\partial \tau_ {x y}} \right] ^ {T} $$ $$ \begin{array}{l} \boldsymbol {d} _ {D} = \boldsymbol {D} _ {f} ^ {\prime} \boldsymbol {a} ^ {\prime} \\ A = - \frac {1}{\lambda} \frac {\dot {c} F}{\dot {c} H} d H \\ \end{array} $$ in which $\lambda$ is the proportionality constant. Here we cater for Von Mises and Tresca materials only. We can thus use a slightly modified version of the coding described in Chapter 7 when evaluating $(D_{rp'})_f$ and when testing for yielding etc. # 9.4.2 Solution of the nonlinear equilibrium equations for layered plates The incremental equilibrium equations for the plate can be written at some stage in the solution (i.e., at an iteration during a load increment) as $$ \psi (\boldsymbol {d} ^ {p}) + \boldsymbol {K} _ {T} (\boldsymbol {d} ^ {p}) \Delta \boldsymbol {d} ^ {p} = 0 \tag {9.21} $$ where $\psi$ is obtained from (9.15) and $K_{T}(d^{p})$ is the tangential stiffness matrix which may be approximated as $$ \boldsymbol {K} _ {T} (\boldsymbol {d} ^ {p}) = \int_ {A} \int_ {- t / 2} ^ {t / 2} \{[ \boldsymbol {B} _ {f} ] ^ {T} [ \boldsymbol {D} _ {e p ^ {\prime}} ] _ {f} \boldsymbol {B} _ {f} + [ \boldsymbol {B} _ {s} ] ^ {T} \boldsymbol {D} _ {s ^ {\prime}} \boldsymbol {B} _ {s} \} d z d A. \tag {9.22} $$ Since $[D_{ep}^{\prime}]_{f}$ is a function of z we may employ a numerical integration technique to evaluate the integral over the thickness of the plate. Here, we divide the plate into layers and use a mid-ordinate rule as described in Chapter 5 for the Timoshenko beam. We use a similar method to evaluate $\psi(d^{p})$ . Thus we have $$ \boldsymbol {K} _ {T} (\boldsymbol {d} ^ {p}) = \int_ {A} \{[ \boldsymbol {B} _ {f} ] ^ {T} [ \hat {\boldsymbol {D}} _ {e p} ] _ {f} \boldsymbol {B} _ {f} \dots [ \boldsymbol {B} _ {s} ] ^ {T} \hat {\boldsymbol {D}} _ {s} \boldsymbol {B} _ {s} ] \} d A \tag {9.23} $$ where $$ [ \hat {D} _ {e p} ] _ {f} = \int_ {- t / 2} ^ {t / 2} [ D _ {e p} ^ {\prime} ] _ {f} d z $$ and $$ \hat {\boldsymbol {D}} _ {s} = \int_ {- t / 2} ^ {t / 2} \boldsymbol {D} _ {s} ^ {\prime} d z. $$ We now use the standard procedure to solve (9.21). Instead of using $K_{T}(d^{p})$ we may use some previously calculated value of $K_{T}$ just as in the other applications. # 9.4.3 Plasticity in nonlayered plates In Chapter 5 we considered the elasto-plastic nonlayered analysis of Timoshenko beams in which we assumed that when the bending moment reaches the yield moment $M_{0}$ , the whole cross-section of the beam becomes plastic instantaneously. We noted that this is a convenient fiction as in reality there is always a gradual spread of plasticity over the depth of the beam. In elasto-plastic nonlayered Mindlin plate analysis we make a similar approximation. Here we assume that the yield function $\hat{F}$ is expressed as a function of the bending moments $\hat{\sigma}_{f}$ , but not of the shear forces $\hat{\sigma}_{s}$ . The yield function is also assumed to be a function of a hardening parameter $\hat{H}$ . During yield it is assumed that the stress resultants $\hat{\sigma}_{f}$ must remain on the yield surface so that $$ \hat {F} (\hat {\sigma} _ {f}, \hat {H}) = 0 \tag {9.24} $$ where for the Tresca and Von Mises materials under consideration $$ \hat {F} (\hat {\sigma} _ {f}, \hat {H}) = \int_ {- t / 2} ^ {t / 2} F (\sigma_ {f}, H) d z. \tag {9.25} $$ Therefore, although $\hat{F}$ replaces $F$ , $(M_x, M_y, M_{xy})$ replace $(\sigma_x, \sigma_y, \tau_{xy})$ and $M_0 = \sigma_0 t^2 / 4$ replaces $\sigma_0$ , everything else remains unchanged and we can again make use of the coding given in Chapter 7. The incremental stress-strain resultant relationships are given as $$ \left[ \begin{array}{c} d \hat {\boldsymbol {\sigma}} _ {f} \\ d \hat {\boldsymbol {\sigma}} _ {s} \end{array} \right] = \left[ \begin{array}{c c} [ \hat {\boldsymbol {D}} _ {e p} ] _ {f} & \mathbf {0} \\ \mathbf {0} & \hat {\boldsymbol {D}} _ {s} \end{array} \right] \left[ \begin{array}{c} d \hat {\boldsymbol {\epsilon}} _ {f} \\ d \hat {\boldsymbol {\epsilon}} _ {s} \end{array} \right] \tag {9.26} $$ in which $$ [ \hat {\boldsymbol {D}} _ {e p} ] _ {f} = \hat {\boldsymbol {D}} _ {f} - \frac {\hat {\boldsymbol {d}} _ {D} \hat {\boldsymbol {d}} _ {D} ^ {T}}{\hat {\boldsymbol {A}} + \hat {\boldsymbol {d}} _ {D} ^ {T} \hat {\boldsymbol {a}}} \tag {9.27} $$ in which $$ \hat {\boldsymbol {a}} = \left[ \frac {\partial \hat {F}}{\partial M _ {x}}, \frac {\partial \hat {F}}{\partial M _ {y}}, \frac {\partial \hat {F}}{\partial M _ {x y}} \right] ^ {T} $$ $$ \hat {\boldsymbol {d}} _ {D} = \hat {\boldsymbol {D}} _ {f} \hat {\boldsymbol {a}} $$ $$ \hat {A} = - \frac {1}{\lambda} \frac {\partial \hat {F}}{\partial \hat {H}} d \hat {H} $$ and $$ \hat {\boldsymbol {D}} _ {f} = \int_ {- t / 2} ^ {t / 2} \boldsymbol {D} _ {f} ^ {\prime} z d z. $$ Note also that $$ \hat {\boldsymbol {D}} _ {s} = \int_ {- t / 2} ^ {t / 2} \boldsymbol {D} _ {s} ^ {\prime} d z. $$ # 9.4.4 Solution of nonlinear equilibrium equations for nonlayered Mindlin plates For the nonlayered plates the equilibrium equations are identical to (9.21). Here, however, the tangential stiffness matrix is given as $$ \boldsymbol {K} _ {T} = \int_ {A} \left\{\left[ \boldsymbol {B} _ {f} \right] ^ {T} \left[ \hat {\boldsymbol {D}} _ {e p} \right] _ {f} \boldsymbol {B} _ {f} + \left[ \boldsymbol {B} _ {s} \right] ^ {T} \hat {\boldsymbol {D}} _ {s} \boldsymbol {B} _ {s} \right\} d A. \tag {9.28} $$ Apart from this modification the solution procedure is unchanged. # 9.4.5 Summary of solution procedures The solution procedures for elasto-plastic Mindlin plate analysis are summarised in Tables 9.1–9.3. The overall process is given in Table 9.1. The iteration loop is shown for the nonlayered and layered plates in Tables 9.2 and 9.3 respectively. Table 9.1 Equation solving technique for layered and nonlayered Mindlin plates

1	Begin new load increment, $f = f + \Delta f$ .
2	Set $\Delta f$ equal to the current load increment vector.
3	Set $d^{0}$ equal to 0 for the first increment or equal to the total displacement vector at the end of the last load increment.
4	Set $\psi^{0}$ equal to the residual force vector at the end of the last increment or equal to 0 for the first load increment.
5	Set $\psi^{0} = \psi^{0} + \Delta f$ .
6	Solve $\Delta d^{0} = -[K_{T}]^{-1} \psi^{0}$ .Use old or updated value $K_{T}$ .
7	Set $d^{1} = d^{0} + \Delta d^{0}$ .
8	Evaluate $\psi^{1}(d^{1})$ .
9	If solution has converged go to 11; otherwise continue.
10	Iterate until solution has converged.
11	If this is not the last increment go to 1; otherwise stop.

Table 9.2 The iteration loop for elasto-plastic nonlayered Mindlin plates.

1	Set iteration number i = 1.
2	Solve Δdi = - [KT]−1ψi. Use old or updated KT.
3	Set di+1 = di + Δdi.
4	For each Gauss point, evaluate the increments in strain resultants

$$ \Delta \hat {\epsilon} _ {f} ^ {i} = B _ {f} \Delta d ^ {i} $$ $$ \Delta \hat {\boldsymbol {\epsilon}} _ {s} ^ {i} = \boldsymbol {B} _ {s} \Delta \boldsymbol {d} ^ {i}. $$ Table 9.2—continued 5 Using the elastic rigidities estimate, at each Gauss point, the increments in stress resultants and hence the total stress resultants $$ \Delta \hat {\boldsymbol {\sigma}} _ {f} ^ {i} = \hat {\boldsymbol {D}} _ {f} \Delta \hat {\boldsymbol {\epsilon}} _ {f} ^ {i} \quad \text { hence } \quad \hat {\boldsymbol {\sigma}} _ {f} ^ {i + 1} = \hat {\boldsymbol {\sigma}} _ {f} ^ {i} + \Delta \hat {\boldsymbol {\sigma}} _ {f} ^ {i} $$ $$ \Delta \hat {\boldsymbol {\sigma}} _ {s} ^ {i} = \hat {\boldsymbol {D}} _ {s} \Delta \hat {\boldsymbol {\epsilon}} _ {s} ^ {i} \quad \text { hence } \quad \hat {\boldsymbol {\sigma}} _ {s} ^ {i + 1} = \hat {\boldsymbol {\sigma}} _ {s} ^ {i} + \Delta \hat {\boldsymbol {\sigma}} _ {s} ^ {i}. $$ 6 At each Gauss point, depending on the states of $\hat{\sigma}_{f}^{i}$ and $\hat{\sigma}_{f}^{i+1}$ , adjust $\hat{\sigma}_{f}^{i+1}$ to satisfy the yield criterion and preserve the normality condition. 7 Evaluate the residual force vector $$ \psi^ {i + 1} = \int_ {A} \left\{\left[ B _ {f} \right] ^ {T} \hat {\sigma} _ {f} + \left[ B _ {s} \right] ^ {T} \hat {\sigma} _ {s} \right\} d A - f. $$ 8 If the solution has converged, continue, otherwise set $i = i + 1$ and go to 2. 9 Move to next load increment. Table 9.3 The iteration loop for elasto-plastic layered Mindlin plates. 1 Set iteration number $i = 1$ . 2 Solve $\Delta d^i = -[K_T]^{-1}\psi^i$ . Use old or updated $K_T$ . 3 Set $d^{i+1} = d^i + \Delta d^i$ . 4 For each Gauss point in each layer evaluate the increment in strain $$ \Delta \epsilon_ {f} ^ {i} = z B _ {f} \Delta d ^ {i} $$ $$ \Delta \epsilon_ {s} ^ {i} = B _ {s} \Delta d ^ {i}. $$ 5 Estimate the increments in stress at each Gauss point in each layer using the elastic stress-strain matrix. Hence evaluate the total stress value. $$ \Delta \sigma_ {f} ^ {i} = \mathbf {D} _ {f} ^ {\prime} \Delta \epsilon_ {f} ^ {\prime}, \quad \sigma_ {f} ^ {i + 1} = \sigma_ {f} ^ {i} + \Delta \sigma_ {f} ^ {i} $$ $$ \Delta \sigma_ {s} ^ {i} = \mathbf {D} _ {s} ^ {\prime} \Delta \epsilon_ {s} ^ {\prime}, \quad \sigma_ {s} ^ {i + 1} = \sigma_ {s} ^ {i} + \Delta \sigma_ {s} ^ {i}. $$ 6 Depending on the states of $\sigma_f^i$ and $\sigma_f^{i+1}$ , adjust $\sigma_f^{i+1}$ to satisfy the yield criterion and preserve the normality condition. 7 Evaluate the stress resultants $\hat{\sigma}_{f}^{i+1}$ and $\hat{\sigma}_{s}^{i+1}$ at each Gauss point. 8 Evaluate the residual force vector $$ \psi^ {i + 1} = \int_ {A} \left\{\left[ \boldsymbol {B} _ {f} \right] ^ {T} \hat {\sigma} _ {f} + \left[ \boldsymbol {B} _ {s} \right] ^ {T} \hat {\sigma} _ {s} \right\} d A - \boldsymbol {f}. $$ 9 If the solution has converged continue, otherwise set $i = i + 1$ and go to 2. 10 Move to next load increment.