-
Notifications
You must be signed in to change notification settings - Fork 13
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* Add `morph()` for JAX-based `Material` * remove Python 3.9 support (EOL) * Update index.rst * update to Python 3.9+ (docs) * Update _morph.py * Update _morph.py * Update _material.py * Create .gitkeep * Update CHANGELOG.md * Delete .gitkeep * Create __init__.py
- Loading branch information
Showing
13 changed files
with
283 additions
and
7 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -14,7 +14,7 @@ authors = [ | |
| {email = "[email protected]"}, | ||
|
|
||
| ] | ||
| requires-python = ">=3.8" | ||
| requires-python = ">=3.9" | ||
| license = {file = "LICENSE"} | ||
| keywords = [ | ||
| "computational-mechanics", | ||
|
|
@@ -38,7 +38,6 @@ classifiers = [ | |
| "Operating System :: OS Independent", | ||
| "Programming Language :: Python", | ||
| "Programming Language :: Python :: 3", | ||
| "Programming Language :: Python :: 3.8", | ||
| "Programming Language :: Python :: 3.9", | ||
| "Programming Language :: Python :: 3.10", | ||
| "Programming Language :: Python :: 3.11", | ||
|
|
||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,5 +1,6 @@ | ||
| from . import models | ||
| from ._hyperelastic import Hyperelastic | ||
| from ._material import Material | ||
| from ._tools import vmap | ||
|
|
||
| __all__ = ["Hyperelastic", "Material", "vmap"] | ||
| __all__ = ["Hyperelastic", "Material", "models", "vmap"] |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,3 @@ | ||
| from . import hyperelastic, lagrange | ||
|
|
||
| __all__ = ["hyperelastic", "lagrange"] |
Empty file.
8 changes: 8 additions & 0 deletions
8
src/felupe/constitution/autodiff/jax/models/lagrange/__init__.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,8 @@ | ||
| from ._morph import morph | ||
|
|
||
| __all__ = [ | ||
| "morph", | ||
| ] | ||
|
|
||
| # default (stable) material parameters | ||
| morph.kwargs = dict(p=[0, 0, 0, 0, 0, 1, 0, 0]) |
228 changes: 228 additions & 0 deletions
228
src/felupe/constitution/autodiff/jax/models/lagrange/_morph.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,228 @@ | ||
| # -*- 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 <http://www.gnu.org/licenses/>. | ||
| """ | ||
|
|
||
|
|
||
| def morph(F, statevars, p): | ||
| r"""Second Piola-Kirchhoff stress tensor of the | ||
| `MORPH <https://doi.org/10.1016/s0749-6419(02)00091-8>`_ model formulation [1]_. | ||
| Parameters | ||
| ---------- | ||
| C : tensortrax.Tensor | ||
| Right Cauchy-Green deformation tensor. | ||
| statevars : array | ||
| Vector of stacked state variables (CTS, C, SA). | ||
| p : list of float | ||
| A list which contains the 8 material parameters. | ||
| Notes | ||
| ----- | ||
| The MORPH material model is implemented as a second Piola-Kirchhoff stress-based | ||
| formulation with automatic differentiation. The Tresca invariant of the distortional | ||
| part of the right Cauchy-Green deformation tensor is used as internal state | ||
| variable, see Eq. :eq:`morph-state`. | ||
| .. warning:: | ||
| While the `MORPH <https://doi.org/10.1016/s0749-6419(02)00091-8>`_-material | ||
| formulation captures the Mullins effect and quasi-static hysteresis effects of | ||
| rubber mixtures very nicely, it has been observed to be unstable for medium- to | ||
| highly-distorted states of deformation. | ||
| .. math:: | ||
| :label: morph-state | ||
| \boldsymbol{C} &= \boldsymbol{F}^T \boldsymbol{F} | ||
| I_3 &= \det (\boldsymbol{C}) | ||
| \hat{\boldsymbol{C}} &= I_3^{-1/3} \boldsymbol{C} | ||
| \hat{\lambda}^2_\alpha &= \text{eigvals}(\hat{\boldsymbol{C}}) | ||
| \hat{C}_T &= \max \left( \hat{\lambda}^2_\alpha - \hat{\lambda}^2_\beta \right) | ||
| \hat{C}_T^S &= \max \left( \hat{C}_T, \hat{C}_{T,n}^S \right) | ||
| A sigmoid-function is used inside the deformation-dependent variables | ||
| :math:`\alpha`, :math:`\beta` and :math:`\gamma`, see Eq. :eq:`morph-sigmoid`. | ||
| .. math:: | ||
| :label: morph-sigmoid | ||
| f(x) &= \frac{1}{\sqrt{1 + x^2}} | ||
| \alpha &= p_1 + p_2 \ f(p_3\ C_T^S) | ||
| \beta &= p_4\ f(p_3\ C_T^S) | ||
| \gamma &= p_5\ C_T^S\ \left( 1 - f\left(\frac{C_T^S}{p_6}\right) \right) | ||
| The rate of deformation is described by the Lagrangian tensor and its Tresca- | ||
| invariant, see Eq. :eq:`morph-rate-of-deformation`. | ||
| .. note:: | ||
| It is important to evaluate the incremental right Cauchy-Green tensor by the | ||
| difference of the final and the previous state of deformation, not by its | ||
| variation with respect to the deformation gradient tensor. | ||
| .. math:: | ||
| :label: morph-rate-of-deformation | ||
| \hat{\boldsymbol{L}} &= \text{sym}\left( | ||
| \text{dev}(\boldsymbol{C}^{-1} \Delta\boldsymbol{C}) | ||
| \right) \hat{\boldsymbol{C}} | ||
| \lambda_{\hat{\boldsymbol{L}}, \alpha} &= \text{eigvals}(\hat{\boldsymbol{L}}) | ||
| \hat{L}_T &= \max \left( | ||
| \lambda_{\hat{\boldsymbol{L}}, \alpha}-\lambda_{\hat{\boldsymbol{L}}, \beta} | ||
| \right) | ||
| \Delta\boldsymbol{C} &= \boldsymbol{C} - \boldsymbol{C}_n | ||
| The additional stresses evolve between the limiting stresses, see Eq. | ||
| :eq:`morph-stresses`. The additional deviatoric-enforcement terms [1]_ are neglected | ||
| in this implementation. | ||
| .. math:: | ||
| :label: morph-stresses | ||
| \boldsymbol{S}_L &= \left( | ||
| \gamma \exp \left(p_7 \frac{\hat{\boldsymbol{L}}}{\hat{L}_T} | ||
| \frac{\hat{C}_T}{\hat{C}_T^S} \right) + | ||
| p8 \frac{\hat{\boldsymbol{L}}}{\hat{L}_T} | ||
| \right) \boldsymbol{C}^{-1} | ||
| \boldsymbol{S}_A &= \frac{ | ||
| \boldsymbol{S}_{A,n} + \beta\ \hat{L}_T\ \boldsymbol{S}_L | ||
| }{1 + \beta\ \hat{L}_T} | ||
| \boldsymbol{S} &= 2 \alpha\ \text{dev}( \hat{\boldsymbol{C}} ) | ||
| \boldsymbol{C}^{-1}+\text{dev}\left(\boldsymbol{S}_A\ \boldsymbol{C}\right) | ||
| \boldsymbol{C}^{-1} | ||
| Examples | ||
| -------- | ||
| .. pyvista-plot:: | ||
| :context: | ||
| >>> import felupe as fem | ||
| >>> | ||
| >>> umat = fem.constitution.autodiff.jax.Material( | ||
| ... fem.constitution.autodiff.jax.models.lagrange.morph, | ||
| ... p=[0.039, 0.371, 0.174, 2.41, 0.0094, 6.84, 5.65, 0.244], | ||
| ... nstatevars=13, | ||
| ... ) | ||
| >>> ax = umat.plot( | ||
| ... incompressible=True, | ||
| ... ux=fem.math.linsteps( | ||
| ... # [1, 2, 1, 2.75, 1, 3.5, 1, 4.2, 1, 4.8, 1, 4.8, 1], | ||
| ... [1, 2.75, 1, 2.75], | ||
| ... num=20, | ||
| ... ), | ||
| ... ps=None, | ||
| ... bx=None, | ||
| ... ) | ||
| .. pyvista-plot:: | ||
| :include-source: False | ||
| :context: | ||
| :force_static: | ||
| >>> import pyvista as pv | ||
| >>> | ||
| >>> fig = ax.get_figure() | ||
| >>> chart = pv.ChartMPL(fig) | ||
| >>> chart.show() | ||
| References | ||
| ---------- | ||
| .. [1] D. Besdo and J. Ihlemann, "A phenomenological constitutive model for | ||
| rubberlike materials and its numerical applications", International Journal | ||
| of Plasticity, vol. 19, no. 7. Elsevier BV, pp. 1019–1036, Jul. 2003. doi: | ||
| `10.1016/s0749-6419(02)00091-8 <https://doi.org/10.1016/s0749-6419(02)00091-8>`_. | ||
| """ | ||
|
|
||
| from jax.numpy import ( | ||
| array, | ||
| concatenate, | ||
| diag, | ||
| eye, | ||
| maximum, | ||
| sqrt, | ||
| trace, | ||
| triu_indices, | ||
| ) | ||
| from jax.numpy.linalg import det, eigvalsh, inv | ||
| from jax.scipy.linalg import expm | ||
|
|
||
| # right Cauchy-Green deformation tensor | ||
| C = F.T @ F | ||
|
|
||
| # extract old state variables | ||
| CTSn = statevars[0] | ||
| from_triu = lambda C: C[array([[0, 1, 2], [1, 3, 4], [2, 4, 5]])] | ||
| Cn = from_triu(statevars[1:7]) | ||
| SAn = from_triu(statevars[7:13]) | ||
|
|
||
| # distortional part of right Cauchy-Green deformation tensor | ||
| I3 = det(C) | ||
| CG = C * I3 ** (-1 / 3) | ||
|
|
||
| # inverse of and incremental right Cauchy-Green deformation tensor | ||
| invC = inv(C) | ||
| dC = C - Cn | ||
|
|
||
| # eigenvalues of right Cauchy-Green deformation tensor (sorted in ascending order) | ||
| eigvalsh2 = lambda C: eigvalsh(C + diag(array([1e-4, -1e-4, 0]))) | ||
| λCG = eigvalsh2(CG) | ||
|
|
||
| # Tresca invariant of distortional part of right Cauchy-Green deformation tensor | ||
| CTG = λCG[-1] - λCG[0] | ||
|
|
||
| # maximum Tresca invariant in load history | ||
| CTS = maximum(CTG, CTSn) | ||
|
|
||
| def sigmoid(x): | ||
| "Algebraic sigmoid function." | ||
| return 1 / sqrt(1 + x**2) | ||
|
|
||
| # material parameters | ||
| α = p[0] + p[1] * sigmoid(p[2] * CTS) | ||
| β = p[3] * sigmoid(p[2] * CTS) | ||
| γ = p[4] * CTS * (1 - sigmoid(CTS / p[5])) | ||
|
|
||
| dev = lambda C: C - trace(C) / 3 * eye(3) | ||
| sym = lambda C: (C + C.T) / 2 | ||
|
|
||
| LG = sym(dev(invC @ dC)) @ CG | ||
| λLG = eigvalsh2(LG) | ||
| LTG = λLG[-1] - λLG[0] | ||
|
|
||
| # limiting stresses "L" and additional stresses "A" | ||
| SL = (γ * expm(p[6] * LG / LTG * CTG / CTS) + p[7] * LG / LTG) @ invC | ||
| SA = (SAn + β * LTG * SL) / (1 + β * LTG) | ||
|
|
||
| # second Piola-Kirchhoff stress tensor | ||
| S = 2 * α * dev(CG) @ invC + dev(SA @ C) @ invC | ||
|
|
||
| i, j = triu_indices(3) | ||
| to_triu = lambda C: C[i, j] | ||
| statevars_new = concatenate([array([CTS]), to_triu(C), to_triu(SA)]) | ||
|
|
||
| return F @ S, statevars_new |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters