Sources of nonlinearity

The discretized principle of virtual work is recalled here, from our blog entitled "Mechanical finite element simulations explained",

\begin{align} \bigg\lbrace \int_{\Omega_h} \mathbf{B}^\text{T} \{{\pmb{\sigma}}\} \, \text{d}V \bigg\rbrace - \bigg\lbrace \int_{\Omega_h} \mathbf{N}^\text{T} \mathbf{b} \, \text{d}V \bigg\rbrace + \bigg\lbrace \int_{\partial\Omega_h,N } \mathbf{N}^\text{T} \overline{\mathbf{t}} \, \text{d}A \bigg\rbrace = \mathbf{\{0\}}, \end{align}

Three major sources of nonlinearity are potentially identified in Eq. 1:

• Material nonlinearity. The relationship between the stresses and the strains is nonlinear, i.e. \pmb{\sigma} \neq \mathbb{D}^\text{e} : \pmb{\varepsilon}, where \mathbb{D}^\text{e} is the constant elastic stiffness tensor; "constant" here meaning that \mathbb{D}^\text{e} does not depend on \pmb{\varepsilon}, or in other words, that it does not follow Hooke's law. A good example is plastic materials.

• Traction nonlinearity. The relationship between the traction and the displacements is nonlinear, i.e. \overline{\mathbf{t}} \neq \mathbf{A} \mathbf{u} + \overline{\mathbf{t}}_0, for some constant second-order tensor \mathbf{A} and constant prescribed traction \overline{\mathbf{t}}_0; constant in this context meaning independent of \mathbf{u}, but potentially dependent on time and/or space. Good examples are springs with non-constant stiffness, or contact.

If none of these nonlinearities is considered, then Eq. 1 is rendered linear, and its solution process is greatly simplified by using \{ \pmb{\sigma}\} = [ \mathbb{D} ] \{ \pmb{\varepsilon}\} = [ \mathbb{D} ] \mathbf{B} \{ \mathbf{u} \}, yielding

\begin{align*}\int_{\Omega_h} ( \mathbf{B}^\text{T} [ \mathbb{D}^\text{e} ] \mathbf{B} ) \, \text{d}V \{ \mathbf{u} \} - \int_{\Omega_h} \mathbf{N}^\text{T} \mathbf{b} \, \text{d}V - \int_{\partial\Omega_{h,N}} \mathbf{N}^\text{T} \overline{\mathbf{t}} \, \text{d}V = \{ \textbf{0} \} ,\end{align*}

and whose solution can be obtained as

\begin{align*} { \mathbf{\{u\}} = \left[ \int_{\Omega_h} ( \mathbf{B}^\text{T} [ \mathbb{D}^\text{e} ] \mathbf{B} ) \, \text{d}V \right]^{-1} \bigg\lbrace \int_{\Omega_h} \mathbf{N}^\text{T} \mathbf{b} \, \text{d}V + \int_{\partial\Omega_{h,N}} \mathbf{N}^\text{T} \overline{\mathbf{t}} \, \text{d}V } \bigg\rbrace. \end{align*}

Material nonlinearity

Linearization of the principle of virtual work

Material nonlinearity is the main source of nonlinearity under the assumption of small strains, and we will focus on it in the following. Assuming that the body \Omega is initially unstressed and that the displacements at a certain load step n are \mathbf{u}_n, the displacements at the load step n+1, \mathbf{u}_{n+1}, are obtained by solving Eq. 1, where now the stresses depend on the strains in a nonlinear fashion. In order to solve the, now nonlinear, principle of virtual work, Eq. 1 needs to be linearized, leading to

\begin{align*} &\int_{\Omega_h} \left( \mathbf{B}^\text{T} \left[ \frac{\partial \pmb{\sigma}}{\partial \pmb{\varepsilon}} \right]_{n+1} \mathbf{B} \right) \, \text{d}V \, \{ \delta\mathbf{u} \} + \int_{\Omega_h} \mathbf{B}^\text{T} \{ \pmb{\sigma} \}_{n+1} \, \text{d}V \cdots \end{align*} \\ \begin{align} &\cdots - \int_{\Omega_h} \mathbf{N}^\text{T} \mathbf{b}_{n+1} \, \text{d}V - \int_{\partial\Omega_{h,N}} \mathbf{N}^\text{T} \overline{\mathbf{t}}_{n+1} \, \text{d}V = \{\textbf{0}\}, \quad \quad \quad \quad \quad \quad \quad \end{align}

where \left[ \mathbb{D} \right] = \left[ \frac{\partial \pmb{\sigma}}{\partial \pmb{\varepsilon}} \right] is the stiffness tensor. Equation 2 for the iterative increment of displacements vector { \delta\mathbf{u}}^{k+1} at the iteration k+1 is solved by using a Newton-Raphson scheme, leading to the system of linear equations

\begin{align*} \mathbf{K}^k \delta\mathbf{u}^{k+1} = - \mathbf{R}^k, \end{align*}

where

\begin{align*} \mathbf{K}^k &= \int_{\Omega_h} \left( \mathbf{B}^\text{T} \left[ \frac{\partial \pmb{\sigma}}{\partial \pmb{\varepsilon}} \right]^k_{n+1} \mathbf{B} \right) \, \text{d}V, \\ \mathbf{R}^k &= \int_{\Omega_h} \mathbf{B} ^\text{T} \{ \pmb{\sigma} \}^k_{n+1} \, \text{d}V - \int_{\Omega_h} \mathbf{N}^\text{T} \mathbf{b}_{n+1} \, \text{d}V - \int_{\partial\Omega_{h,N}} \mathbf{N}^\text{T} \overline{\mathbf{t}}_{n+1} \, \text{d}V, \\ \{ \delta\mathbf{u} \}^{k+1} &= - ( \mathbf{K}^k )^{-1} \mathbf{R}^k, \end{align*}

and

\begin{align*} \mathbf{u}^{k+1}_{n+1} = \mathbf{u}^k_{n+1} + \delta\mathbf{u}^{k+1}. \end{align*}

Plasticity

The archetype for material nonlinearity under a small strain setting is plasticity. Plasticity is the capability of a material to undergo irreversible deformations under the effect of external loads. In this section we will consider an associative plastic model with a von Mises yield surface and linear isotropic hardening (these concepts will be explained in the following). The uniaxial behavior of such model can be seen in Fig. 1. Note that, even if the hardening behaviour is linear, the overall stress-strain behaviour is nonlinear.

Figure 1. Unaxial stress-strain plot depicting plasticity with linear hardening \\ (\sigma_{y} = 240 MPa, H = 1000 MPa).

\\

In a small strain setting, plasticity leads to an additive decomposition of the strain (or strain rate \dot{\pmb{\varepsilon}} in the following) in elastic and plastic components, such that [2]

\begin{align} \dot{\pmb{\varepsilon}} = \dot{\pmb{\varepsilon}}^\text{e} + \dot{\pmb{\varepsilon}}^\text{p}. \end{align}

It is also convenient to define a series of internal variables, which in the case contemplated here is only a single scalar variable, \kappa, but for more complicated models it could include additional variables and/or of higher tensorial order (e.g. vectors, second-order tensors,     \dots). Internal variables are quantities that describe dissipative mechanisms, which in the general case of plasticity lead to the irreversible transformation of mechanical energy, causing permanent deformation. In the example we are considering, since only isotropic hardening is considered, \kappa describes the evolution of the radius of the yield surface (Eq. 4). Note that also kinematic hardening could be considered in a plastic model, which would imply a translation of the yield surface, and not just a dilation.

The limit in which the considered material stops being elastic and becomes plastic is the yield surface, and in the von Mises case it is defined as

\begin{align} \phi = \phi ( \pmb{\sigma}, \sigma_y) = q ( \pmb{\sigma} ) - \sigma_y, \end{align}

which is a function of both the stress and the yield stress. q = \sqrt{\frac{3}{2} \, \mathbf{s} : \mathbf{s}} is the von Mises stress, \mathbf{s}, is the deviatoric component of \pmb{\sigma},

\begin{align}\sigma_y = H \kappa,\end{align}

is the yield stress, and H is the isotropic constant hardening modulus. Equation 5 describes linear isotropic hardening, or as previously mentioned, the dilation of the yield surface. A graphic representation of the von Mises yield surface can be seen in Fig. 2.

Figure 2. von Mises yield surface in principal stress space. The dashed red line is the hydrostatic axis (i.e. \sigma_{1} = \sigma_{2} = \sigma_{3} ).

The evolution in time of both \pmb{\varepsilon}^\text{p} and \kappa, i.e. \dot{\pmb{\varepsilon}}^\text{p} and \dot{\kappa}, must also be defined to complete the definition of a plasticity model. If these evolution equations are chosen such that

\begin{align} \dot{\pmb{\varepsilon}}^\text{p} &= \dot{\gamma} \frac{\partial \phi}{\partial \pmb{\sigma}} = \dot{\gamma} \sqrt{\frac{3}{2}} \frac{\mathbf{s}}{|| \mathbf{s} ||}, \end{align}

\begin{align*}\dot{\kappa} &= - \dot{\gamma} \frac{\partial \phi}{\partial \sigma_y} = \dot{\gamma},\end{align*}

i.e. the plastic strain rate is normal to the yield surface in stress space, then the model is said to be associative; \dot{\gamma} in Eqs. 6 is the plastic multiplier, which determines whether plastic flow is occurring or not.

Numerical solution of the plasticity model

In order to solve Eqs. 3, 4, and 6, an elastic predictor/plastic corrector algorithm is considered. Applying this two-step algorithm to the time increment n+1 consists of the two following steps:

• Elastic trial step : In this step it is assumed that \dot{\gamma} = 0, and thus that the time increment n+1 is purely elastic, then for a given total strain increment \Delta\pmb{\varepsilon},

  \begin{align*} \pmb{\varepsilon}^{\text{e}\star}_{n+1} & = \pmb{\varepsilon}^\text{e}_n + \Delta\pmb{\varepsilon}, \\ \kappa^\star_{n+1} & = \kappa_n,\end{align*}

  where the \star superscript denotes a trial quantity. These trial quantities can be used to compute the trial stress and the trial yield stress, such that

  \begin{align*} \pmb{\sigma}^\star_{n+1} &= \mathbb{D}^\text{e} :\pmb{\varepsilon}^{\text{e}\star}_{n+1}, \\ \sigma^\star_{y,n+1} &= H \kappa^\star_{n+1}.\end{align*}

  The trial step is deemed valid if it lies within the yield surface, meaning that the time step is completely elastic, as evaluated by

  \begin{align*}\phi^\star = \phi ( \pmb{\sigma}^\star_{n+1}, \sigma^\star_{y,n+1} ) = q^\star_{n+1} - \sigma^\star_{y,n+1} \leq 0,\end{align*}

  leading to the accepted solution for time increment n+1,

  \begin{align*}\pmb{\varepsilon}^\text{e}_{n+1} &= \pmb{\varepsilon}^{\text{e}\star}_{n+1}, \\ \kappa_{n+1} &= \kappa^\star_{n+1}.\end{align*}

• Plastic corrector step: If the trial stress computed assuming fully elastic behaviour, as in the first step, is greater than the trial yield stress, then the material undergoes plastic deformation and the trial stress needs to be projected back to the updated yield surface. Mathematically speaking, if \phi^\star = \phi \left( \pmb{\sigma}^\star_{n+1}, \sigma^\star_{y,n+1} \right) = q^\star_{n+1} - \sigma^\star_{y,n+1} > 0, and thus \dot{\gamma} \neq 0, the following system of algebraic equations must be solved (a backward Euler time integration approach has been applied to Eqs. 3, 4, and 6 here):

  \begin{align} \begin{Bmatrix} \pmb{\varepsilon}^\text{e} - \pmb{\varepsilon}^{\text{e}\star}_{n+1} + \Delta\gamma \frac{\partial \phi}{\partial \pmb{\sigma}} \Big|_{n+1} \\ \kappa_{n+1} - \kappa^\star_{n+1} - \Delta\gamma \\ \phi \left( \pmb{\sigma}_{n+1}, \sigma_{y,n+1} \right) \end{Bmatrix} = \begin{Bmatrix} \pmb{\varepsilon}^\text{e} - \pmb{\varepsilon}^{\text{e}\star}_{n+1} + \Delta\gamma \sqrt{\frac{3}{2}} \frac{\mathbf{s}}{||\mathbf{s} ||} \\ \kappa_{n+1} - \kappa^\star_{n+1} - \Delta\gamma \\ q_{n+1} - H \kappa_{n+1} \end{Bmatrix} = \begin{Bmatrix} \mathbf{0} \\ 0 \\ 0\end{Bmatrix}, \end{align}

  for the unknowns

  \begin{align*} \begin{Bmatrix} \pmb{\varepsilon}^\text{e}_{n+1} \\ \kappa_{n+1} \\ \Delta\gamma \end{Bmatrix}, \end{align*}

• and which, in the associative linear isotropic hardening von Mises case, has the explicit solution

  \begin{align}\Delta\gamma = \frac{\phi^\star}{3G + H}.\end{align}

The elastic predictor/plastic corrector algorithm is summarized graphically for a general plastic material with positive hardening in Fig. 3. The updated elastic domain is larger than the initial elastic domain, and thus the material has hardened.

Figure 3. Elastic Predictor / plastic corrector scheme, where \mathbf{A} = \sigma_{y} in this specific case.

Consistent tangent modulus

The consistent tangent modulus is the expression

\begin{align*}\mathbb{D} = \frac{\partial \pmb{\sigma}_{n+1}}{\partial \pmb{\varepsilon}_{n+1}},\end{align*}

which in the case of the elastic trial step is the same as the elastic modulus, such that

\begin{align*} \mathbb{D} = \frac{\partial \pmb{\sigma}^\star_{n+1}}{\partial \pmb{\varepsilon}^\star_{n+1}} = \mathbb{D}^\text{e} = 2G \, \mathbb{P}_\text{dev} + K \, \mathbf{I} \otimes \mathbf{I} \end{align*}

where G and K are, respectively, the shear and bulk moduli; \mathbb{P}_\text{dev} is the fourth-order deviatoric projection tensor, such that, for instance, \mathbf{s} = \mathbb{P}_\text{dev} : \pmb{\sigma}; \mathbf{I} is the second-order identity tensor; and \otimes is the tensor product.

If plastic correction is needed, the consistent tangent modulus is generally obtained as the derivative of an implicit function defined by Eqs. 7. In the specific example considered here, the consistent tangent modulus can be obtained from the first equation in Eqs. 7 by applying the elastic modulus \mathbb{D}^\text{e} on both sides, and Eq. 8, leading to

\begin{align*} \mathbb{D} = \frac{\partial \pmb{\sigma}_{n+1}}{\partial \pmb{\varepsilon}^\star_{n+1}} &= \frac{\partial}{\partial \pmb{\varepsilon}^\star_{n+1}} \left( \pmb{\sigma}^\star_{n+1} - \Delta\gamma \, \mathbb{D}^\text{e} : \frac{\partial \phi}{\partial \pmb{\sigma}} \bigg|_{n+1} \right) \\ &= \frac{\partial}{\partial \pmb{\varepsilon}^\star_{n+1}} \left[ \left( \mathbb{D}^\text{e} - 6 G^2 \frac{\Delta\gamma}{q^\star_{n+1}} \mathbb{P}_\text{dev} \right) : \pmb{\varepsilon}^\star_{n+1} \right] \\ & = \mathbb{D}^\text{e} - 6 G^2 \frac{\Delta\gamma}{q^\star_{n+1}} \mathbb{P}_\text{dev} - 6 G^2 \frac{\partial}{\partial \pmb{\varepsilon}^\star_{n+1}} \left(\frac{\Delta\gamma}{q^\star_{n+1}} \right) \mathbb{P}_\text{dev} : \pmb{\varepsilon}^\star_{n+1} \\ &= \mathbb{D}^\text{e} - 6 G^2 \frac{\Delta\gamma}{q^\star_{n+1}} \mathbb{P}_\text{dev} \cdots \\ &\quad \quad \cdots - 6 G^2 ( \mathbb{P}_\text{dev} : \pmb{\varepsilon}^\star_{n+1} ) \otimes \left( \frac{1} {q^\star_{n+1}} \frac{\partial \Delta\gamma}{\partial \pmb{\varepsilon}^\star_{n+1}} - \frac{\Delta\gamma}{( q^\star_{n+1} )^2} \frac{\partial q^\star_{n+1}}{\partial \pmb{\varepsilon}^\star_{n+1}} \right),\end{align*}

where

\begin{align*} \frac{\partial \Delta\gamma}{\partial \pmb{\varepsilon}^\star_{n+1}} &= \frac{\partial}{\partial \pmb{\varepsilon}^\star_{n+1}} \left( \frac{q^\star_{n+1} - H \kappa^\star_{n+1}}{3G + H} \right) = \frac{1}{3G + H} \frac{\partial q^\star_{n+1}}{\partial \pmb{\varepsilon}^\star_{n+1}} = \frac{1}{3G + H} \sqrt{\frac{3}{2}} \frac{\mathbf{s}^\star_{n+1}}{||\mathbf{s}^\star_{n+1} ||}, \\ \frac{\partial q^\star_{n+1}}{\partial \pmb{\varepsilon}^\star_{n+1}} &= 2G \sqrt{\frac{3}{2}} \frac{\mathbf{s}^\star_{n+1}}{| |\mathbf{s}^\star_{n+1} ||}, \end{align*}

finally leading to

\begin{align*} \mathbb{D} = \mathbb{D}^\text{e} &- 6 G^2 \frac{\Delta\gamma}{q^\star_{n+1}} \mathbb{P}_\text{dev} \cdots \\ &\cdots - \frac{6 G^2}{q^\star_{n+1}} \sqrt{\frac{3}{2}} \left( \frac{1}{3G + H} - \frac{2G \Delta\gamma}{q^\star_{n+1}} \right) ( \mathbb{P}_\text{dev} : \pmb{\varepsilon}^\star_{n+1} ) \otimes \frac{\mathbf{s}^\star_{n+1}}{|| \mathbf{s}^\star_{n+1} ||}. \end{align*}

Case studies

The two case studies presented here are extracted from the books by Hill and de Souza et al.¹²

Internally pressurized cylinder

This case consists of a long metallic thick-walled cylinder subjected to a linearly increasing internal pressure. Since this problem is symmetric in three planes, only an eight of the total geometry needs to be used (Fig. 4).

Figure 4. Geometry of the internally pressurized cylinder.

The length of the cylinder is 1 m, with internal an external radii of, respectively, a = 0.1 and b = 0.2 m. The applied boundary conditions are the corresponding symmetries where the geometry has been cut plus an internal pressure that is linearly increased, from 0 MPa to 190 MPa. The material used in this case study is a linear elastic material and a von Mises plastic surface with perfect plasticity, i.e. no hardening; the parameters are given in Table 1 and a graphical depiction in Fig. 5.

Figure 5. Uniaxial stress-strain plot depicting perfect plasticity; the material parameters can be seen in Table 1.

In this problem, plastic yield propagates radially from r = a to r = b, and plastic collapse occurs when r = b, i.e. when the entire cylinder has yielded. The zero hardening modulus means that the cylinder expands indefinitely under no further pressure at this point. Hill derived a closed-form solution for this problem, which establishes a relationship between the applied internal pressure and the radius c of the advancing front of plastic yield.¹ This expression is

\begin{align}P = \frac{2 \sigma_y}{\sqrt{3}} \left[ \text{ln} \left( \frac{c}{a} \right) + \frac{1}{2} \left( 1 - \frac{c^2}{b^2} \right) \right]\end{align}

Substituting c = a and c = b in Eq. 9 gives, respectively, the initiation of plastic yield and complete plastic collapse of the cylinder:

\begin{align}P_0 &= \frac{\sigma_y}{\sqrt{3}} \left( 1 - \frac{a^2}{b^2} \right), \end{align}

\begin{align*}P_\text{lim} &= \frac{2 \sigma_y}{\sqrt{3}} \text{ln} \left( \frac{b}{a} \right).\end{align*}

Applying Eqs. 10 to the current case results in P_0 \approx 140 MPa and P_\text{lim} \approx 332.7 MPa. The radial displacement at the outer surface u_b is given by

\begin{align}u_b = \begin{cases}\frac{2 P b}{E ( \frac{b^2}{a^2} - 1 )} ( 1 - \nu^2 ) \quad &\text{if } P <; P_0, \\ \frac{2 \sigma_y c^2}{\sqrt{3} E b} ( 1 - \nu^2 ) \quad & \text{otherwise}. \end{cases} \end{align}

Equation 11 is compared to a finite element simulation performed in Solve by using a P2 tetrahedral mesh. The corresponding results can be seen in Fig. 6.

Figure 6. Pressure vs radial displacement at the outer face of the internally pressurized cylinder.

As one would expect, the finite element simulation is highly accurate, and thus totally consistent with the analytical results derived by Hill.¹ The von Mises distribution across the thickness is seen in Fig. 7. This von Mises stress distribution corresponds to                   P = 190 MPa, so the plastic yield front is almost at r=b, and the von Mises stresses correspond to the yield value.

Figure 7. von Mises stress distribution across the thickness of the internally pressurized cylinder.

Internally pressurized spherical shell

This case consists of a metallic thick-walled spherical shell subjected to a linearly increasing internal pressure. Since this problem is symmetric in three planes, only an eight of the total geometry needs to be used (Fig. 8).

Figure 8. Geometry of the internally pressurized spherical shell.

The internal an external radii are, respectively, a = 0.1 and b = 0.2 m. The applied boundary conditions are the corresponding symmetries where the geometry has been cut plus an internal pressure that is linearly increased, from 0 MPa to 320 MPa. The material used in this case study is the same as in the internally pressurized cylinder (Table 1).

Again in this problem, plastic yield propagates radially from the r = a to r = b, and plastic collapse occurs when r = b. Hill derived a closed-form solution for this problem, which establishes a relationship between the applied internal pressure and the radius c of the advancing front of plastic yield.¹ This expression is

\begin{align}P = 2 \sigma_y \text{ln} \left( \frac{c}{a} \right) + \frac{2 \sigma_y}{3} \left( 1 - \frac{c^3}{b^3} \right).\end{align}

Substituting c = a and c = b in Eq. 12 gives, respectively, the initiation of plastic yield and complete plastic collapse of the spherical shell:

\begin{align}P_0 & = \frac{2 \sigma_y}{3} \left( 1 - \frac{a^3}{b^3} \right), \end{align}

\begin{align*} P_\text{lim} & = 2 \sigma_y \text{ln} \left( \frac{b}{a} \right).\end{align*}

Applying Eq. 13 to the current case at hand results in P_0 = 140 MPa and                     P_\text{lim} \approx 332.7 MPa. The radial displacement at the outer surface u_b is given by

\begin{align} u_b = \begin{cases}\frac{3 P b}{2 E ( \frac{b^3}{a^3} - 1 )} ( 1 - \nu ) \quad &\text{if } P < P_0, \\ \frac{\sigma_y c^3}{E b^2} ( 1 - \nu ) \quad & \text{otherwise}. \end{cases} \end{align}

Equation 14 is compared to a finite element simulation performed in Solve by using a P2 tetrahedral mesh. The corresponding results can be seen in Fig. 9.

Figure 9. Pressure vs radial displacement at the outer face of the internally pressurized spherical shell.

The von Mises distribution across the thickness is seen in Fig. 10. This von Mises stress distribution corresponds to P = 320 MPa, so the plastic yield front is almost at r=b, and the von Mises stresses correspond to the yield value.

Figure 10. Von Mises stress distribution across the thickness of the internally pressurized spherical shell.

In both Fig. 6 and Fig. 9 the curve would plateau at the value of P_\text{lim} due to the lack of hardening. Since in both these simulations we are prescribing pressure, if a value           P > P_\text{lim} were applied, the simulations would fail for any time step corresponding to such load as there would not be a possible solution to the system.

Defining a plastic material in Solve

In order to define the material used in these case studies in Solve we need to modify the material component of the SimAPI Python code, as following:

 
material = on.CustomMaterial('my_material')
elastic_model = on.material_models.ElasticIsotropicLinear(
youngs_modulus = 210000000000, poissons_ratio =0.3)
plastic_model = on.material_models.PlasticVonMisesLinearHardening(
yield_stress = 240000000, hardening = 0.0)
material.set('elastic_model', elastic_model)
material.set('plastic_model', plastic_model)
material >> geometry.parts[0]

Coming up next!

The next blog will introduce the necessary background in nonlinear continuum mechanics in order to understand geometrical nonlinearity, and also expand on the latter.

References

• R. Hill. The Mathematical Theory of Plasticity. Oxford engineering science series. Clarendon Pres, 1983. ISBN: 9780198561620.
• E.A. de Souza Neta, D. Peric, and D.R.J. Owen. Computational Methods for Plasticity: Theory and Applications. Wiley, 2011. ISBN: 9781119964544.