# SPDX-License-Identifier: GPL-3.0-or-later
import numpy as np
from ...math import cdya_ik, ddot, det, dot, dya, identity, inv, trace, transpose
from .._base import ConstitutiveMaterial
from ._volumetric_strain_energy import _volumetric_strain_energy
[docs]
class NeoHooke(ConstitutiveMaterial):
r"""Nearly-incompressible isotropic hyperelastic Neo-Hookean material formulation.
The strain energy density function of the Neo-Hookean material formulation is a
linear function of the trace of the isochoric part of the right Cauchy-Green
deformation tensor.
Parameters
----------
mu : float or None, optional
Shear modulus (default is None)
bulk : float or None, optional
Bulk modulus of the volumetric part of the strain energy density function, if
``volumetric`` is a kind. If None, there is no volumetric part. It must be None
if ``volumetric`` is an instance, which holds its own parameters (default is
None).
parallel : bool, optional
A flag to evaluate the gradient and the hessian in chunks of cells by a pool of
threads (default is False).
volumetric : str or VolumetricStrainEnergy, optional
The volumetric part of the strain energy density function :math:`U(J)`, given by
a kind with the bulk modulus, ``"volume-ratio"`` for :math:`U(J)=K(J-1)^2/2`,
``"stretch"`` for :math:`U(J)=9K(J^{1/3}-1)^2/2` or ``"log"`` for
:math:`U(J)=K\ln(J)^2/2`, or by an instance of
:class:`~felupe.VolumetricStrainEnergy`, e.g. of a custom subclass. Default is
``"volume-ratio"``.
Notes
-----
.. note::
At least one of the two material parameters must not be None.
In a nearly-incompressible constitutive framework the strain energy
density is an additive composition of an isochoric and a volumetric
part. While the isochoric part is defined on the distortional part of
the deformation gradient, the volumetric part of the strain
energy function is defined on the determinant of the deformation
gradient.
.. math::
\psi &= \hat{\psi}(\hat{\boldsymbol{C}}) + U(J)
\hat\psi(\hat{\boldsymbol{C}}) &= \frac{\mu}{2}
(\text{tr}(\hat{\boldsymbol{C}}) - 3)
with
.. math::
J &= \text{det}(\boldsymbol{F})
\hat{\boldsymbol{F}} &= J^{-1/3} \boldsymbol{F}
\hat{\boldsymbol{C}} &= \hat{\boldsymbol{F}}^T \hat{\boldsymbol{F}}
The volumetric part of the strain energy density function is a function of the
volume ratio, given by ``volumetric``, see :class:`~felupe.VolumetricStrainEnergy`,
e.g. for ``volumetric="volume-ratio"`` (default)
.. math::
U(J) = \frac{K}{2} (J - 1)^2
The first Piola-Kirchhoff stress tensor is evaluated as the gradient
of the strain energy density function. The hessian of the strain
energy density function enables the corresponding elasticity tensor.
.. math::
\boldsymbol{P} &= \frac{\partial \psi}{\partial \boldsymbol{F}}
\mathbb{A} &= \frac{\partial^2 \psi}{\partial \boldsymbol{F}\ \partial
\boldsymbol{F}}
A chain rule application leads to the following expression for the stress tensor.
It is formulated as a sum of the **physical**-deviatoric (not the mathematical
deviator!) and the physical-hydrostatic stress tensors.
.. math::
\boldsymbol{P} &= \boldsymbol{P}' + \boldsymbol{P}_U
\boldsymbol{P}' &= \frac{\partial \hat{\psi}}{\partial \hat{\boldsymbol{F}}} :
\frac{\partial \hat{\boldsymbol{F}}}{\partial \boldsymbol{F}} =
\bar{\boldsymbol{P}} - \frac{1}{3} (\bar{\boldsymbol{P}} : \boldsymbol{F})
\boldsymbol{F}^{-T}
\boldsymbol{P}_U &= \frac{\partial U(J)}{\partial J} \frac{\partial J}{\partial
\boldsymbol{F}} = U'(J) J \boldsymbol{F}^{-T}
with
.. math::
\frac{\partial \hat{\boldsymbol{F}}}{\partial \boldsymbol{F}} &= J^{-1/3} \left(
\boldsymbol{I} \overset{ik}{\otimes} \boldsymbol{I} - \frac{1}{3} \boldsymbol{F}
\otimes \boldsymbol{F}^{-T} \right)
\frac{\partial J}{\partial \boldsymbol{F}} &= J \boldsymbol{F}^{-T}
\bar{\boldsymbol{P}} &= J^{-1/3} \frac{\partial \hat{\psi}}{\partial
\hat{\boldsymbol{F}}}
With the above partial derivatives the first Piola-Kirchhoff stress
tensor of the Neo-Hookean material model takes the following form.
.. math::
\boldsymbol{P} = \mu J^{-2/3} \left( \boldsymbol{F} - \frac{1}{3} (
\boldsymbol{F} : \boldsymbol{F}) \boldsymbol{F}^{-T} \right) + U'(J) J
\boldsymbol{F}^{-T}
Again, a chain rule application leads to an expression for the elasticity tensor.
.. math::
\mathbb{A} &= \mathbb{A}' + \mathbb{A}_{U}
\mathbb{A}' &= \bar{\mathbb{A}} - \frac{1}{3} \left( (\bar{\mathbb{A}} :
\boldsymbol{F}) \otimes \boldsymbol{F}^{-T} + \boldsymbol{F}^{-T} \otimes
(\boldsymbol{F} : \bar{\mathbb{A}}) \right ) + \frac{1}{9} (\boldsymbol{F} :
\bar{\mathbb{A}} : \boldsymbol{F}) \boldsymbol{F}^{-T} \otimes
\boldsymbol{F}^{-T}
\mathbb{A}_{U} &= (U''(J) J + U'(J)) J \boldsymbol{F}^{-T} \otimes
\boldsymbol{F}^{-T} - U'(J) J \boldsymbol{F}^{-T} \overset{il}{\otimes}
\boldsymbol{F}^{-T}
with
.. math::
\bar{\mathbb{A}} = J^{-1/3} \frac{\partial^2 \hat\psi}{\partial
\hat{\boldsymbol{F}}\ \partial \hat{\boldsymbol{F}}} J^{-1/3}
With the above partial derivatives the (physical-deviatoric and
-hydrostatic) parts of the elasticity tensor associated
to the first Piola-Kirchhoff stress tensor of the Neo-Hookean
material model takes the following form, here for ``volumetric="volume-ratio"``
with :math:`U'(J)=K(J-1)` and :math:`U''(J)=K`.
.. math::
\mathbb{A} &= \mathbb{A}' + \mathbb{A}_{U}
\mathbb{A}' &= J^{-2/3} \left(\boldsymbol{I} \overset{ik}{\otimes} \boldsymbol{I}
- \frac{1}{3} \left( \boldsymbol{F} \otimes \boldsymbol{F}^{-T} +
\boldsymbol{F}^{-T} \otimes \boldsymbol{F} \right ) + \frac{1}{9} (\boldsymbol{F}
: \boldsymbol{F}) \boldsymbol{F}^{-T} \otimes \boldsymbol{F}^{-T} \right)
\mathbb{A}_{U} &= K J \left( (2J - 1) \boldsymbol{F}^{-T} \otimes
\boldsymbol{F}^{-T} - (J - 1) \boldsymbol{F}^{-T} \overset{il}{\otimes}
\boldsymbol{F}^{-T} \right)
Examples
--------
.. plot::
>>> import felupe as fem
>>>
>>> umat = fem.NeoHooke(mu=1.0, bulk=2.0)
>>> ax = umat.plot()
The volumetric part of the strain energy density function is given by a kind or by
an instance of :class:`~felupe.VolumetricStrainEnergy`.
.. plot::
>>> import felupe as fem
>>>
>>> umat = fem.NeoHooke(mu=1.0, bulk=2.0, volumetric="stretch")
>>> ax = umat.plot()
"""
def __init__(self, mu=None, bulk=None, parallel=False, volumetric="volume-ratio"):
self.parallel = parallel
self.mu = mu
self.bulk = bulk
# volumetric part of the strain energy density function U(J), given by a kind
# (with the bulk modulus) or by an instance
self._volumetric = volumetric
_volumetric_strain_energy(volumetric, bulk)
self.kwargs = {}
if self.mu is not None:
self.kwargs["mu"] = self.mu
if self.bulk is not None:
self.kwargs["bulk"] = self.bulk
# aliases for function, gradient and hessian
self.energy = self.function
self.stress = self.gradient
self.elasticity = self.hessian
# initial variables for calling
# ``self.gradient(self.x)`` and ``self.hessian(self.x)``
self.x = [np.eye(3), np.zeros(0)]
@property
def volumetric(self):
"""The volumetric part of the strain energy density function, see
:class:`~felupe.VolumetricStrainEnergy`. For a kind, it is created with the
current bulk modulus and it is None without a bulk modulus.
"""
return _volumetric_strain_energy(self._volumetric, self.bulk)
[docs]
def function(self, x):
"""Strain energy density function per unit undeformed volume of the Neo-Hookean
material formulation.
Parameters
----------
x : list of ndarray
List with the Deformation gradient ``F`` (3x3) as first item
"""
F = x[0]
mu = self.mu
volumetric = self.volumetric
J = det(F)
C = dot(transpose(F), F)
W = np.zeros_like(J)
if mu is not None:
W += mu / 2 * (J ** (-2 / 3) * trace(C) - 3)
if volumetric is not None:
W += volumetric(J)
return [W]
def _gradient(self, x):
"""Gradient of the strain energy density function per unit undeformed volume of
the Neo-Hookean material formulation, evaluated for a chunk of cells by
:meth:`~felupe.ConstitutiveMaterial.gradient`.
Parameters
----------
x : list of ndarray
List with the Deformation gradient ``F`` (3x3) as first item
"""
F, statevars = x[0], x[-1]
mu = self.mu
volumetric = self.volumetric
J = det(F)
iFT = transpose(inv(F, J))
P = np.zeros_like(F)
if mu is not None:
# "physical"-deviatoric (not math-deviatoric!) part of P
trC = ddot(F, F)
trC_3 = np.divide(trC, 3, out=trC)
np.multiply(trC_3, iFT, out=P)
np.add(F, -P, out=P)
Jm23 = np.power(J, -2 / 3, out=trC)
np.multiply(P, Jm23, out=P)
np.multiply(P, mu, out=P)
if volumetric is not None:
# "physical"-volumetric (not math-volumetric!) part of P
dUdJ = volumetric.dUdJ(J)
JiFT = np.multiply(J, iFT, out=iFT)
dUdF = np.multiply(dUdJ, JiFT, out=JiFT)
np.add(P, dUdF, out=P)
return [P, statevars]
def _hessian(self, x):
"""Hessian of the strain energy density function per unit undeformed volume of
the Neo-Hookean material formulation, evaluated for a chunk of cells by
:meth:`~felupe.ConstitutiveMaterial.hessian`.
Parameters
----------
x : list of ndarray
List with the Deformation gradient ``F`` (3x3) as first item
"""
F = x[0]
mu = self.mu
volumetric = self.volumetric
J = det(F)
iFT = transpose(inv(F, J))
A4 = np.zeros((*F.shape[:2], *F.shape[:2], *F.shape[-2:]), dtype=F.dtype)
trC = None
A4b = None
A4c = None
if mu is not None:
# "physical"-deviatoric (not math-deviatoric!) part of A4
eye = identity(F)
trC = ddot(F, F)
np.add(A4, cdya_ik(eye, eye), out=A4)
A4b = dya(F, iFT)
np.multiply(-2 / 3, A4b, out=A4b)
np.add(A4, A4b, out=A4)
np.add(A4, np.transpose(A4b, [2, 3, 0, 1, 4, 5]), out=A4)
A4b = dya(iFT, iFT, out=A4b)
trC_3 = np.divide(trC, 3, out=trC)
A4c = np.multiply(A4b, trC_3, out=A4c)
np.add(A4, np.transpose(A4c, [0, 3, 2, 1, 4, 5]), out=A4)
np.multiply(A4c, 2 / 3, out=A4c)
np.add(A4, A4c, out=A4)
np.multiply(mu, A4, out=A4)
Jm23 = np.power(J, -2 / 3, out=trC)
np.multiply(Jm23, A4, out=A4)
if volumetric is not None:
# "physical"-volumetric (not math-volumetric!) part of A4
if A4b is None:
A4b = dya(iFT, iFT, out=A4b)
# first and second derivatives of U(J), evaluated before J is overwritten
dUdJ = volumetric.dUdJ(J)
d2UdJdJ = volumetric.d2UdJdJ(J)
pJ = np.multiply(dUdJ, J)
J2 = np.power(J, 2, out=J)
d2UdJdJ_J2 = np.multiply(J2, d2UdJdJ, out=J2)
qJ = np.add(pJ, d2UdJdJ_J2, out=d2UdJdJ_J2)
A4c = np.multiply(qJ, A4b, out=A4c)
np.add(A4, A4c, out=A4)
np.multiply(-pJ, np.transpose(A4b, [0, 3, 2, 1, 4, 5]), out=A4c)
np.add(A4, A4c, out=A4)
return [A4]
