Implementation of finite Strain Viscoelastic Model with Mooney-Rivlin Hyperelasticity

README.md

@@ -203,6 +203,7 @@ Included [pseudo-elastic material models](https://github.com/adtzlr/matadi/blob/
 
 Included [viscoelastic material models](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_viscoelasticity.py):
 - [Finite-Strain-Viscoelastic](https://doi.org/10.1016/j.cma.2013.07.004) ([code](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_viscoelasticity.py#L4-L18))
+- [Finite-Strain-Viscoelastic (Mooney-Rivlin)](https://doi.org/10.1002/nme.5724) ([code](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_viscoelasticity.py#L21-L51))
 
 Included [other material models](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_misc.py):
 - [MORPH](https://doi.org/10.1016/S0749-6419(02)00091-8) ([code](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_misc.py#L19-L75))
@@ -259,8 +260,9 @@ d2WdFdF, d2WdFdp, d2Wdpdp = NH.hessian([defgrad, pressure, statevars])
 The Neo-Hooke, the MORPH and the Finite-Strain-Viscoelastic [[4](https://doi.org/10.1016/j.cma.2013.07.004)] material model formulations are available as ready-to-go materials in `matadi.models` as:
 
 * [`NeoHookeOgdenRoxburgh()`](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_templates.py), 
-* [`Morph()`](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_templates.py) and
-* [`Viscoelastic()`](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_templates.py).
+* [`Morph()`](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_templates.py),
+* [`Viscoelastic()`](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_templates.py) and
+* [`ViscoelasticMooneyRivlin()`](https://github.com/adtzlr/matadi/blob/main/src/matadi/models/_templates.py).
 
 **Hint**: *The state variable concept is also implemented for the `Material` class.*
 

src/matadi/_templates.py

@@ -204,8 +204,9 @@ def _fun_wrapper(self, x, **kwargs):
 
 
 class MaterialComposite:
+    "Composite Material as a sum of a list of hyperelastic materials."
+
     def __init__(self, materials):
-        "Composite Material as a sum of a list of hyperelastic materials."
         self.materials = materials
         self.fun = self.composite
 
@@ -231,10 +232,11 @@ def hessian(self, x, **kwargs):
 
 
 class MaterialTensorGeneral(MaterialTensor):
-    def __init__(self, fun, statevars_shape=(1, 1), x=None, triu=True, **kwargs):
-        """A (first Piola-Kirchhoff stress) tensor-based material definition with
-        state variables of a given shape."""
+    """A (first Piola-Kirchhoff stress) tensor-based material definition with
+    state variables of a given shape.
+    """
 
+    def __init__(self, fun, statevars_shape=(1, 1), x=None, triu=True, **kwargs):
         if x is None:
             x = [Variable("F", 3, 3)]
 

src/matadi/math.py

@@ -257,3 +257,29 @@ def astensor(A, scale=1):
 
     else:
         raise ValueError("Unknown shape of input.")
+
+
+def unimodular(T):
+    """
+    Compute the unimodular part of a tensor.
+
+    The unimodular part of a tensor is a modified version of the tensor where
+    the determinant is raised to the power of (-1/3) and multiplied to the tensor.
+    This operation preserves the isochoric (volume-preserving) part of the tensor
+    while removing the volumetric part.
+    """
+    return (det(T) ** (-1 / 3)) * T
+
+
+def sqrtm(C, eps=8e-5):
+    """
+    Compute the matrix square root of a tensor C using eigendecomposition.
+    """
+    w = eigvals(C, eps=eps)
+    eye = SX.eye(3)
+
+    M1 = (C - w[1] * eye) * (C - w[2] * eye) / (w[0] - w[1]) / (w[0] - w[2])
+    M2 = (C - w[2] * eye) * (C - w[0] * eye) / (w[1] - w[2]) / (w[1] - w[0])
+    M3 = (C - w[0] * eye) * (C - w[1] * eye) / (w[2] - w[0]) / (w[2] - w[1])
+
+    return sqrt(w[0]) * M1 + sqrt(w[1]) * M2 + sqrt(w[2]) * M3

src/matadi/models/__init__.py

@@ -19,8 +19,13 @@
 )
 from ._misc import morph
 from ._pseudo_elasticity import ogden_roxburgh
-from ._templates import Morph, NeoHookeOgdenRoxburgh, Viscoelastic
-from ._viscoelasticity import finite_strain_viscoelastic
+from ._templates import (
+    Morph,
+    NeoHookeOgdenRoxburgh,
+    Viscoelastic,
+    ViscoelasticMooneyRivlin,
+)
+from ._viscoelasticity import finite_strain_viscoelastic, finite_strain_viscoelastic_mr
 from .microsphere.nonaffine import miehe_goektepe_lulei
 
 __all__ = [
@@ -47,5 +52,6 @@
     "NeoHookeOgdenRoxburgh",
     "Viscoelastic",
     "finite_strain_viscoelastic",
+    "finite_strain_viscoelastic_mr",
     "miehe_goektepe_lulei",
 ]

src/matadi/models/_templates.py

@@ -3,7 +3,7 @@
 from ._hyperelasticity_isotropic import neo_hooke
 from ._misc import morph
 from ._pseudo_elasticity import ogden_roxburgh
-from ._viscoelasticity import finite_strain_viscoelastic
+from ._viscoelasticity import finite_strain_viscoelastic, finite_strain_viscoelastic_mr
 
 
 class NeoHookeOgdenRoxburgh(MaterialTensorGeneral):
@@ -81,3 +81,27 @@ def fun(x, mu, eta, dtime):
             eta=eta,
             dtime=dtime,
         )
+
+
+class ViscoelasticMooneyRivlin(MaterialTensorGeneral):
+    "Finite strain viscoelastic material formulation with Mooney-Rivlin hyperelasticity."
+
+    def __init__(
+        self,
+        c10=1,
+        c01=1,
+        eta=1,
+        dtime=1,
+    ):
+        def fun(x, c10, c01, eta, dtime):
+            P, statevars = finite_strain_viscoelastic_mr(x, c10, c01, eta, dtime)
+            return P, statevars
+
+        super().__init__(
+            fun=fun,
+            statevars_shape=(6, 1),
+            c10=c10,
+            c01=c01,
+            eta=eta,
+            dtime=dtime,
+        )

src/matadi/models/_viscoelasticity.py

@@ -1,4 +1,4 @@
-from ..math import astensor, asvoigt, det, gradient, inv, trace
+from ..math import astensor, asvoigt, det, eye, gradient, inv, sqrtm, trace, unimodular
 
 
 def finite_strain_viscoelastic(x, mu, eta, dtime):
@@ -16,3 +16,36 @@ def finite_strain_viscoelastic(x, mu, eta, dtime):
 
     # first Piola-Kirchhoff stress tensor and state variable
     return gradient(mu / 2 * (I1 - 3), F), asvoigt(Ci)
+
+
+def finite_strain_viscoelastic_mr(x, c10, c01, eta, dtime):
+    """
+    Finite strain viscoelastic material formulation with Mooney-Rivlin hyperelasticity
+    (Shutov 2018) https://doi.org/10.1002/nme.5724
+    """
+
+    # Split input into the deformation gradient and the vector of state variables
+    F, Cin = x[0], x[-1]
+
+    # Right cauchy-green deformation tensor
+    C = F.T @ F
+
+    # Based on
+    # <<TABLE 1: Iteration-free Euler backward method on the reference configuration>>
+    A = (
+        unimodular(sqrtm(inv(C)))
+        @ (astensor(Cin) + (dtime / eta) * c10 * unimodular(C))
+        @ unimodular(sqrtm(inv(C)))
+    )
+    eps = c01 * (dtime / eta)
+    phi0 = det(A) ** (1 / 3)
+    phi = phi0 - (trace(A) / (3 * phi0)) * eps
+    X = 2 * A @ inv(sqrtm(phi * phi * eye(3) + 4 * eps * A) + phi * eye(3))
+
+    Ci = unimodular(sqrtm(C) @ X @ sqrtm(C))
+
+    I1 = trace(unimodular(C @ inv(Ci)))
+    I2 = trace(unimodular(Ci @ inv(C)))
+
+    # First Piola-Kirchhoff stress tensor and state variable
+    return gradient(c10 / 2 * (I1 - 3) + c01 / 2 * (I2 - 3), F), asvoigt(Ci)

tests/test_template-materials.py

@@ -6,6 +6,7 @@
     Morph,
     NeoHookeOgdenRoxburgh,
     Viscoelastic,
+    ViscoelasticMooneyRivlin,
     neo_hooke,
     volumetric,
 )
@@ -47,7 +48,12 @@ def test_templates():
 def test_templates_models():
     # Material as a function of `F`
     # with additional state variables `z`
-    for M in [NeoHookeOgdenRoxburgh(), Morph(), Viscoelastic()]:
+    for M in [
+        NeoHookeOgdenRoxburgh(),
+        Morph(),
+        Viscoelastic(),
+        ViscoelasticMooneyRivlin(),
+    ]:
         FF = (np.random.rand(3, 3, 8, 100) - 0.5) / 2
         zz = np.random.rand(*M.x[-1].shape, 8, 100)