From 6f5d5424013be0ad69db0440ab9de2ff07ba1458 Mon Sep 17 00:00:00 2001 From: Sad-Abd Date: Fri, 10 May 2024 21:26:55 +0330 Subject: [PATCH 1/5] Added `unimodular` to calculate unimodular part of a tensor --- src/matadi/math.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/matadi/math.py b/src/matadi/math.py index 7ef7b04..7c0f180 100644 --- a/src/matadi/math.py +++ b/src/matadi/math.py @@ -257,3 +257,14 @@ 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 \ No newline at end of file From d992ae8b59a850a38d08e3966e654ce054cf185e Mon Sep 17 00:00:00 2001 From: Sad-Abd Date: Fri, 10 May 2024 21:30:53 +0330 Subject: [PATCH 2/5] Added `sqrtm` to compute matrix square root --- src/matadi/math.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/matadi/math.py b/src/matadi/math.py index 7c0f180..3489c53 100644 --- a/src/matadi/math.py +++ b/src/matadi/math.py @@ -267,4 +267,17 @@ def unimodular(T): This operation preserves the isochoric (volume-preserving) part of the tensor while removing the volumetric part. """ - return (det(T) ** (-1 / 3)) * T \ No newline at end of file + 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 \ No newline at end of file From f0c91a3c49e901649e1e57ed864917c44be21984 Mon Sep 17 00:00:00 2001 From: Sad-Abd Date: Fri, 10 May 2024 21:48:59 +0330 Subject: [PATCH 3/5] Added finite strain viscoelasticity with Mooney-Rivlin --- src/matadi/models/_viscoelasticity.py | 34 ++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/src/matadi/models/_viscoelasticity.py b/src/matadi/models/_viscoelasticity.py index 0576c23..ce54522 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, gradient, inv, trace, eye, unimodular, sqrtm def finite_strain_viscoelastic(x, mu, eta, dtime): @@ -16,3 +16,35 @@ 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) From 1faa9ccd25752a28119b45b30d45b461a6bfb028 Mon Sep 17 00:00:00 2001 From: Sad-Abd Date: Fri, 10 May 2024 21:54:28 +0330 Subject: [PATCH 4/5] Added `Viscoelastic_MR` to templates --- src/matadi/models/_templates.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/src/matadi/models/_templates.py b/src/matadi/models/_templates.py index d982050..83add21 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,26 @@ def fun(x, mu, eta, dtime): eta=eta, dtime=dtime, ) + +class Viscoelastic_MR(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, + ) From fac7b14bf1c1583db407c3c47121e91e9fcf3a33 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Sat, 11 May 2024 00:41:47 +0200 Subject: [PATCH 5/5] Finalize the viscoelastic Mooney-Rivlin model --- README.md | 6 ++++-- src/matadi/_templates.py | 10 ++++++---- src/matadi/math.py | 22 ++++++++++++---------- src/matadi/models/__init__.py | 10 ++++++++-- src/matadi/models/_templates.py | 19 ++++++++++--------- src/matadi/models/_viscoelasticity.py | 7 ++++--- tests/test_template-materials.py | 8 +++++++- 7 files changed, 51 insertions(+), 31 deletions(-) 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 3489c53..7e42ca5 100644 --- a/src/matadi/math.py +++ b/src/matadi/math.py @@ -258,18 +258,20 @@ def astensor(A, scale=1): else: raise ValueError("Unknown shape of input.") + def unimodular(T): - """ - Compute the unimodular part of a tensor. + """ + 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 - 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): +def sqrtm(C, eps=8e-5): """ Compute the matrix square root of a tensor C using eigendecomposition. """ @@ -280,4 +282,4 @@ def sqrtm(C, eps = 8e-5): 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 \ No newline at end of file + 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 83add21..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, finite_strain_viscoelastic_MR +from ._viscoelasticity import finite_strain_viscoelastic, finite_strain_viscoelastic_mr class NeoHookeOgdenRoxburgh(MaterialTensorGeneral): @@ -81,26 +81,27 @@ def fun(x, mu, eta, dtime): eta=eta, dtime=dtime, ) - -class Viscoelastic_MR(MaterialTensorGeneral): + + +class ViscoelasticMooneyRivlin(MaterialTensorGeneral): "Finite strain viscoelastic material formulation with Mooney-Rivlin hyperelasticity." def __init__( self, - c10 = 1, - c01 = 1, - eta = 1, + 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) + 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, + c10=c10, + c01=c01, eta=eta, dtime=dtime, ) diff --git a/src/matadi/models/_viscoelasticity.py b/src/matadi/models/_viscoelasticity.py index ce54522..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, eye, unimodular, sqrtm +from ..math import astensor, asvoigt, det, eye, gradient, inv, sqrtm, trace, unimodular def finite_strain_viscoelastic(x, mu, eta, dtime): @@ -18,7 +18,7 @@ def finite_strain_viscoelastic(x, mu, eta, dtime): return gradient(mu / 2 * (I1 - 3), F), asvoigt(Ci) -def finite_strain_viscoelastic_MR(x, c10, c01, eta, dtime): +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 @@ -30,7 +30,8 @@ def finite_strain_viscoelastic_MR(x, c10, c01, eta, dtime): # Right cauchy-green deformation tensor C = F.T @ F - # Based on <
> + # Based on + # <
> A = ( unimodular(sqrtm(inv(C))) @ (astensor(Cin) + (dtime / eta) * c10 * unimodular(C)) 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)