#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# ---------------------------------------------------------------------
#
# _____ _ _ _ _____ _____
# | ____|__| | ___| |_ _____(_)___ ___| ___| ____|
# | _| / _` |/ _ \ \ \ /\ / / _ \ / __/ __| |_ | _|
# | |__| (_| | __/ |\ V V / __/ \__ \__ \ _| | |___
# |_____\__,_|\___|_| \_/\_/ \___|_|___/___/_| |_____|
#
#
# Unit of Strength of Materials and Structural Analysis
# University of Innsbruck,
# 2017 - today
#
# Daniel Reitmair daniel.reitmair@uibk.ac.at
#
# This file is part of EdelweissFE.
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# The full text of the license can be found in the file LICENSE.md at
# the top level directory of EdelweissFE.
# ---------------------------------------------------------------------
import numpy as np
import numpy.linalg as lin
from edelweissfe.materials.base.basehyperelasticmaterial import BaseHyperElasticMaterial
from edelweissfe.utils.exceptions import CutbackRequest
from edelweissfe.utils.plasticityutils import (
deviatoric,
deviatoricNorm,
dTensorExp_dA,
dTensorLog_dA,
tensorExp,
tensorLogarithmEig,
)
from edelweissfe.utils.voigtnotation import (
doVoigtMatrix,
doVoigtStrain,
doVoigtStress,
undoVoigtMatrix,
undoVoigtStrain,
)
# define variables and functions
maxIter = 20 # max inner Newton cycles
tol = 1e-12 # tolerance
# deviatoric multiplicator
IDev = np.array(
[
[2 / 3, -1 / 3, -1 / 3, 0, 0, 0],
[-1 / 3, 2 / 3, -1 / 3, 0, 0, 0],
[-1 / 3, -1 / 3, 2 / 3, 0, 0, 0],
[0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 1],
]
)
IDev4D = undoVoigtMatrix(3, IDev)
Ihs = np.diag([1, 1, 1, 0.5, 0.5, 0.5])
Ids = np.diag([1, 1, 1, 2, 2, 2])
Ieye = np.eye(3)
Ieye4D = undoVoigtMatrix(3, np.eye(6))
dR_dETrial = np.vstack([-np.eye(6), [np.zeros([6])]])
def _elasticWa(mu, K, Btr, dF, invFold, Bold):
"""Energy, Kirchhoff stress and elastic tangent for formulation Wa."""
J = np.sqrt(lin.det(Btr))
I1 = np.trace(Btr)
invB = lin.inv(Btr)
lambdaBar = (K - 2 / 3 * mu) * (J**2 - J) - mu
muBar = (K / 2 - mu / 3) * (2 * J**2 - J)
T = mu * Btr + lambdaBar * np.eye(3)
dB_dF = np.einsum("ik,np,jp,ln->ijkl", Ieye, Bold, dF, invFold) + np.einsum(
"io,on,jk,ln->ijkl", dF, Bold, Ieye, invFold
)
CTauF = mu * dB_dF + muBar * np.einsum("ij,mn,mnkl->ijkl", np.eye(3), invB.T, dB_dF)
energy = mu / 2 * (I1 - 3) + (K / 2 - mu / 3) * (J - 1) ** 2 - mu * np.log(J)
return (energy, T, CTauF)
def _kirchhoffStressWa(mu, K, e):
"""Kirchhoff stress for formulation Wa."""
B, _ = tensorExp(2 * e)
J = np.sqrt(lin.det(B))
muBar = (K - 2 / 3 * mu) * (J**2 - J) - mu
T = mu * B + muBar * Ieye
return T
def _dKirchhoff_dEWa(mu, K, e):
"""Derivative of the Kirchhoff stress w.r.t. the spatial logarithmic Hencky strain for formulation Wa."""
B, n = tensorExp(2 * e)
J = np.sqrt(lin.det(B))
dExpE = np.einsum("ijmn,mnkl->ijkl", dTensorExp_dA(2 * e, n), Ieye4D)
return 2 * mu * dExpE + (K - 2 / 3 * mu) * (2 * J**2 - J) * np.einsum("ij,kl->ijkl", Ieye, Ieye)
def _elasticWb(mu, K, Btr, dF, invFold, Bold):
"""Energy, Kirchhoff stress and elastic tangent for formulation Wb."""
J = np.sqrt(lin.det(Btr))
I1 = np.trace(Btr)
invB = lin.inv(Btr)
lambdaHat = K / 4 * (J**2 + 1 / J**2) + mu * I1 / (9 * J ** (2 / 3))
muBar = mu / (3 * J ** (2 / 3))
lambdaBar = K / 4 * (J**2 - 1 / J**2) - muBar * I1
T = 3 * muBar * Btr + lambdaBar * np.eye(3)
dB_dF = np.einsum("ik,np,jp,ln->ijkl", Ieye, Bold, dF, invFold) + np.einsum(
"io,on,jk,ln->ijkl", dF, Bold, Ieye, invFold
)
CTauF = muBar * (
-np.einsum("ij,mn,mnkl->ijkl", Btr, invB.T, dB_dF)
+ 3 * dB_dF
- np.einsum("ij,mn,mnkl->ijkl", Ieye, Ieye, dB_dF)
) + lambdaHat * np.einsum("ij,mn,mnkl->ijkl", np.eye(3), invB.T, dB_dF)
energy = mu / 2 * (I1 / J ** (2 / 3) - 3) + K / 8 * (J**2 + 1 / J**2 - 2)
return (energy, T, CTauF)
def _kirchhoffStressWb(mu, K, e):
"""Kirchhoff stress for formulation Wb."""
B, _ = tensorExp(2 * e)
I1 = B.trace()
J = np.sqrt(lin.det(B))
muBar = K / 4 * (J**2 - 1 / J**2) - mu * I1 / (3 * J ** (2 / 3))
T = mu / J ** (2 / 3) * B + muBar * Ieye
return T
def _dKirchhoff_dEWb(mu, K, e):
"""Derivative of the Kirchhoff stress w.r.t. the spatial logarithmic Hencky strain for formulation Wb."""
B, n = tensorExp(2 * e)
J = np.sqrt(lin.det(B))
dExpE = np.einsum("ijmn,mnkl->ijkl", dTensorExp_dA(2 * e, n), Ieye4D)
I1 = B.trace()
p = 2 * mu / J ** (2 / 3)
lambdaBar = K / 2 * (J**2 + 1 / J**2) + p * I1 / 9
dT_de = (
p * dExpE
- p / 3 * np.einsum("ij,kl->ijkl", B, Ieye)
+ lambdaBar * np.einsum("ij,kl->ijkl", Ieye, Ieye)
- p / 3 * np.einsum("ij,mn,mnkl->ijkl", Ieye, Ieye, dExpE)
)
return dT_de
def _elasticWc(mu, K, Btr, dF, invFold, Bold):
"""Energy, Kirchhoff stress and elastic tangent for formulation Wc."""
J = np.sqrt(lin.det(Btr))
I1 = Btr.trace()
invB = lin.inv(Btr)
muBar = mu * J ** (2 / 3 - K / mu)
lambdaBar = (K / mu - 2 / 3) * muBar / 2
T = mu * Btr - muBar * np.eye(3)
dB_dF = np.einsum("ik,np,jp,ln->ijkl", Ieye, Bold, dF, invFold) + np.einsum(
"io,on,jk,ln->ijkl", dF, Bold, Ieye, invFold
)
CTauF = mu * dB_dF + lambdaBar * np.einsum("ij,mn,mnkl->ijkl", Ieye, invB.T, dB_dF)
energy = mu / 2 * (I1 - 3) + 3 * mu**2 / (2 * K - 2 * mu) * (J ** (2 / 3 - K / mu) - 1)
return (energy, T, CTauF)
def _kirchhoffStressWc(mu, K, e):
"""Kirchhoff stress for formulation Wc."""
B, _ = tensorExp(2 * e)
muBar = np.exp((2 / 3 - K / mu) * e.trace())
T = mu * (B - muBar * Ieye)
return T
def _dKirchhoff_dEWc(mu, K, e):
"""Derivative of the Kirchhoff stress w.r.t. the spatial logarithmic Hencky strain for formulation Wc."""
B, n = tensorExp(2 * e)
J = np.sqrt(lin.det(B))
dExpE = np.einsum("ijmn,mnkl->ijkl", dTensorExp_dA(2 * e, n), Ieye4D)
return 2 * mu * dExpE + (K - 2 / 3 * mu) * J ** (2 / 3 - K / mu) * np.einsum("ij,kl->ijkl", Ieye, Ieye)
_FORMULATIONS = {
1: (_elasticWa, _kirchhoffStressWa, _dKirchhoff_dEWa),
2: (_elasticWb, _kirchhoffStressWb, _dKirchhoff_dEWb),
3: (_elasticWc, _kirchhoffStressWc, _dKirchhoff_dEWc),
}
[docs]
class NeoHookeanPlasticMaterial(BaseHyperElasticMaterial):
"""Shared J2-plasticity return-mapping implementation for the compressible neo-Hookean formulations
Wa/Wb/Wc according to [1], combined with isotropic hardening.
[1] Pence, T. J., & Gou, K. (2015). On compressible versions of the incompressible neo-Hookean material.
Mathematics and Mechanics of Solids, 20(2), 157–182. https://doi.org/10.1177/1081286514544258
Parameters
----------
materialProperties
The numpy array containing the material properties for the requested material, in order:
formulation selector (1=Wa, 2=Wb, 3=Wc), mu, K, yieldStress, HLin, deltaYieldStress, delta,
and optionally the density."""
[docs]
def getNumberOfRequiredStateVars(self) -> int:
"""Returns number of needed material state Variables per integration point in the material.
Returns
-------
int
Number of needed material state Vars."""
return 20
def __init__(self, materialProperties: np.ndarray):
self._materialProperties = materialProperties
formulation = int(materialProperties[0])
if formulation not in _FORMULATIONS:
raise Exception(
"Unknown Neo-Hookean Pence-Gou formulation selector {:}; must be 1 (Wa), 2 (Wb), or 3 (Wc).".format(
formulation
)
)
self._elasticImpl, self._kirchhoffStressImpl, self._dKirchhoff_dEImpl = _FORMULATIONS[formulation]
# elasticity parameters
self._mu = materialProperties[1]
self._K = materialProperties[2]
# plasticity parameters
self.yieldStress = materialProperties[3]
self.HLin = materialProperties[4]
self.deltaYieldStress = materialProperties[5]
self.delta = materialProperties[6]
# isotropic hardening
self._fy = (
lambda kappa_: self.yieldStress
+ self.HLin * kappa_
+ self.deltaYieldStress * (1.0 - np.exp(-self.delta * kappa_))
)
# yield function
self._f = lambda devNorm, kappa_: devNorm - np.sqrt(2 / 3) * self._fy(kappa_)
# derivative of fy dKappa
self._dfy_ddKappa = lambda kappa_: self.HLin + self.deltaYieldStress * self.delta * np.exp(-self.delta * kappa_)
if len(materialProperties) > 7:
self._density = materialProperties[7]
def _energyAndStressAndElasticTangent(self, Btr, dF, invFold, Bold):
return self._elasticImpl(self._mu, self._K, Btr, dF, invFold, Bold)
def _kirchhoffStress(self, e):
return self._kirchhoffStressImpl(self._mu, self._K, e)
def _dKirchhoff_dE(self, e):
return self._dKirchhoff_dEImpl(self._mu, self._K, e)
[docs]
def assignCurrentStateVars(self, currentStateVars: np.ndarray):
"""Assign new current state vars.
Parameters
----------
currentStateVars
Array containing the material state vars."""
self._energy = np.reshape(currentStateVars[0:1], (1))
self._kappaOld = np.reshape(currentStateVars[1:2], (1))
self._Fold = np.reshape(currentStateVars[2:11], (3, 3))
self._ee = np.reshape(currentStateVars[11:], (3, 3))
if np.all(self._Fold == 0):
self._Fold[:] = np.eye(3)
[docs]
def computeKirchhoff(
self,
stress: np.ndarray,
dStress_dDeformationGradient: np.ndarray,
deformationGradient: np.ndarray,
time: float,
dTime: float,
):
"""Computes the stresses for a 3D material/2D material with plane strain.
Parameters
----------
stress
Vector containing the stresses.
dStress_dDeformationGradient
Matrix containing dStress/dStrain.
deformationGradient
The deformation gradient at time step t.
time
Array of step time and total time.
dTime
Current time step size."""
F = deformationGradient
invFold = lin.inv(self._Fold)
dF = F @ invFold
# update to new values
Bold, _ = tensorExp(2 * self._ee)
Btrial = dF @ Bold @ dF.T
e_eTrial = 1 / 2 * tensorLogarithmEig(Btrial)
self._energy[0], trialKirchhoff, CTauF = self._energyAndStressAndElasticTangent(Btrial, dF, invFold, Bold)
# handle zero strain
if np.all(np.abs(dF - Ieye) < 1e-14):
dStress_dDeformationGradient[:] = CTauF
return
# calculate deviatoric norm of Kirchhoff stress
devNormKirchhoff = deviatoricNorm(deviatoric(trialKirchhoff))
if self._f(devNormKirchhoff, self._kappaOld[0]) > 0.0: # plastic step
# small strain algorithm to solve for kappa and epsilon
def residuals(X):
"""Compute the residual and its differential for the plastic solving process.
Parameters
----------
X
The current solution vector.
Returns
-------
np.ndarray
The residual vector.
np.ndarray
The derivative of the residual vector w.r.t. X."""
e_e = undoVoigtStrain(3, X[:6])
dKappa = X[6]
# compute Kirchhoff stress and Tensor exponential row number
kirchhoff = self._kirchhoffStress(e_e)
devKirchhoff = deviatoric(kirchhoff)
devNormKirchhoff = deviatoricNorm(devKirchhoff)
n = devKirchhoff / devNormKirchhoff # projection direction
# compute the residuals
Reps = doVoigtStrain(3, e_e - e_eTrial + np.sqrt(3 / 2) * dKappa * n)
Rk = devNormKirchhoff - np.sqrt(2 / 3) * self._fy(self._kappaOld[0] + dKappa)
# compute the residual derivations
dT_de = doVoigtMatrix(3, np.einsum("ijkl,klmn->ijmn", IDev4D, self._dKirchhoff_dE(e_e))) @ Ihs
nv = doVoigtStress(3, n) # make Voigt
dReps_dEe = np.eye(6) + np.sqrt(3 / 2) * dKappa / devNormKirchhoff * (
dT_de - Ids @ np.outer(nv, nv) @ dT_de
)
dReps_dKappa = Ids @ np.vstack(np.sqrt(3 / 2) * nv)
dRk_dKappa = -np.sqrt(2 / 3) * self._dfy_ddKappa(self._kappaOld[0] + dKappa)
# combine derivations into one matrix
dR_dX = np.vstack([np.hstack([dReps_dEe, dReps_dKappa]), np.hstack([nv @ dT_de, [dRk_dKappa]])])
return (np.hstack([Reps, [Rk]]), dR_dX)
# start values
X = np.hstack([doVoigtStrain(3, e_eTrial), 0])
counter = 0
R, dR_dX = residuals(X)
relTol = tol * self._fy(self._kappaOld[0])
while lin.norm(R) > relTol: # check if residuals fulfill tolerances
if counter == maxIter:
raise CutbackRequest("J2 plasticity solving failed.", 0.5)
# compute residual and its derivation
R, dR_dX = residuals(X)
X -= lin.solve(dR_dX, R)
counter += 1
# update kappa
self._kappaOld[:] += X[6]
# calculation of the elastoplastic consistent tangent
def elastoPlasticKirchhoffTangent(e):
"""Compute the elastoplastic Kirchhoff tangent dKirchhoff/dEtrial.
Parameters
----------
e
The elastic spatial Hencky strain.
Returns
-------
np.ndarray
The elastoplastic tangent dKirchhoff/dEtrial."""
dX_dEtrial = -lin.solve(dR_dX, dR_dETrial)
dT_dE = doVoigtMatrix(3, self._dKirchhoff_dE(e)) @ Ihs
return undoVoigtMatrix(3, dT_dE @ dX_dEtrial[:6, :6])
# update to back projected state
self._ee[:] = undoVoigtStrain(3, X[:6])
dT_dE = elastoPlasticKirchhoffTangent(self._ee)
# calculate back projected kirchhoff stress and convert to Voigt
T = self._kirchhoffStress(self._ee)
stress[:] = doVoigtStress(3, T)
# compute elastoplastic consistent tangent
dEe_dB = 1 / 2 * dTensorLog_dA(Btrial)
dB_dF = np.einsum("ik,np,jp,ln->ijkl", Ieye, Bold, dF, invFold) + np.einsum(
"io,on,jk,ln->ijkl", dF, Bold, Ieye, invFold
)
dEe_dF = np.einsum("ijmn,mnkl->ijkl", dEe_dB, dB_dF)
dStress_dDeformationGradient[:] = np.einsum("ijop,opkl->ijkl", dT_dE, dEe_dF)
else: # elastic step
stress[:] = doVoigtStress(3, trialKirchhoff)
dStress_dDeformationGradient[:] = CTauF
self._ee[:] = e_eTrial
self._Fold[:] = F
[docs]
def getResult(self, result: str) -> float:
"""Get the result, as a persistent view which is continiously
updated by the material.
Parameters
----------
result
The name of the result.
Returns
-------
float
The result.
"""
if result == "energy":
return self._energy
elif result == "kappa":
return self._kappaOld
else:
raise Exception("This result doesn't exist for the current material.")