.. _models-foam_plasticity:

*********************
Foam Plasticity Model
*********************

Theory
======

The foam plasticity model was developed to describe the response of porous materials (like 
closed-cell polyurethane foams) exhibiting irreversible, elastic-plastic like responses through 
large deformations.  Such foams can exhibit significant plastic deviatoric and volumetric 
strains leading to permanent shape and volume changes, respectively.  The former behavior is quite typical 
of metals and corresponding theories are well established.  The latter response, however, is not typical of metals 
and a theory combining these two behaviors is needed.  Given these responses of interest, the foam plasticity 
model is well suited to use with metal foams and many closed-cell polymeric foams (e.g. polyurethane, polystyrene 
bead, etc.) subjected to large deformations.  As permanent strains are of interest, this model is not appropriate 
for use with flexible foams that return to their undeformed shape after loads are removed.

Specifically, the model developed by Neilsen et al. [:footcite:`mat:ref:neilsen4`] seeks to capture the response 
associated with three distinct deformation regimes.  First, when foams are initially compressed, they typically 
exhibit an elastic response.  After sufficient load is applied, a plateau of nearly constant stress over 
a large deformation region is noted as pores start to compress and cell walls undergo substantial deformation.  
Eventually, the various collapsed cells and walls begin to interact and a densification response with substantial 
hardening is observed.  Details of these deformation processes may be found in the text of Gibson and 
Ashby [:footcite:`mat:ref:gibson`].

Like other plasticity-based models, the incremental constitutive law for the foam plasticity model is written as,

.. math::

   \stackrel{\circ}{\sigma}_{ij}=\mathbb{C}_{ijkl}D^{\text{e}}_{kl},

where an additive decomposition of the strain rates such that :math:`D_{ij}=D^{\text{e}}_{ij}+D^{\text{p}}_{ij}` is assumed.  
To describe the inelastic response of the foams of interest, Neilsen and coworkers [:footcite:`mat:ref:neilsen4`] proposed 
a yield function of the form

.. math::
   :label: eqn:fp:theory:yield

   f=\frac{\bar{\sigma}^2}{a^2}+\frac{p^2}{b^2}-1,

where :math:`\bar{\sigma}` is the von Mises effective stress (:math:`\bar{\sigma}=\sqrt{\left(3/2\right)s_{ij}s_{ij}}` with :math:`s_{ij}` being 
the deviatoric stress) and :math:`p` being the hydrostatic pressure (:math:`p=\left(1/3\right)\sigma_{kk}`).  
In such a form, the initial yield surface forms an ellipsoid about the hydrostat.  The two denominators, :math:`a` and :math:`b`, are state variables 
capturing hardening effects and have the functional form of,

.. math::
   :label: eqn:fp-theory:aHard

   a = A_0+A_1\phi^{A_2},

.. math::
   :label: eqn:fp-theory:bHard

   b = \left\{\begin{array}{cc} B_0+B_1\phi^{B_2} & p \geq 0 \\ B_0 & p < 0 \end{array} \right.,

with :math:`A_0,~A_1,~A_2,~B_0,~B_1,` and :math:`B_2` being model parameters and :math:`\phi` being the maximum volume fraction of solid material 
obtained through the loading history and is defined as,

.. math::

   \phi = \max\left(\phi_c\right)

where,

.. math::

   \phi_c=\phi_0\frac{V_0}{V},

in which :math:`\phi_0` is the initial volume fraction of solid material and :math:`V_0` and :math:`V` are the initial volume and current volume, respectively, 
of the foam.  Put in terms of the deformation,

.. math::

   \phi_c = \phi_0\frac{1}{1+\varepsilon_{\text{V}}},

where :math:`\varepsilon_{\text{V}}` is the engineering volume strain.

To describe the inelastic plastic deformation, a non-associated flow rule is used.  Specifically,

.. math::

   D^{\text{p}}_{ij}=\dot{\gamma}g_{ij},

where :math:`\dot{\gamma}` is the consistency multiplier found by enforcing the corresponding condition and 

.. math::
   :label: eqn:fp-theory:flowRule

   g_{ij}=\frac{\left(1-\beta\right)g^{\text{a}}_{ij}+\beta g^{\text{r}}_{ij}}{|\left(1-\beta\right)g^{\text{a}}_{ij}+\beta g^{\text{r}}_{ij}|},

with the superscripts a and r being used to denote associated and radial flow directions, respectively.  
The model parameter :math:`\beta` is introduced in :eq:`eqn:fp-theory:flowRule` to enable associated (:math:`\beta=0`), radial (:math:`\beta=1`), or a 
linear combination of the two flow rules (:math:`0<\beta<1`) to be used.  The two direction vectors may written as,

.. math::
   :label: eqn:fp:theory:associated

   g^{\text{a}}_{ij} = \frac{\frac{\partial f}{\partial \sigma_{ij}}}{|\frac{\partial f}{\partial\sigma_{kl}}|} = \frac{\frac{3}{a^2}s_{ij}+\frac{2}{3b^2}p\delta_{ij}}{|\frac{3}{a^2}s_{ij}+\frac{2}{3b^2}p\delta_{ij}|},


.. math::
   :label: eqn:fp:theory:rad

   g^{\text{r}}_{ij} =\frac{\sigma_{ij}}{|\sigma_{kl}|}=\frac{\sigma_{ij}}{\sqrt{\sigma_{kl}\sigma_{kl}}}.

In alternative models of foam plasticity, the yield function is offset in the
:math:`\bar{\sigma} - p` plane along the pressure axis. To generalize the model,
an offset yield surface centered at pressure :math:`p_0` may be used in the form

.. math::
   :label: eqn:fp:theory:yield_offset

   \hat{f}=\frac{\bar{\sigma}^2}{a^2}+\frac{(p-p_0)^2}{b^2}-1.

Note that with this form of the yield function, the component of the flow
direction :math:`g^a_{ij}` is now non-associative. This component is associative
to the yield surface centered at the origin, but non-associative to the
yield surface centered at :math:`p_0`. 

Implementation
==============

Like other more classical rate-independent plasticity models (e.g. :numref:`models-elastic_plastic-implementation`), the foam 
plasticity model is implemented in a hypoelastic fashion using an elastic predictor-inelastic 
corrector scheme.  As such, a trial material state is calculated by assuming purely elastic deformations.  
The trial stress is given by,

.. math::

   T^{tr}_{ij}=T_{ij}^n+\Delta t\left(\lambda\delta_{ij}d_{kk}+2\mu d_{ij}\right),

and an updated logarithmic volume strain is given by,

.. math::

   \varepsilon^{n+1}_{kk}=\varepsilon^n_{kk}+\Delta td_{kk}.

The engineering volume strain may then be readily computed via :math:`\varepsilon^{n+1}_{\text{V}}=\exp\left(\varepsilon^{n+1}_{kk}\right)-1`.  
A trial solid volume fraction is then calculated, :math:`\phi^{tr}=\phi_0\frac{1}{1+\varepsilon^{n+1}_{\text{V}}}`, and compared 
to the previous maximum to obtain the maximum solid volume fraction over the loading history,

.. math::
   :label: eqn:fp-imp:phi

   \phi_c^{n+1}=\max\left(\phi^n,\phi^{tr}\right).

Equations :eq:`eqn:fp-theory:aHard` and :eq:`eqn:fp-theory:bHard` are evaluated using the volume fraction found in 
:eq:`eqn:fp-imp:phi`.  Using invariants of the trial stress state, the yield function :eq:`eqn:fp:theory:yield` is calculated.  
If :math:`f\leq0`, the loading is elastic and the trial solution is correct.  On the other hand, if :math:`f>0` a correction scheme is necessary 
to iterate and determine the inelastic solution.  To that end, by noting
:math:`\Delta T_{ij}=-\mathbb{C}_{ijkl}\Delta d^{\text{P}}_{kl}=-\Delta \gamma\mathbb{C}_{ijkl}g_{kl}` (with :math:`\Delta` being a correction 
increment), the consistency condition may be used to find,

.. math::

   \Delta\gamma=\frac{f}{\frac{\partial f}{\partial \sigma_{ij}}\mathbb{C}_{ijkl}g_{kl}},

where the fact that the strain (and therefore :math:`a` and :math:`b`) do not change over an increment.  The correction is repeated until :math:`f<\text{tol}`.

Verification
============

The foam plasticity model is verified through a hydrostatic compression tests.  Material properties used 
for this test are presented in :numref:`tab-fp-ver-props` and correspond to room temperature 
properties of the PMDI20 rigid polyurethane foam characterized in [:footcite:`mat:ref:neilsen4`].

.. _tab-fp-ver-props:

.. csv-table:: Material properties and model parameters for the foam plasticity model used during verification testing.
   :align: center
   :delim: &

   :math:`E` & 22,600 psi & :math:`\nu` & 0.343
   :math:`A_0` & 513.3 psi & :math:`A_1` & 4,629 psi
   :math:`A_2` & 2.90 & :math:`\phi_0` & 0.238
   :math:`B_0` & 971 psi & :math:`B_1` & 7,377.5 psi
   :math:`B_2` & 4.89 & :math:`\beta` & 0.95

Hydrostatic Compression
-----------------------

The response of the foam plasticity model to hydrostatic compression is investigated here.  Specifically, 
a displacement of the form :math:`u_i=\lambda` is imposed resulting in a total strain field of 
:math:`\varepsilon_{11}=\varepsilon_{22}=\varepsilon_{33}=\ln\left(1+\lambda\right)` and the engineering 
volume strain is simply :math:`\varepsilon_{\text{V}}=\left(1+\lambda\right)^3-1`.  Furthermore, the maximum solid 
volume fraction monotonically increases and may be found to be :math:`\phi=\phi_0\frac{1}{\left(1+\lambda\right)^3}`.  
The stress state undergoes a similar reduction and is given to :math:`\sigma_{ij}=-p\delta_{ij}` and :math:`s_{ij}=0`.  
Note, in the theory section pressure is defined positive in tension whereas in the current exercise 
compression is positive.  There is no practical impact of this selection and it is done for simplicity of 
presentation.  This simplification leads to a reduced yield function of the form,

.. math::

   f=\frac{p^2}{b^2}-1,
 

where :math:`b` is evaluated via :eq:`eqn:fp-theory:bHard` and is a function of the strain.  The model may then be simply solved as,

.. math::

   p=\left\{\begin{array}{cc} -3K\ln\left(1+\lambda\right) & f \leq 0 \\ b & f > 0 \end{array}\right. .


The elastic strains then reduce to :math:`\varepsilon^{\text{e}}_{ij}=-\frac{p}{3K}\delta_{ij}` and the plastic strains 
computed as :math:`\varepsilon^{\text{p}}_{ij}=\left(\ln\left(1+\lambda\right)+\frac{p}{3K}\right)\delta_{ij}`.  The resulting 
engineering strain vs. pressure results determined numerically and analytically are presented in :numref:`fig-fp-pressureStrain`.

.. _fig-fp-pressureStrain:

.. figure:: ../../_static/figures/hydroCompression.png
   :align: center
   :scale: 25 %

   Pressure vs. engineering volume strain (:math:`\varepsilon_{\text{V}}`) response of the foam plasticity model through a hydrostatic compression cycle.  For convenience, pressure is defined as positive in compression in the figure.

User Guide
==========

.. code-block:: sierrainput

   BEGIN PARAMETERS FOR MODEL FOAM_PLASTICITY
     #
     # Elastic constants
     #
     YOUNGS MODULUS = <real>
     POISSONS RATIO = <real>
     SHEAR MODULUS  = <real>
     BULK MODULUS   = <real>
     LAMBDA         = <real>
     TWO MU         = <real>
     #
     #
     #
     PHI             = <real>
     SHEAR STRENGTH  = <real>
     SHEAR HARDENING = <real>
     SHEAR EXPONENT  = <real>
     HYDRO STRENGTH  = <real>
     HYDRO HARDENING = <real>
     HYDRO EXPONENT  = <real>
     BETA            = <real>
     P0              = <real>
   END [PARAMETERS FOR MODEL FOAM_PLASTICITY]


In the above command blocks:

.. - See :numref:`mat:elasticConstants` for more information on elastic constants input.

- The initial volume fraction of solid material in the foam, :math:`\varphi`, is defined with the ``PHI`` command line. For example, solid polyurethane weighs 75 pounds per cubic foot (pcf); uncompressed 10 pcf polyurethane foam would have a :math:`\varphi` of :math:`0.133 = 10/75`.
- The shear (deviatoric) strength of uncompressed foam is defined with the ``SHEAR STRENGTH`` command line.
- The shear hardening modulus for the foam is defined with the ``SHEAR HARDENING`` command line.
- The shear hardening exponent is defined with the ``SHEAR EXPONENT`` command line. The deviatoric strength is given by ``(SHEAR STRENGTH)`` + ``(SHEAR HARDENING) * PHI**(SHEAR EXPONENT)``.
- The hydrostatic (volumetric) strength of the uncompressed foam is defined with the ``HYDRO STRENGTH`` command line.
- The hydrodynamic hardening modulus is defined with the ``HYDRO HARDENING`` command line.
- The hydrodynamic hardening exponent is defined with the ``HYDRO EXPONENT`` command line. The hydrostatic strength is given by ``(HYDRO STRENGTH)`` + ``(HYDRO HARDENING) * PHI**(HYDRO EXPONENT)``.
- The prescription for non associated flow, :math:`\beta`, is defined with the ``BETA`` command line. When :math:`\beta = 0.0`, the flow direction is given by the normal to the yield surface (associated flow). When :math:`\beta = 1.0`, the flow direction is given by the stress tensor. It is recommended that associated flow :math:`\beta = 0.0` be used for most analyses.

Output variables available for this model are listed in:numref:`out-tab-foamplasstvar`.

.. _out-tab-foamplasstvar:

.. csv-table:: State Variables for FOAM PLASTICITY Model
   :align: center
   :delim: &
   :header: Name, Description

   ``ITER`` & iterations 
   ``EVOL`` & volumetric strain 
   ``PHI``  & phi, :math:`\phi` 
   ``EQPS`` & equivalent plastic strain, :math:`\bar{\varepsilon}^{p}` 
   ``PA``   & :math:`A` 
   ``PB``   & :math:`B` 

.. raw::
   html

   <hr>

.. footbibliography::
