Add the Storåkers hyperelastic foam model (#900)

* Add the hyperelastic foam model (tensortrax)

* Update _foam.py

* Update _foam.py

* Check success of transverse-constraint in `ViewMaterial.plot()`

* rename the Foam material to Storakers

* Update _storakers.py

* Update __init__.py

* Update __init__.py

* Update _storakers.py

* Add the hyperelastic foam model also for JAX

* add tests for hyperelastic foam model

* Update _storakers.py

* Update jax.rst

* Update tensortrax.rst

* Test the foam model with other parameters

which will fail with incompressibility as initial guess in plotting

* Update ex08_shear.py

* Test value error in `umat.plot()`

* Update test_constitution.py

* Update _storakers.py

docs/felupe/constitution/autodiff/jax.rst

@@ -49,6 +49,7 @@ These material model formulations are defined by a strain energy density functio
 
    miehe_goektepe_lulei
    mooney_rivlin
+   storakers
    third_order_deformation
    yeoh
 
@@ -93,6 +94,8 @@ formulations in :class:`~felupe.constitution.jax.Material`.
 
 .. autofunction:: felupe.constitution.jax.models.hyperelastic.mooney_rivlin
 
+.. autofunction:: felupe.constitution.jax.models.hyperelastic.storakers
+
 .. autofunction:: felupe.constitution.jax.models.hyperelastic.third_order_deformation
 
 .. autofunction:: felupe.constitution.jax.models.hyperelastic.yeoh

docs/felupe/constitution/autodiff/tensortrax.rst

@@ -41,6 +41,7 @@ These material model formulations are defined by a strain energy density functio
    ogden
    ogden_roxburgh
    saint_venant_kirchhoff
+   storakers
    third_order_deformation
    van_der_waals
    yeoh
@@ -113,6 +114,8 @@ and :func:`~felupe.constitution.tensortrax.updated_lagrange` material formulatio
 
 .. autofunction:: felupe.constitution.tensortrax.models.hyperelastic.saint_venant_kirchhoff
 
+.. autofunction:: felupe.constitution.tensortrax.models.hyperelastic.storakers
+
 .. autofunction:: felupe.constitution.tensortrax.models.hyperelastic.third_order_deformation
 
 .. autofunction:: felupe.constitution.tensortrax.models.hyperelastic.van_der_waals
@@ -121,4 +124,4 @@ and :func:`~felupe.constitution.tensortrax.updated_lagrange` material formulatio
 
 .. autofunction:: felupe.constitution.tensortrax.models.lagrange.morph
 
-.. autofunction:: felupe.constitution.tensortrax.models.lagrange.morph_representative_directions
+.. autofunction:: felupe.constitution.tensortrax.models.lagrange.morph_representative_directions

examples/ex08_shear.py

@@ -72,6 +72,7 @@
 # solid body.
 
 import felupe.constitution.tensortrax as mat
+
 # import felupe.constitution.jax as mat
 
 umat = mat.Hyperelastic(

src/felupe/__init__.py

@@ -262,6 +262,7 @@
         ogden,
         ogden_roxburgh,
         saint_venant_kirchhoff,
+        storakers,
         third_order_deformation,
         total_lagrange,
         updated_lagrange,
@@ -289,6 +290,7 @@
         "ogden",
         "ogden_roxburgh",
         "saint_venant_kirchhoff",
+        "storakers",
         "third_order_deformation",
         "van_der_waals",
         "yeoh",

src/felupe/constitution/__init__.py

@@ -67,6 +67,7 @@
         ogden,
         ogden_roxburgh,
         saint_venant_kirchhoff,
+        storakers,
         third_order_deformation,
         van_der_waals,
         yeoh,
@@ -89,6 +90,7 @@
         "ogden",
         "ogden_roxburgh",
         "saint_venant_kirchhoff",
+        "storakers",
         "third_order_deformation",
         "van_der_waals",
         "yeoh",

src/felupe/constitution/_view.py

@@ -28,7 +28,7 @@ def check_stretches(stretches):
     "Check if any stretch is lower or equal zero."
 
     if np.any(stretches <= 0.0):
-        raise ValueError("All stretches must greater than 0.")
+        raise ValueError("All stretches must be greater than 0.")
 
     return
 
@@ -193,7 +193,7 @@ def __init__(
         self.statevars_included = np.prod(self.umat.x[-1].shape) > 0
         self.statevars = statevars
 
-    def uniaxial(self, stretches=None):
+    def uniaxial(self, stretches=None, **kwargs):
         """Normal force per undeformed area vs. stretch curve for a uniaxial
         deformation.
 
@@ -237,7 +237,16 @@ def fun(λ3):
 
         from scipy.optimize import root
 
-        λ2 = λ3 = root(fun, λ3).x
+        res = root(fun, λ3, **kwargs)
+
+        if not res.success:
+            λ2 = λ3 = np.ones_like(λ1)
+            res = root(fun, λ3, **kwargs)
+
+            if not res.success:
+                raise ValueError("Uniaxial deformation failed.")
+
+        λ2 = λ3 = res.x
         F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)
 
         if self.statevars_included:
@@ -258,9 +267,9 @@ def fun(λ3):
 
         return λ1, P[0, 0].ravel(), "Uniaxial"
 
-    def planar(self, stretches=None):
+    def planar(self, stretches=None, **kwargs):
         """Normal force per undeformed area vs stretch curve for a planar shear
-        incompressible deformation.
+        deformation.
 
         Parameters
         ----------
@@ -301,7 +310,16 @@ def fun(λ3):
 
         from scipy.optimize import root
 
-        λ3 = root(fun, λ3).x
+        res = root(fun, λ3, **kwargs)
+
+        if not res.success:
+            λ3 = np.ones_like(λ1)
+            res = root(fun, λ3, **kwargs)
+
+            if not res.success:
+                raise ValueError("Planar deformation failed.")
+
+        λ3 = res.x
         F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)
 
         if self.statevars_included:
@@ -322,9 +340,9 @@ def fun(λ3):
 
         return λ1, P[0, 0].ravel(), "Planar Shear"
 
-    def biaxial(self, stretches=None):
+    def biaxial(self, stretches=None, **kwargs):
         """Normal force per undeformed area vs stretch curve for a equi-biaxial
-        incompressible deformation.
+        deformation.
 
         Parameters
         ----------
@@ -364,7 +382,16 @@ def fun(λ3):
 
         from scipy.optimize import root
 
-        λ3 = root(fun, λ3).x
+        res = root(fun, λ3, **kwargs)
+
+        if not res.success:
+            λ3 = np.ones_like(λ1)
+            res = root(fun, λ3, **kwargs)
+
+            if not res.success:
+                raise ValueError("Biaxial deformation failed.")
+
+        λ3 = res.x
         F = eye * np.array([λ1, λ2, λ3]).reshape(1, 3, 1, -1)
 
         if self.statevars_included:

src/felupe/constitution/jax/models/hyperelastic/__init__.py

@@ -1,17 +1,20 @@
 from ._miehe_goektepe_lulei import miehe_goektepe_lulei
 from ._mooney_rivlin import mooney_rivlin
+from ._storakers import storakers
 from ._third_order_deformation import third_order_deformation
 from ._yeoh import yeoh
 
 __all__ = [
     "miehe_goektepe_lulei",
     "mooney_rivlin",
+    "storakers",
     "third_order_deformation",
     "yeoh",
 ]
 
 # default (stable) material parameters
 miehe_goektepe_lulei.kwargs = dict(mu=0, N=100, U=0, p=2, q=2)
 mooney_rivlin.kwargs = dict(C10=0, C01=0)
+storakers.kwargs = dict(mu=[0], alpha=[2], beta=[1])
 third_order_deformation.kwargs = dict(C10=0, C01=0, C11=0, C20=0, C30=0)
 yeoh.kwargs = dict(C10=0, C20=0, C30=0)

src/felupe/constitution/jax/models/hyperelastic/_storakers.py

@@ -0,0 +1,36 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FElupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
+"""
+from functools import wraps
+
+from jax.numpy import array, sqrt
+from jax.numpy import sum as asum
+from jax.numpy.linalg import eigvalsh
+
+from ....tensortrax.models.hyperelastic import storakers as storakers_docstring
+
+
+@wraps(storakers_docstring)
+def storakers(C, mu, alpha, beta):
+    λ1, λ2, λ3 = sqrt(eigvalsh(C))
+    J = λ1 * λ2 * λ3
+
+    μ = array(mu)
+    α = array(alpha)
+    β = array(beta)
+
+    return asum(2 * μ / α**2 * (λ1**α + λ2**α + λ3**α - 3 + (J ** (-α * β) - 1) / β))

src/felupe/constitution/tensortrax/models/hyperelastic/__init__.py

@@ -22,6 +22,7 @@
 from ._ogden import ogden
 from ._ogden_roxburgh import ogden_roxburgh
 from ._saint_venant_kirchhoff import saint_venant_kirchhoff
+from ._storakers import storakers
 from ._third_order_deformation import third_order_deformation
 from ._van_der_waals import van_der_waals
 from ._yeoh import yeoh
@@ -40,6 +41,7 @@
     "ogden",
     "ogden_roxburgh",
     "saint_venant_kirchhoff",
+    "storakers",
     "third_order_deformation",
     "van_der_waals",
     "yeoh",
@@ -58,6 +60,7 @@
 ogden.kwargs = dict(mu=[0, 0], alpha=[2, -2])
 ogden_roxburgh.kwargs = dict(r=100, m=1, beta=0, material=neo_hooke, mu=0)
 saint_venant_kirchhoff.kwargs = dict(mu=0.0, lmbda=0.0)
+storakers.kwargs = dict(mu=[0], alpha=[2], beta=[1])
 third_order_deformation.kwargs = dict(C10=0, C01=0, C11=0, C20=0, C30=0)
 van_der_waals.kwargs = dict(mu=0, beta=0, a=0, limit=100)
 yeoh.kwargs = dict(C10=0, C20=0, C30=0)

src/felupe/constitution/tensortrax/models/hyperelastic/_storakers.py

@@ -0,0 +1,122 @@
+# -*- coding: utf-8 -*-
+"""
+This file is part of FElupe.
+
+FElupe is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+FElupe is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with FElupe.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+from tensortrax.math import sum as tsum
+from tensortrax.math.linalg import det, eigvalsh
+
+
+def storakers(C, mu, alpha, beta):
+    r"""Strain energy function of the Storåkers isotropic hyperelastic
+    `foam <https://doi.org/10.1016/0022-5096(86)90033-5>`_ material formulation [1]_.
+
+    Parameters
+    ----------
+    C : tensortrax.Tensor or jax.Array
+        Right Cauchy-Green deformation tensor.
+    mu : list of float
+        List of moduli.
+    alpha : list of float
+        List of stretch exponents.
+    beta : list of float
+        List of coefficients for the degree of compressibility.
+
+    Notes
+    -----
+    The strain energy function is given in Eq. :eq:`psi-foam`
+
+    ..  math::
+        :label: psi-ogden
+
+        \psi = \sum_i \frac{2 \mu_i}{\alpha^2_i} \left[
+            \hat{\lambda}_1^{\alpha_i} +
+            \hat{\lambda}_2^{\alpha_i} +
+            \hat{\lambda}_3^{\alpha_i} - 3
+            + \frac{1}{\beta_i} \left( J^{-\alpha \beta} - 1 \right)
+        \right]
+
+    The sum of the moduli :math:`\mu_i` is equal to the initial shear modulus
+    :math:`\mu`, see Eq. :eq:`shear-modulus-foam`,
+
+    ..  math::
+        :label: shear-modulus-ogden
+
+        \mu = \sum_i \mu_i
+
+    and the initial bulk modulus is given in Eq. :eq:`bulk-modulus-foam`.
+
+    ..  math::
+        :label: bulk-modulus-ogden
+
+        K = \sum_i 2 \mu_i \left( \frac{1}{3} + \beta_i \right)
+
+    Examples
+    --------
+    First, choose the desired automatic differentiation backend
+
+    ..  pyvista-plot::
+        :context:
+
+        >>> # import felupe.constitution.jax as mat
+        >>> import felupe.constitution.tensortrax as mat
+
+    and create the hyperelastic material.
+
+    ..  pyvista-plot::
+        :context:
+
+        >>> import felupe as fem
+        >>>
+        >>> umat = mat.Hyperelastic(
+        ...     mat.models.hyperelastic.storakers,
+        ...     mu=[4.5 * (1.85 / 2), -4.5 * (-9.2 / 2)],
+        ...     alpha=[1.85, -9.2],
+        ...     beta=[0.92, 0.92],
+        ... )
+        >>> ax = umat.plot(
+        ...     ux=fem.math.linsteps([1, 2], 15),
+        ...     ps=fem.math.linsteps([1, 1], 15),
+        ...     bx=fem.math.linsteps([1, 1], 9),
+        ... )
+
+    ..  pyvista-plot::
+        :include-source: False
+        :context:
+        :force_static:
+
+        >>> import pyvista as pv
+        >>>
+        >>> fig = ax.get_figure()
+        >>> chart = pv.ChartMPL(fig)
+        >>> chart.show()
+
+    References
+    ----------
+    .. [1] B. Storåkers, "On material representation and constitutive branching in
+       finite compressible elasticity", Journal of the Mechanics and Physics of Solids,
+       vol. 34, no. 2. Elsevier BV, pp. 125–145, Jan. 1986. doi:
+       10.1016/0022-5096(86)90033-5.
+    """
+
+    λ2 = eigvalsh(C)
+
+    return tsum(
+        [
+            2 * μ / α**2 * (tsum(λ2 ** (α / 2)) - 3 + (det(C) ** (-α * β / 2) - 1) / β)
+            for μ, α, β in zip(mu, alpha, beta)
+        ]
+    )