diff --git a/README.md b/README.md index e2600e1..58d23de 100644 --- a/README.md +++ b/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.* diff --git a/src/matadi/_templates.py b/src/matadi/_templates.py index 2622f2c..d5ed075 100644 --- a/src/matadi/_templates.py +++ b/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)] diff --git a/src/matadi/math.py b/src/matadi/math.py index 7ef7b04..7e42ca5 100644 --- a/src/matadi/math.py +++ b/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 diff --git a/src/matadi/models/__init__.py b/src/matadi/models/__init__.py index f3d6d5a..2ade3df 100644 --- a/src/matadi/models/__init__.py +++ b/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", ] diff --git a/src/matadi/models/_templates.py b/src/matadi/models/_templates.py index d982050..89f9700 100644 --- a/src/matadi/models/_templates.py +++ b/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, + ) diff --git a/src/matadi/models/_viscoelasticity.py b/src/matadi/models/_viscoelasticity.py index 0576c23..b49f41a 100644 --- a/src/matadi/models/_viscoelasticity.py +++ b/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 + # <> + 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) diff --git a/tests/test_template-materials.py b/tests/test_template-materials.py index 213af5c..ed0ef15 100644 --- a/tests/test_template-materials.py +++ b/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)