Hello everyone, i am trying to implement the elastic plastic analysis code Elasto-plastic analysis of a 2D von Mises material — Computational Mechanics Numerical Tours with FEniCSx
In here I planned to visualise the cumulative plastic strain, however upon visualising i am looking at the results of cumulative plastic strain, cumulative total strain were looking the same. Even the maximum value of magnitude of strain tensor is also in same range as this cumulative plastic strain. So i am surprised, whether i am making any blender in the post processing as the value of plastic strain is looking same as total strain.
My code is
import ufl
import gmsh
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import basix
from dolfinx import fem, io, mesh
import dolfinx.fem.petsc
from dolfinx.io import gmshio
# ------------------------------------------------------------
# Read mesh from gmsh (UNCHANGED PART)
# ------------------------------------------------------------
msh, cell_markers, facet_markers = gmshio.read_from_msh(
"2D_spot_welding_geometry.msh", MPI.COMM_WORLD, gdim=2
)
gdim = 2
# ------------------------------------------------------------
# Function space (same order as nonlinear code)
# ------------------------------------------------------------
deg_u = 2
V = fem.functionspace(msh, ("P", deg_u, (gdim,)))
Vx, \_ = V.sub(0).collapse()
Vy, \_ = V.sub(1).collapse()
# ------------------------------------------------------------
# Boundary conditions (FROM YOUR SECOND CODE)
# ------------------------------------------------------------
bottom_surface_dofs = fem.locate_dofs_topological(
(V.sub(0), Vx), msh.topology.dim - 1, facet_markers.find(14)
)
bottom_surface_dofsy = fem.locate_dofs_topological(
(V.sub(1), Vy), msh.topology.dim - 1, facet_markers.find(14)
)
u0x = fem.Function(Vx)
u0y = fem.Function(Vy)
bc_fixed_x = fem.dirichletbc(u0x, bottom_surface_dofs, V.sub(0))
bc_fixed_y = fem.dirichletbc(u0y, bottom_surface_dofsy, V.sub(1))
bcs = \[bc_fixed_x, bc_fixed_y\]
# ------------------------------------------------------------
# Traction loading (facet 12 from your mesh)
# ------------------------------------------------------------
n = ufl.FacetNormal(msh)
loading = fem.Constant(msh, 0.0)
ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_markers)
# ------------------------------------------------------------
# Material parameters (FROM FIRST CODE STYLE)
# ------------------------------------------------------------
E_steel = 209040.0
E_copper = 87800.0
# Create a function for Young's modulus and Poisson's ratio
E = fem.Function(fem.functionspace(msh, ("DG", 0)))
E.x.array\[cell_markers.find(17)\] = E_steel
E.x.array\[cell_markers.find(18)\] = E_copper
nu_copper = 0.326
nu_steel = 0.2897
nu = fem.Function(fem.functionspace(msh, ("DG", 0)))
nu.x.array\[cell_markers.find(17)\] = nu_steel
nu.x.array\[cell_markers.find(18)\] = nu_copper
lmbda = E \* nu / (1 + nu) / (1 - 2 \* nu)
mu = E / 2.0 / (1 + nu)
sig0 = fem.Constant(msh, 275.0) # Yield stress
Et = E / 100.0
H = E \* Et / (E - Et)
# ------------------------------------------------------------
# Strain and stress operators
# ------------------------------------------------------------
def eps(v):
e = ufl.sym(ufl.grad(v))
return ufl.as_tensor(\[\[e\[0, 0\], e\[0, 1\], 0\],
\[e\[0, 1\], e\[1, 1\], 0\],
\[0, 0, 0\]\])
def elastic_behavior(eps_el):
return lmbda \* ufl.tr(eps_el) \* ufl.Identity(3) + 2 \* mu \* eps_el
def as_3D_tensor(X):
return ufl.as_tensor(\[\[X\[0\], X\[3\], 0\],
\[X\[3\], X\[1\], 0\],
\[0, 0, X\[2\]\]\])
def to_vect(X):
return ufl.as_vector(\[X\[0, 0\], X\[1, 1\], X\[2, 2\], X\[0, 1\]\])
def equiv_strain(eps):
eps_dev = ufl.dev(eps)
return ufl.sqrt(2.0/3.0 \* ufl.inner(eps_dev, eps_dev))
ppos = lambda x: ufl.max_value(x, 0)
# ------------------------------------------------------------
# Quadrature spaces for internal variables
# ------------------------------------------------------------
deg_quad = 2
W0e = basix.ufl.quadrature_element(
msh.basix_cell(), value_shape=(), scheme="default", degree=deg_quad
)
We = basix.ufl.quadrature_element(
msh.basix_cell(), value_shape=(4,), scheme="default", degree=deg_quad
)
W = fem.functionspace(msh, We)
W0 = fem.functionspace(msh, W0e)
sig = fem.Function(W)
sig_old = fem.Function(W)
n_elas = fem.Function(W)
beta = fem.Function(W0)
p = fem.Function(W0)
dp = fem.Function(W0)
u = fem.Function(V, name="Displacement")
du = fem.Function(V)
Du = fem.Function(V)
delta_eps_p = fem.Function(W)
eps_p = fem.Function(W)
eps_p_old = fem.Function(W)
P0 = fem.functionspace(msh, ("DG", 0))
p_avg = fem.Function(P0, name="Plastic_strain")
eps_eq = fem.Function(W0) # cumulative equivalent total strain
deps_eq = fem.Function(W0) # increment
eps_eq.vector.set(0.0)
# ------------------------------------------------------------
# Constitutive update (UNCHANGED FROM FIRST CODE)
# ------------------------------------------------------------
def constitutive_update(Δε, old_sig, old_p, eps_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + elastic_behavior(Δε)
s = ufl.dev(sig_elas)
sig_eq = ufl.sqrt(3 / 2.0 \* ufl.inner(s, s))
f_elas = sig_eq - sig0 - H \* old_p
dp = ppos(f_elas) / (3 \* mu + H)
n_elas = s / sig_eq \* ppos(f_elas) / f_elas
beta = 3 \* mu \* dp / sig_eq
new_sig = sig_elas - beta \* s
delta_eps_p = dp \* (3 / 2) \* (s / sig_eq)
return to_vect(new_sig), to_vect(n_elas), beta, dp, to_vect(delta_eps_p)
def sigma_tang(eps):
N_elas = as_3D_tensor(n_elas)
return (
elastic_behavior(eps)
- 3 \* mu \* (3 \* mu / (3 \* mu + H) - beta) \* ufl.inner(N_elas, eps) \* N_elas
- 2 \* mu \* beta \* ufl.dev(eps)
)
# ------------------------------------------------------------
# Weak forms (NONLINEAR FORM WITH TRACTION ON FACET 12)
# ------------------------------------------------------------
dx = ufl.Measure(
"dx", domain=msh,
metadata={"quadrature_degree": deg_quad, "quadrature_scheme": "default"}
)
Residual = ufl.inner(eps(ufl.TestFunction(V)), as_3D_tensor(sig)) \* dx \\
- ufl.inner(-loading \* n, ufl.TestFunction(V)) \* ds(12)
tangent_form = ufl.inner(eps(ufl.TrialFunction(V)), sigma_tang(eps(ufl.TestFunction(V)))) \* dx
basix_celltype = getattr(basix.CellType, msh.topology.cell_name())
quadrature_points, weights = basix.make_quadrature(basix_celltype, deg_quad)
map_c = msh.topology.index_map(msh.topology.dim)
num_cells = map_c.size_local + map_c.num_ghosts
cells = np.arange(0, num_cells, dtype=np.int32)
def interpolate_quadrature(ufl_expr, function):
expr_expr = fem.Expression(ufl_expr, quadrature_points)
expr_eval = expr_expr.eval(msh, cells)
function.x.array\[:\] = expr_eval.flatten()\[:\]
# ------------------------------------------------------------
# Newton linear problem (FROM FIRST CODE)
# ------------------------------------------------------------
class CustomLinearProblem(fem.petsc.LinearProblem):
def assemble_rhs(self, u=None):
with self.\_b.localForm() as b_loc:
b_loc.set(0)
fem.petsc.assemble_vector(self.\_b, self.\_L)
x0 = \[\] if u is None else \[u.vector\]
fem.petsc.apply_lifting(self.\_b, \[self.\_a\], bcs=\[self.bcs\], x0=x0, scale=1.0)
self.\_b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
x0 = None if u is None else u.vector
fem.petsc.set_bc(self.\_b, self.bcs, x0, scale=1.0)
def assemble_lhs(self):
self.\_A.zeroEntries()
fem.petsc.assemble_matrix_mat(self.\_A, self.\_a, bcs=self.bcs)
self.\_A.assemble()
def solve_system(self):
self.\_solver.solve(self.\_b, self.\_x)
self.u.x.scatter_forward()
tangent_problem = CustomLinearProblem(
tangent_form,
-Residual,
u=du,
bcs=bcs,
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
},
)
# Load stepping
Nitermax, tol = 100, 1e-6
Nincr = 20
load_steps = np.linspace(0, 1.0, Nincr + 1)\[1:\]
sig.vector.set(0.0)
sig_old.vector.set(0.0)
eps_p.vector.set(0.0)
eps_p_old.vector.set(0.0)
p.vector.set(0.0)
u.vector.set(0.0)
n_elas.vector.set(0.0)
beta.vector.set(0.0)
Δε = eps(Du)
sig\_, n_elas\_, beta\_, dp\_, delta_eps_p\_ = constitutive_update(Δε, sig_old, p, eps_p_old)
deps_eq\_ = equiv_strain(Δε)
def as_2D_tensor(X):
return ufl.as_tensor(\[\[X\[0\], X\[3\]\], \[X\[3\], X\[1\] \] \])
def project(function, space):
p = ufl.TrialFunction(space)
q = ufl.TestFunction(space)
dz = ufl.Measure('dx', domain = msh, metadata={"quadrature_degree": 2, "quadrature_scheme": "default"})
a = ufl.inner(p, q) \* dz
L = ufl.inner(function, q) \* dz
problem = fem.petsc.LinearProblem(a, L, bcs = \[\])
return problem.solve()
V_strain = fem.functionspace(msh, ("DG", 0, (3,3)))
strain_output = io.XDMFFile(MPI.COMM_WORLD, "strain.xdmf", "w")
strain_output.write_mesh(msh)
cu_strainp_output = io.XDMFFile(MPI.COMM_WORLD, "Cumulative_Plastic_strain.xdmf", "w")
cu_strainp_output.write_mesh(msh)
totalstrain_output = io.XDMFFile(MPI.COMM_WORLD, "Cumulative_Total_Strain.xdmf", "w")
totalstrain_output.write_mesh(msh)
for t in load_steps:
loading.value = t \* 500.0 # traction magnitude (same scale as your old T)
tangent_problem.assemble_rhs()
nRes0 = tangent_problem.\_b.norm()
nRes = nRes0
Du.x.array\[:\] = 0
niter = 0
while nRes / nRes0 > tol and niter < Nitermax:
tangent_problem.assemble_lhs()
tangent_problem.solve_system()
Du.vector.axpy(1, du.vector)
Du.x.scatter_forward()
# interpolate the new stresses and internal state variables
interpolate_quadrature(sig\_, sig)
interpolate_quadrature(n_elas\_, n_elas)
interpolate_quadrature(beta\_, beta)
tangent_problem.assemble_rhs()
nRes = tangent_problem.\_b.norm()
niter += 1
# Update global state
u.vector.axpy(1, Du.vector)
u.x.scatter_forward()
interpolate_quadrature(dp\_, dp)
p.vector.axpy(1, dp.vector)
interpolate_quadrature(deps_eq\_, deps_eq)
eps_eq.vector.axpy(1, deps_eq.vector)
# Update the previous stress
sig_old.x.array\[:\] = sig.x.array\[:\]
strain_expr = fem.Expression(eps(u), V_strain.element.interpolation_points())
strain = fem.Function(V_strain)
strain.interpolate(strain_expr)
#p_avg.interpolate(project(p, P0))
p_avg = project(p, P0)
eps_eq_avg = project(eps_eq, P0)
# Save results
strain_output.write_function(strain,t)
cu_strainp_output.write_function(p_avg,t)
totalstrain_output.write_function(eps_eq_avg, t)
print("max cumulative total strain :", np.max(eps_eq_avg.x.array))
print("max cumulative plastic strain :", np.max(p_avg.x.array))
Any suggestion in this is appreciated, whether the visualisation that i am carrying out is wrong regarding plastic strain or any other issues that causing this. Thanks for all your time and considerations.