From 98f42f3cb418042956495c2a380c3e24ce4d68e8 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Tue, 19 Nov 2024 21:34:17 +0100 Subject: [PATCH] Add `constitution.lame_converter_orthotropic()` (#904) --- docs/felupe/constitution/tools.rst | 3 + src/felupe/constitution/__init__.py | 2 + .../linear_elasticity/__init__.py | 3 +- .../linear_elasticity/_lame_converter.py | 71 +++++++++++++++++ .../_linear_elastic_orthotropic.py | 76 +++++-------------- .../_saint_venant_kirchhoff_orthotropic.py | 18 +++-- tests/test_constitution.py | 22 +++--- 7 files changed, 120 insertions(+), 75 deletions(-) diff --git a/docs/felupe/constitution/tools.rst b/docs/felupe/constitution/tools.rst index ed6f658f..7614342e 100644 --- a/docs/felupe/constitution/tools.rst +++ b/docs/felupe/constitution/tools.rst @@ -29,6 +29,7 @@ This page contains tools and helpers for constitutive material formulations. .. autosummary:: constitution.lame_converter + constitution.lame_converter_orthotropic **Detailed API Reference** @@ -68,3 +69,5 @@ This page contains tools and helpers for constitutive material formulations. :inherited-members: .. autofunction:: felupe.constitution.lame_converter + +.. autofunction:: felupe.constitution.lame_converter_orthotropic \ No newline at end of file diff --git a/src/felupe/constitution/__init__.py b/src/felupe/constitution/__init__.py index 13d0f113..1907f408 100644 --- a/src/felupe/constitution/__init__.py +++ b/src/felupe/constitution/__init__.py @@ -12,6 +12,7 @@ LinearElasticPlaneStress, LinearElasticTensorNotation, lame_converter, + lame_converter_orthotropic, ) from .poisson import Laplace from .small_strain import MaterialStrain @@ -32,6 +33,7 @@ "LinearElasticPlaneStress", "LinearElasticTensorNotation", "lame_converter", + "lame_converter_orthotropic", "OgdenRoxburgh", "LinearElasticPlasticIsotropicHardening", "Material", diff --git a/src/felupe/constitution/linear_elasticity/__init__.py b/src/felupe/constitution/linear_elasticity/__init__.py index ca1ce89f..3e0f8ef9 100644 --- a/src/felupe/constitution/linear_elasticity/__init__.py +++ b/src/felupe/constitution/linear_elasticity/__init__.py @@ -1,4 +1,4 @@ -from ._lame_converter import lame_converter +from ._lame_converter import lame_converter, lame_converter_orthotropic from ._linear_elastic import ( LinearElastic, LinearElasticPlaneStrain, @@ -10,6 +10,7 @@ __all__ = [ "lame_converter", + "lame_converter_orthotropic", "LinearElastic", "LinearElasticLargeStrain", "LinearElasticOrthotropic", diff --git a/src/felupe/constitution/linear_elasticity/_lame_converter.py b/src/felupe/constitution/linear_elasticity/_lame_converter.py index 611a84a5..3e688876 100644 --- a/src/felupe/constitution/linear_elasticity/_lame_converter.py +++ b/src/felupe/constitution/linear_elasticity/_lame_converter.py @@ -15,6 +15,7 @@ You should have received a copy of the GNU General Public License along with FElupe. If not, see . """ +from numpy.linalg import inv def lame_converter(E, nu): @@ -50,3 +51,73 @@ def lame_converter(E, nu): mu = E / (2 * (1 + nu)) return lmbda, mu + + +def lame_converter_orthotropic(E, nu, G): + r"""Convert elastic orthotropic material parameters to Lamé parameter. + + Parameters + ---------- + E : list of float + List of the three elastic moduli :math:`E_1, E_2, E_3`. + nu : list of float + List of three poisson ratios :math:`\nu_{12}, \nu_{23}, \nu_{31}`. + G : list of float + List of three shear moduli :math:`G_{12}, G_{23}, G_{31}`. + + Returns + ------- + lmbda : list of float + List of six (upper triangle) first Lamé parameters + :math:`\lambda_11, \lambda_12, \lambda_13, \lambda_22, \lambda_23, \lambda_33`. + mu : list of float + List of the three second Lamé parameters :math:`\mu_1,\mu_2, \mu3`. + + """ + + # unpack orthotropic elastic material parameters + E1, E2, E3 = E + ν12, ν23, ν31 = nu + G12, G23, G31 = G + + # orthotropic symmetry + ν21 = ν12 * E2 / E1 + ν32 = ν23 * E3 / E2 + ν13 = ν31 * E1 / E3 + + C = inv( + [ + [1 / E1, -ν21 / E2, -ν31 / E3, 0, 0, 0], + [-ν12 / E1, 1 / E2, -ν32 / E3, 0, 0, 0], + [-ν13 / E1, -ν23 / E2, 1 / E3, 0, 0, 0], + [0, 0, 0, 1 / G12, 0, 0], + [0, 0, 0, 0, 1 / G23, 0], + [0, 0, 0, 0, 0, 1 / G31], + ] + ) + + # take the components from this matrix + # + # [λ11 + 2 * μ1, λ12, λ13, 0, 0, 0] + # [λ12, λ22 + 2 * μ2, λ23, 0, 0, 0] + # [λ13, λ23, λ33 + 2 * μ3, 0, 0, 0] + # [0, 0, 0, (μ1 + μ2) / 2, 0, 0] + # [0, 0, 0, 0, (μ2 + μ3) / 2, 0] + # [0, 0, 0, 0, 0, (μ1 + μ3) / 2] + + λ12 = C[0, 1] + λ23 = C[1, 2] + λ13 = C[0, 2] + + μ1 = C[3, 3] - C[4, 4] + C[5, 5] + μ2 = C[3, 3] + C[4, 4] - C[5, 5] + μ3 = -C[3, 3] + C[4, 4] + C[5, 5] + + λ11 = C[0, 0] - 2 * μ1 + λ22 = C[1, 1] - 2 * μ2 + λ33 = C[2, 2] - 2 * μ3 + + lmbda = [λ11, λ12, λ13, λ22, λ23, λ33] + mu = [μ1, μ2, μ3] + + return lmbda, mu diff --git a/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py b/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py index b6ff19ec..62967864 100644 --- a/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py +++ b/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py @@ -27,24 +27,12 @@ class LinearElasticOrthotropic(ConstitutiveMaterial): Parameters ---------- - E1 : float - Young's modulus. - E2 : float - Young's modulus. - E3 : float - Young's modulus. - nu12 : float - Poisson ratio. - nu23 : float - Poisson ratio. - nu13 : float - Poisson ratio. - G12 : float - Shear modulus. - G23 : float - Shear modulus. - G13 : float - Shear modulus. + E : float + Young's modulus (E1, E2, E3). + nu : float + Poisson ratio (nu12, nu23, n31). + G : float + Shear modulus (G12, G23, G31). Examples -------- @@ -54,8 +42,8 @@ class LinearElasticOrthotropic(ConstitutiveMaterial): >>> import felupe as fem >>> >>> umat = fem.LinearElasticOrthotropic( - >>> E1=1, E2=1, E3=1, nu12=0.3, nu23=0.3, nu13=0.3, G12=0.4, G23=0.4, G13=0.4 - >>> ) + ... E=[1, 1, 1], nu=[0.3, 0.3, 0.3], G=[0.4, 0.4, 0.4] + ... ) >>> ax = umat.plot() .. pyvista-plot:: @@ -71,30 +59,12 @@ class LinearElasticOrthotropic(ConstitutiveMaterial): """ - def __init__(self, E1, E2, E3, nu12, nu23, nu13, G12, G23, G13): - self.E1 = E1 - self.E2 = E1 - self.E3 = E3 - - self.nu12 = nu12 - self.nu23 = nu23 - self.nu13 = nu13 - - self.G12 = G12 - self.G23 = G23 - self.G13 = G13 - - self.kwargs = { - "E1": self.E1, - "E2": self.E2, - "E3": self.E3, - "nu12": self.nu12, - "nu23": self.nu23, - "nu13": self.nu13, - "G12": self.G12, - "G23": self.G23, - "G13": self.G13, - } + def __init__(self, E, nu, G): + self.E = E + self.nu = nu + self.G = G + + self.kwargs = {"E": self.E, "nu": self.nu, "G": self.G} # aliases for gradient and hessian self.stress = self.gradient @@ -149,17 +119,9 @@ def hessian(self, x=None, shape=(1, 1), dtype=None): """ - E1 = self.E1 - E2 = self.E2 - E3 = self.E3 - - nu12 = self.nu12 - nu23 = self.nu23 - nu13 = self.nu13 - - G12 = self.G12 - G23 = self.G23 - G13 = self.G13 + E1, E2, E3 = self.E + nu12, nu23, nu31 = self.nu + G12, G23, G31 = self.G if x is not None: dtype = x[0].dtype @@ -168,7 +130,7 @@ def hessian(self, x=None, shape=(1, 1), dtype=None): nu21 = nu12 * E2 / E1 nu32 = nu23 * E3 / E2 - nu31 = nu13 * E3 / E1 + nu13 = nu31 * E1 / E3 J = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13) @@ -192,7 +154,7 @@ def hessian(self, x=None, shape=(1, 1), dtype=None): [2, 0, 2, 0], [0, 0, 2, 2], [2, 2, 0, 0], - ] = G13 + ] = G31 elast[ [1, 2, 1, 2], diff --git a/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py b/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py index 63a0774a..f79546ea 100644 --- a/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py +++ b/src/felupe/constitution/tensortrax/models/hyperelastic/_saint_venant_kirchhoff_orthotropic.py @@ -30,9 +30,10 @@ def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3): C : tensortrax.Tensor Right Cauchy-Green deformation tensor. mu : list of float - List of the three shear moduli. + List of the three second Lamé parameters :math:`\mu_1,\mu_2, \mu3`. lmbda : list of float - The six (upper triangle) orthotropic moduli. + List of six (upper triangle) first Lamé parameters + :math:`\lambda_11, \lambda_12, \lambda_13, \lambda_22, \lambda_23, \lambda_33`. r1 : list of float First normal vector of planes of symmetry. r2 : list of float @@ -65,16 +66,21 @@ def saint_venant_kirchhoff_orthotropic(C, mu, lmbda, r1, r2, r3): >>> import felupe.constitution.tensortrax as mat >>> import felupe as fem >>> - >>> r = fem.math.rotation_matrix(30, axis=2) + >>> r = fem.math.rotation_matrix(0, axis=2) + >>> lmbda, mu = fem.constitution.lame_converter_orthotropic( + ... E=[10, 10, 10], + ... nu=[0.3, 0.3, 0.3], + ... G=[1, 1, 1], + ... ) >>> umat = mat.Hyperelastic( ... mat.models.hyperelastic.saint_venant_kirchhoff_orthotropic, - ... mu=[1, 1, 1], - ... lmbda=[20, 20, 20, 20, 20, 20], + ... mu=mu, + ... lmbda=lmbda, ... r1=r[:, 0], ... r2=r[:, 1], ... r3=r[:, 2], ... ) - >>> ax = umat.plot(incompressible=False) + >>> ax = umat.plot(ux=fem.math.linsteps([1, 1.1], num=10), ps=None, bx=None) .. pyvista-plot:: :include-source: False diff --git a/tests/test_constitution.py b/tests/test_constitution.py index 7f021187..7c127d44 100644 --- a/tests/test_constitution.py +++ b/tests/test_constitution.py @@ -136,22 +136,22 @@ def test_linear(): def test_linear_orthotropic(): r, F = pre(sym=False, add_identity=True) + lmbda, mu = fem.constitution.lame_converter_orthotropic( + E=[6919, 271, 450], + nu=[0.388, 0.278, 0.375], + G=[262, 34, 354], + ) + for Material in [ (fem.constitution.LinearElasticOrthotropic, {}), ]: - LinearElastic, kwargs = Material + LinearElasticOrtho, kwargs = Material # doi.org/10.2478/ace-2018-0027 (pine wood) - le = LinearElastic( - E1=6919, - E2=271, - E3=450, - nu12=0.388, - nu23=0.278, - nu13=0.375, - G12=262, - G23=34, - G13=354, + le = LinearElasticOrtho( + E=[6919, 271, 450], + nu=[0.388, 0.278, 0.375], + G=[262, 34, 354], **kwargs, )