From a567a817e5027a1b96a10c2420a6664f45b84506 Mon Sep 17 00:00:00 2001 From: Andreas Dutzler Date: Thu, 29 Aug 2024 23:22:48 +0200 Subject: [PATCH] Add `LinearElasticOrthotropic` (#840) * Add `LinearElasticOrthotropic` * Update _linear_elastic_orthotropic.py * Update test_constitution.py * Update CHANGELOG.md * Update test_constitution.py * Update _linear_elastic_orthotropic.py * Add `LinearElasticOrthotropic` to the tests --- CHANGELOG.md | 1 + docs/felupe/constitution/core.rst | 6 + src/felupe/__init__.py | 2 + src/felupe/constitution/__init__.py | 2 + .../linear_elasticity/__init__.py | 2 + .../_linear_elastic_orthotropic.py | 204 ++++++++++++++++++ tests/test_constitution.py | 34 +++ 7 files changed, 251 insertions(+) create mode 100644 src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 6ec87e2e..12cd0983 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ All notable changes to this project will be documented in this file. The format - Add an optional `mechanics.Assemble(..., multiplier=None)` argument which is used in external items like `SolidBodyForce`, `SolidBodyGravity` and `PointLoad` and is applied in `newtonrhapson(items, ...)`. - Add a new submodule `view` which contains the `View...` classes like `ViewSolid` or `ViewField`, previously located at `tools._plot`. - Add the `grad`- and `hess`-arguments to the `reload()`- and `copy()`-methods of a `Region`, i.e. `Region.reload(grad=None, hess=None)`. +- Add `LinearElasticOrthotropic`. ### Changed - Change the internal initialization of `field = Field(region, values=1, dtype=None)` values from `field.values = np.ones(shape) * values` to `field = np.full(shape, fill_value=values, dtype=dtype)`. This enforces `field = Field(region, values=1)` to return the gradient array with data-type `int` which was of type `float` before. diff --git a/docs/felupe/constitution/core.rst b/docs/felupe/constitution/core.rst index 23bb1d86..b28a5398 100644 --- a/docs/felupe/constitution/core.rst +++ b/docs/felupe/constitution/core.rst @@ -22,6 +22,7 @@ This page contains the core (hard-coded) constitutive material model formulation LinearElasticPlaneStress constitution.LinearElasticTensorNotation LinearElasticLargeStrain + LinearElasticOrthotropic **Plasticity** @@ -93,6 +94,11 @@ This page contains the core (hard-coded) constitutive material model formulation .. autofunction:: felupe.linear_elastic_plastic_isotropic_hardening +.. autoclass:: felupe.LinearElasticOrthotropic + :members: + :undoc-members: + :inherited-members: + .. autoclass:: felupe.MaterialStrain :members: :undoc-members: diff --git a/src/felupe/__init__.py b/src/felupe/__init__.py index 3a0edec9..82629352 100644 --- a/src/felupe/__init__.py +++ b/src/felupe/__init__.py @@ -23,6 +23,7 @@ Laplace, LinearElastic, LinearElasticLargeStrain, + LinearElasticOrthotropic, LinearElasticPlaneStrain, LinearElasticPlaneStress, LinearElasticPlasticIsotropicHardening, @@ -171,6 +172,7 @@ "Laplace", "LinearElastic", "LinearElasticLargeStrain", + "LinearElasticOrthotropic", "LinearElasticPlaneStrain", "LinearElasticPlaneStress", "LinearElasticPlasticIsotropicHardening", diff --git a/src/felupe/constitution/__init__.py b/src/felupe/constitution/__init__.py index 673b0981..1fa185f3 100644 --- a/src/felupe/constitution/__init__.py +++ b/src/felupe/constitution/__init__.py @@ -33,6 +33,7 @@ from .linear_elasticity import ( LinearElastic, LinearElasticLargeStrain, + LinearElasticOrthotropic, LinearElasticPlaneStrain, LinearElasticPlaneStress, LinearElasticTensorNotation, @@ -70,6 +71,7 @@ "Laplace", "LinearElastic", "LinearElasticLargeStrain", + "LinearElasticOrthotropic", "LinearElasticPlaneStrain", "LinearElasticPlaneStress", "LinearElasticTensorNotation", diff --git a/src/felupe/constitution/linear_elasticity/__init__.py b/src/felupe/constitution/linear_elasticity/__init__.py index e7777677..ca1ce89f 100644 --- a/src/felupe/constitution/linear_elasticity/__init__.py +++ b/src/felupe/constitution/linear_elasticity/__init__.py @@ -6,11 +6,13 @@ LinearElasticTensorNotation, ) from ._linear_elastic_large_strain import LinearElasticLargeStrain +from ._linear_elastic_orthotropic import LinearElasticOrthotropic __all__ = [ "lame_converter", "LinearElastic", "LinearElasticLargeStrain", + "LinearElasticOrthotropic", "LinearElasticPlaneStrain", "LinearElasticPlaneStress", "LinearElasticTensorNotation", diff --git a/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py b/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py new file mode 100644 index 00000000..b6ff19ec --- /dev/null +++ b/src/felupe/constitution/linear_elasticity/_linear_elastic_orthotropic.py @@ -0,0 +1,204 @@ +# -*- coding: utf-8 -*- +""" +This file is part of FElupe. + +FElupe is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +FElupe is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with FElupe. If not, see . +""" + +import numpy as np + +from ...math import ddot, identity, transpose +from .._base import ConstitutiveMaterial + + +class LinearElasticOrthotropic(ConstitutiveMaterial): + r"""Orthotropic linear-elastic material formulation. + + 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. + + Examples + -------- + .. pyvista-plot:: + :context: + + >>> 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 + >>> ) + >>> ax = umat.plot() + + .. pyvista-plot:: + :include-source: False + :context: + :force_static: + + >>> import pyvista as pv + >>> + >>> fig = ax.get_figure() + >>> chart = pv.ChartMPL(fig) + >>> chart.show() + + """ + + 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, + } + + # aliases for gradient and hessian + 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)] + + def gradient(self, x): + r"""Evaluate the stress tensor (as a function of the deformation gradient). + + Parameters + ---------- + x : list of ndarray + List with Deformation gradient :math:`\boldsymbol{F}` (3x3) as first item. + + Returns + ------- + ndarray of shape (3, 3, ...) + Stress tensor + + """ + + F, statevars = x[0], x[-1] + + # convert the deformation gradient to strain + H = F - identity(F) + strain = (H + transpose(H)) / 2 + + # init stress + elast = self.hessian(x=x)[0] + + return [ddot(elast, strain, mode=(4, 2)), statevars] + + def hessian(self, x=None, shape=(1, 1), dtype=None): + r"""Evaluate the elasticity tensor. The Deformation gradient is only + used for the shape of the trailing axes. + + Parameters + ---------- + x : list of ndarray, optional + List with Deformation gradient :math:`\boldsymbol{F}` (3x3) as first item + (default is None). + shape : tuple of int, optional + Tuple with shape of the trailing axes (default is (1, 1)). + + Returns + ------- + ndarray of shape (3, 3, 3, 3, ...) + elasticity tensor + + """ + + 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 + + if x is not None: + dtype = x[0].dtype + + elast = np.zeros((3, 3, 3, 3, *shape), dtype=dtype) + + nu21 = nu12 * E2 / E1 + nu32 = nu23 * E3 / E2 + nu31 = nu13 * E3 / E1 + + J = 1 / (1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13) + + elast[0, 0, 0, 0] = E1 * (1 - nu23 * nu32) * J + elast[1, 1, 1, 1] = E2 * (1 - nu13 * nu31) * J + elast[2, 2, 2, 2] = E3 * (1 - nu12 * nu21) * J + + elast[0, 0, 1, 1] = elast[1, 1, 0, 0] = E1 * (nu21 + nu31 * nu23) * J + elast[0, 0, 2, 2] = elast[2, 2, 0, 0] = E1 * (nu31 + nu21 * nu32) * J + elast[1, 1, 2, 2] = elast[2, 2, 1, 1] = E2 * (nu32 + nu12 * nu31) * J + + elast[ + [0, 1, 0, 1], + [1, 0, 1, 0], + [0, 0, 1, 1], + [1, 1, 0, 0], + ] = G12 + + elast[ + [0, 2, 0, 2], + [2, 0, 2, 0], + [0, 0, 2, 2], + [2, 2, 0, 0], + ] = G13 + + elast[ + [1, 2, 1, 2], + [2, 1, 2, 1], + [1, 1, 2, 2], + [2, 2, 1, 1], + ] = G23 + + return [elast] diff --git a/tests/test_constitution.py b/tests/test_constitution.py index 15e2fe95..898a1289 100644 --- a/tests/test_constitution.py +++ b/tests/test_constitution.py @@ -150,6 +150,39 @@ def test_linear(): assert np.allclose(*check_dsde) +def test_linear_orthotropic(): + r, F = pre(sym=False, add_identity=True) + + for Material in [ + (fem.constitution.LinearElasticOrthotropic, {}), + ]: + LinearElastic, 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, + **kwargs, + ) + + stress = le.gradient(F)[:-1] + dsde = le.hessian(F) + + assert le.elasticity()[0].shape[-2:] == (1, 1) + + assert stress[0].shape == (3, 3, *F[0].shape[-2:]) + assert dsde[0].shape == (3, 3, 3, 3, 1, 1) + + assert np.allclose(stress, 0) + + def test_linear_planestress(): r, F = pre(sym=False, add_identity=True) F = [F[0][:2][:, :2]] @@ -708,6 +741,7 @@ def test_laplace(): if __name__ == "__main__": test_nh() test_linear() + test_linear_orthotropic() test_linear_planestress() test_linear_planestrain() test_kinematics()