Skip to content

Commit

Permalink
Plane Strain models from 3D models and some bug fixing (#60)
Browse files Browse the repository at this point in the history 
* Add conversion classes

* Tests for wrappers

* test with linear elastic model

* Inlude update method of const-model in solver

* REfine error message

* Remove check for memory address

* Order of shear components

* Fix test

* Docstring and formatting

* Comments for the tests

* Work in suggestion from aradermacher
  • Loading branch information
@srosenbu
srosenbu authored Aug 23, 2024
1 parent e91a6f3 commit 2bc7987
Show file tree
Hide file tree
Showing 11 changed files with 753 additions and 349 deletions.
83 changes: 82 additions & 1 deletion examples/linear_elasticity/test_elasticity.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,12 @@
from linear_elasticity_model import LinearElasticityModel
from mpi4py import MPI

from fenics_constitutive import Constraint, IncrSmallStrainProblem
from fenics_constitutive import (
Constraint,
IncrSmallStrainProblem,
PlaneStrainFrom3D,
UniaxialStrainFrom3D,
)

youngs_modulus = 42.0
poissons_ratio = 0.3
Expand Down Expand Up @@ -173,6 +178,45 @@ def right_boundary(x):
analytical_stress
)

# test the converter from 3D model to uniaxial strain model
law_3d = LinearElasticityModel(
parameters={"E": youngs_modulus, "nu": poissons_ratio},
constraint=Constraint.FULL,
)
wrapped_1d_law = UniaxialStrainFrom3D(law_3d)
u_3d = df.fem.Function(V)
problem_3d = IncrSmallStrainProblem(
wrapped_1d_law,
u_3d,
[bc_left, bc_right],
1,
)
solver_3d = NewtonSolver(MPI.COMM_WORLD, problem_3d)
n, converged = solver_3d.solve(u_3d)
problem_3d.update()

# test that sigma_11 is the same as the analytical solution
assert abs(problem_3d.stress_0.x.array[0] - analytical_stress) < 1e-10 / (
analytical_stress
)
# test that the stresses of the problem with uniaxial strain model
# are the same as the stresses of the problem with the converted 3D model
assert (
np.linalg.norm(problem_3d.stress_0.x.array - problem.stress_0.x.array)
/ np.linalg.norm(problem.stress_0.x.array)
< 1e-14
)

# test that the shear stresses are zero. Since this is uniaxial strain, the
# stress can have nonzero \sigma_22 and \sigma_33 components
assert np.linalg.norm(wrapped_1d_law.stress_3d[3:6]) < 1e-14
# test that the displacement is the same in both versions
assert (
np.linalg.norm(problem_3d._u.x.array - problem._u.x.array)
/ np.linalg.norm(problem._u.x.array)
< 1e-14
)


def test_plane_strain():
# sanity check if out of plane stress is NOT zero
Expand Down Expand Up @@ -213,6 +257,43 @@ def right_boundary(x):
)
> 1e-2
)
# test the model conversion from 3D to plane strain
law_3d = LinearElasticityModel(
parameters={"E": youngs_modulus, "nu": poissons_ratio},
constraint=Constraint.FULL,
)
wrapped_2d_law = PlaneStrainFrom3D(law_3d)
u_3d = df.fem.Function(V)
problem_3d = IncrSmallStrainProblem(
wrapped_2d_law,
u_3d,
[bc_left, bc_right],
1,
)
solver_3d = NewtonSolver(MPI.COMM_WORLD, problem_3d)
n, converged = solver_3d.solve(u_3d)
problem_3d.update()
# test that the stress is nonzero in 33 direction
assert (
np.linalg.norm(
problem_3d.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[
:, 2
]
)
> 1e-2
)
# test that the displacement is the same in both versions
assert (
np.linalg.norm(problem_3d._u.x.array - problem._u.x.array)
/ np.linalg.norm(problem._u.x.array)
< 1e-14
)
# test that the stresses are the same in both versions
assert (
np.linalg.norm(problem_3d.stress_0.x.array - problem.stress_0.x.array)
/ np.linalg.norm(problem.stress_0.x.array)
< 1e-14
)


def test_plane_stress():
Expand Down
95 changes: 64 additions & 31 deletions examples/mises_plasticity/mises_plasticity_isotropic_hardening.py
View file
Original file line number Diff line number Diff line change
@@ -1,32 +1,42 @@
from __future__ import annotations

import numpy as np

from fenics_constitutive import Constraint, IncrSmallStrainModel, strain_from_grad_u

class VonMises3D(IncrSmallStrainModel):

class VonMises3D(IncrSmallStrainModel):
"""Von Mises Plasticity model with linear isotropic hardening.
Computation of trial stress state is entirely deviatoric. Volumetric part is added later
when the stress increment for the current time step is calculated """
when the stress increment for the current time step is calculated"""

def __init__(self, param: dict[str, float]):

self.xioi = np.array([[1, 1, 1, 0, 0, 0], #1dyadic1 tensor
[1, 1, 1, 0, 0, 0],
[1, 1, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0]])
self.I2 = np.zeros(self.stress_strain_dim, dtype=np.float64) # Identity of rank 2 tensor
self.xioi = np.array(
[
[1, 1, 1, 0, 0, 0], # 1dyadic1 tensor
[1, 1, 1, 0, 0, 0],
[1, 1, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
]
)
self.I2 = np.zeros(
self.stress_strain_dim, dtype=np.float64
) # Identity of rank 2 tensor
self.I2[0] = 1.0
self.I2[1] = 1.0
self.I2[2] = 1.0
self.I4 = np.eye(self.stress_strain_dim, dtype=np.float64) # Identity of rank 4 tensor
self.xpp = self.I4 - (1 / 3) * self.xioi # Projection tensor of rank 4
self.I4 = np.eye(
self.stress_strain_dim, dtype=np.float64
) # Identity of rank 4 tensor
self.xpp = self.I4 - (1 / 3) * self.xioi # Projection tensor of rank 4

self.p_ka = param["p_ka"] # kappa - bulk modulus
self.p_mu = param["p_mu"] # mu - shear modulus
self.p_y0 = param["p_y0"] # y0 - initial yield stress
self.p_y00 = param["p_y00"] # y00 - final yield stress
self.p_w = param["p_w"] # w - saturation parameter
self.p_y0 = param["p_y0"] # y0 - initial yield stress
self.p_y00 = param["p_y00"] # y00 - final yield stress
self.p_w = param["p_w"] # w - saturation parameter

def evaluate(
self,
Expand All @@ -37,28 +47,35 @@ def evaluate(
tangent: np.ndarray,
history: np.ndarray | dict[str, np.ndarray],
) -> None:

stress_view = mandel_stress.reshape(-1, self.stress_strain_dim)
tangent_view = tangent.reshape(-1, self.stress_strain_dim**2)
strain_increment = strain_from_grad_u(grad_del_u, self.constraint).reshape(-1, self.stress_strain_dim)
strain_increment = strain_from_grad_u(grad_del_u, self.constraint).reshape(
-1, self.stress_strain_dim
)
eps_n = history["eps_n"].reshape(-1, self.stress_strain_dim)
alpha = history["alpha"]

for n, eps in enumerate(strain_increment):

tr_eps = np.sum(eps[:3]) # trace of strain
eps_dev = eps - tr_eps * self.I2 / 3 # deviatoric strain
tr_eps = np.sum(eps[:3]) # trace of strain
eps_dev = eps - tr_eps * self.I2 / 3 # deviatoric strain

# deviatoric trial internal forces (trial stress)
del_sigtr = 2 * self.p_mu * eps_dev # incremental trial stress
stress_n_dev = stress_view[n] - np.sum(stress_view[n][:3]) * self.I2 / 3 # total deviatoric stress in last time step
sigtr = stress_n_dev + del_sigtr # total (deviatoric) trial stress in current time step
del_sigtr = 2 * self.p_mu * eps_dev # incremental trial stress
stress_n_dev = (
stress_view[n] - np.sum(stress_view[n][:3]) * self.I2 / 3
) # total deviatoric stress in last time step
sigtr = (
stress_n_dev + del_sigtr
) # total (deviatoric) trial stress in current time step

# norm of deviatoric trial internal forces
sigtrn = np.sqrt(np.dot(sigtr, sigtr))

# trial yield criterion
phitr = sigtrn - np.sqrt(2 / 3) * (self.p_y0 + (self.p_y00 - self.p_y0) * (1 - np.exp(-self.p_w * alpha[n])))
phitr = sigtrn - np.sqrt(2 / 3) * (
self.p_y0
+ (self.p_y00 - self.p_y0) * (1 - np.exp(-self.p_w * alpha[n]))
)

# CHECK YIELD CRITERION
# elastic - plastic step
Expand All @@ -71,12 +88,24 @@ def evaluate(
nmax = 100
# flow direction
xn = sigtr / sigtrn

# UPDATE PLASTIC MULTIPLIER VIA NEWTON-RAPHSON SCHEME
def f(x):
return sigtrn - 2 * self.p_mu * x - np.sqrt(2 / 3) * ( self.p_y0 + (self.p_y00 - self.p_y0) * (1 - np.exp(-self.p_w * (alpha[n] + np.sqrt(2 / 3) * x))))
return (
sigtrn
- 2 * self.p_mu * x
- np.sqrt(2 / 3)
* (
self.p_y0
+ (self.p_y00 - self.p_y0)
* (1 - np.exp(-self.p_w * (alpha[n] + np.sqrt(2 / 3) * x)))
)
)

def df(x):
return - 2 * self.p_mu - (2 / 3) * (self.p_y00 - self.p_y0) * self.p_w * np.exp(-self.p_w * (alpha[n] + np.sqrt(2 / 3) * x))
return -2 * self.p_mu - (2 / 3) * (
self.p_y00 - self.p_y0
) * self.p_w * np.exp(-self.p_w * (alpha[n] + np.sqrt(2 / 3) * x))

# start Newton iteration
while np.abs(xr) > tol:
Expand All @@ -89,16 +118,16 @@ def df(x):
gamma = gamma - xr / xg
# exit Newton algorithm for iteration > nmax
if it > nmax:
#print('No Convergence in Newton Raphson Iteration')
print('gamma is', xr)
# print('No Convergence in Newton Raphson Iteration')
print("gamma is", xr)
break
# end of Newton iterration

# compute tangent with converged gamma
xg = df(gamma)

# algorithmic parameters
xc1 = - 1 / xg
xc1 = -1 / xg
xc2 = gamma / sigtrn

# ELASTIC STEP
Expand All @@ -118,7 +147,11 @@ def df(x):
stress_view[n] += sh

# determine elastic-plastic moduli
aah = self.p_ka * self.xioi + 2 * self.p_mu * (1 - 2 * self.p_mu * xc2) * self.xpp + 4 * self.p_mu * self.p_mu * (xc2 - xc1) * np.outer(xn,xn)
aah = (
self.p_ka * self.xioi
+ 2 * self.p_mu * (1 - 2 * self.p_mu * xc2) * self.xpp
+ 4 * self.p_mu * self.p_mu * (xc2 - xc1) * np.outer(xn, xn)
)
tangent_view[n] = aah.flatten()

def update(self) -> None:
Expand All @@ -130,4 +163,4 @@ def constraint(self) -> Constraint:

@property
def history_dim(self) -> int:
return {"eps_n": self.constraint.stress_strain_dim(), "alpha": 1}
return {"eps_n": self.constraint.stress_strain_dim(), "alpha": 1}
Loading

0 comments on commit 2bc7987

Please sign in to comment.