4.27. Foam Plasticity Model

4.27.1. 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. [[1]] 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 [[2]].

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

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

where an additive decomposition of the strain rates such that \(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 [[1]] proposed a yield function of the form

(4.125)\[f=\frac{\bar{\sigma}^2}{a^2}+\frac{p^2}{b^2}-1,\]

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

(4.126)\[a = A_0+A_1\phi^{A_2},\]

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

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

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

where,

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

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

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

where \(\varepsilon_{\text{V}}\) is the engineering volume strain.

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

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

where \(\dot{\gamma}\) is the consistency multiplier found by enforcing the corresponding condition and

(4.128)\[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 \(\beta\) is introduced in (4.128) to enable associated (\(\beta=0\)), radial (\(\beta=1\)), or a linear combination of the two flow rules (\(0<\beta<1\)) to be used. The two direction vectors may written as,

(4.129)\[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}|},\]

(4.130)\[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 \(\bar{\sigma} - p\) plane along the pressure axis. To generalize the model, an offset yield surface centered at pressure \(p_0\) may be used in the form

(4.131)\[\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 \(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 \(p_0\).

4.27.2. Implementation

Like other more classical rate-independent plasticity models (e.g. Section 4.7.2), 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,

\[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,

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

The engineering volume strain may then be readily computed via \(\varepsilon^{n+1}_{\text{V}}=\exp\left(\varepsilon^{n+1}_{kk}\right)-1\). A trial solid volume fraction is then calculated, \(\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,

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

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

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

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

4.27.3. Verification

The foam plasticity model is verified through a hydrostatic compression tests. Material properties used for this test are presented in Table 4.41 and correspond to room temperature properties of the PMDI20 rigid polyurethane foam characterized in [[1]].

Table 4.41 Material properties and model parameters for the foam plasticity model used during verification testing.

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

4.27.3.1. Hydrostatic Compression

The response of the foam plasticity model to hydrostatic compression is investigated here. Specifically, a displacement of the form \(u_i=\lambda\) is imposed resulting in a total strain field of \(\varepsilon_{11}=\varepsilon_{22}=\varepsilon_{33}=\ln\left(1+\lambda\right)\) and the engineering volume strain is simply \(\varepsilon_{\text{V}}=\left(1+\lambda\right)^3-1\). Furthermore, the maximum solid volume fraction monotonically increases and may be found to be \(\phi=\phi_0\frac{1}{\left(1+\lambda\right)^3}\). The stress state undergoes a similar reduction and is given to \(\sigma_{ij}=-p\delta_{ij}\) and \(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,

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

where \(b\) is evaluated via (4.127) and is a function of the strain. The model may then be simply solved as,

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

The elastic strains then reduce to \(\varepsilon^{\text{e}}_{ij}=-\frac{p}{3K}\delta_{ij}\) and the plastic strains computed as \(\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 Fig. 4.123.

Fig. 4.123 Pressure vs. engineering volume strain (\(\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.

4.27.4. User Guide

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:

• The initial volume fraction of solid material in the foam, \(\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 \(\varphi\) of \(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, \(\beta\), is defined with the BETA command line. When \(\beta = 0.0\), the flow direction is given by the normal to the yield surface (associated flow). When \(\beta = 1.0\), the flow direction is given by the stress tensor. It is recommended that associated flow \(\beta = 0.0\) be used for most analyses.

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

Table 4.42 State Variables for FOAM PLASTICITY Model

Name	Description
ITER	iterations
EVOL	volumetric strain
PHI	phi, \(\phi\)
EQPS	equivalent plastic strain, \(\bar{\varepsilon}^{p}\)
PA	\(A\)
PB	\(B\)

---

[1] (1,2,3)

M. K. Neilsen, W. Y. Lu, B. Olsson, and T. Hinnerichs. A viscoplastic constitutive model for polyurethane foams. In Proceedings ASME 2006 International Mechanical Engineering Congress and Exposition. Chicago, IL, 2006. ASME.

[2]

L.J. Gibson and M.F. Ashby. Cellular Solids: Structure and Properties. Cambridge Solid State Science Series. Cambridge University Press, Cambridge, UK, second edition, 1997.