Hi! I am trying to modify the code by @bleyerj for von Mises elasto-plastic example Elasto-plastic analysis of a 2D von Mises material in order to apply Dirichlet BCs. While testing under a single element setting, the code runs fine and there are no convergence issues. Analyzing the displacement, stress, strain, and plastic strain values, everything is consistent.
Currently, I am trying to run the double notched tensile specimen test from the book Computational Methods for Plasticity (section 7.5.5). Considering symmetry, only the top half of the domain is considered. While running the code, the solver converges for the first few load steps, but then diverges. Under current setting, the solver diverges after the load is greater than 0.0187 mm. While testing the code under different settings, my guess is that in case of a non-homogeneous geometry, the solver slowly diverges during the plastic steps.
Is there something wrong with the formulation, or some changes need to be made for the solver? Thanks for the help!
import gmsh
from dolfinx import mesh, fem, io
from mpi4py import MPI
from petsc4py import PETSc
import ufl
import basix
import numpy as np
from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, set_bc
model_rank = 0
gdim = 2
if MPI.COMM_WORLD.rank == 0:
# Initialize GMSH
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
# Create a new model
gmsh.model.add("model1")
# Set OpenCASCADE kernel
gmsh.model.geo.addPoint(0, 0, 0, 1.0, 1)
gmsh.model.geo.addPoint(0, 15, 0, 1.0, 2)
gmsh.model.geo.addPoint(5, 15, 0, 1.0, 3)
gmsh.model.geo.addPoint(5, 0.001, 0, 1.0, 4)
gmsh.model.geo.addPoint(0.5, 0, 0, 1.0, 5)
# Create lines
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 5, 4)
gmsh.model.geo.addLine(5, 1, 5)
# Create curve loop and surface
gmsh.model.geo.addCurveLoop([1, 2, 3, 4, 5], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
# Synchronize to create geometry
gmsh.model.geo.synchronize()
# Physical groups
gmsh.model.addPhysicalGroup(2, [1], 6)
gmsh.model.setPhysicalName(2, 6, "surface")
gmsh.model.addPhysicalGroup(1, [1], 7)
gmsh.model.setPhysicalName(1, 7, "left")
gmsh.model.addPhysicalGroup(1, [5], 8)
gmsh.model.setPhysicalName(1, 8, "bot")
gmsh.model.addPhysicalGroup(1, [2], 9)
gmsh.model.setPhysicalName(1, 9, "top")
# Transfinite meshing constraints
gmsh.model.mesh.setTransfiniteCurve(1, 51, coef=1.0)
gmsh.model.mesh.setTransfiniteCurve(3, 51, coef=1.0)
gmsh.model.mesh.setTransfiniteCurve(2, 21, coef=1.0)
gmsh.model.mesh.setTransfiniteCurve(5, 6, coef=1.0)
gmsh.model.mesh.setTransfiniteCurve(4, 16, coef=1.0)
gmsh.model.mesh.setRecombine(2, 1)
# Generate 2D mesh
gmsh.model.mesh.generate(2)
# Write mesh file (optional)
domain, _, ft = io.gmshio.model_to_mesh(
gmsh.model, MPI.COMM_WORLD, model_rank, gdim=gdim
)
# Finalize GMSH
gmsh.finalize()
E = fem.Constant(domain, 206.9) # in MPa
nu = fem.Constant(domain, 0.29)
lmbda = E * nu / (1 + nu) / (1 - 2 * nu)
mu = E / 2.0 / (1 + nu)
K = E / (3*(1-2*nu))
sig0 = fem.Constant(domain, 0.45) # yield strength in MPa
H = fem.Constant(domain, 0.0) # hardening modulus
deg_u = 2
shape = (gdim,)
V = fem.functionspace(domain, ("P", deg_u, shape))
bottom_dofsy = fem.locate_dofs_topological(V.sub(1), 1, ft.find(8))
left_dofsx = fem.locate_dofs_topological(V.sub(0), 1, ft.find(7))
top_dofsy = fem.locate_dofs_topological(V.sub(1), 1, ft.find(9))
# Specify the displacement value
disp = fem.Constant(domain, 0.0)
# Set the boundary conditions
bcs = [
fem.dirichletbc(PETSc.ScalarType(0.0), left_dofsx, V.sub(0)),
fem.dirichletbc(PETSc.ScalarType(0.0), bottom_dofsy, V.sub(1)),
fem.dirichletbc(disp, top_dofsy, V.sub(1))
]
deg_quad = 4 # quadrature degree for internal state variable representation
W0e = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(), scheme="default", degree=deg_quad)
# Quadrature element to store 4x1 vector sized values
We = basix.ufl.quadrature_element(domain.basix_cell(), value_shape=(4,), scheme="default", degree=deg_quad)
# Function spaces are created for the two previous quadrature elements
W = fem.functionspace(domain, We)
W0 = fem.functionspace(domain, W0e)
sig = fem.Function(W)
n_elas = fem.Function(W)
beta = fem.Function(W0)
sig_eq = fem.Function(W0)
p = fem.Function(W0, name="Cumulative_plastic_strain")
dp = fem.Function(W0)
u = fem.Function(V, name="Total_displacement")
du = fem.Function(V, name="Iteration_correction")
Du = fem.Function(V, name="Current_increment")
delta_eps_p = fem.Function(W)
eps_p = fem.Function(W)
eps_p_old = fem.Function(W)
v = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
ds_top = ufl.Measure(integral_type='ds', domain=domain, subdomain_data=ft, subdomain_id=8)
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]])
# Returns the stress using lame paramters for total strain
def elastic_behavior(eps_el):
return lmbda * ufl.tr(eps_el) * ufl.Identity(3) + 2 * mu * eps_el
# Returns the stress using lame paramters for total strain minus plastic strain
def elpl_behavior(eps_el, eps_pl):
return lmbda * ufl.tr(eps_el-eps_pl) * ufl.Identity(3) + 2 * mu * (eps_el-eps_pl)
# Function to convert a vector of size 1x4 to a 2nd order tensor of size 3x3
def as_3D_tensor(X):
return ufl.as_tensor([[X[0], X[3], 0], [X[3], X[1], 0], [0, 0, X[2]]])
# Function to convert a 2nd order tensor of size 3x3 to a vector of size 1x4
def to_vect(X):
return ufl.as_vector([X[0, 0], X[1, 1], X[2, 2], X[0, 1]])
#Function returns a 2x2 matrix from a 3x3 matrix
def as_2D_tensor(X):
return ufl.as_tensor([[X[0], X[3]], [X[3], X[1] ] ])
def threeD_as_2D_tensor(X):
return ufl.as_tensor([[X[0, 0], X[0, 1]], [X[1, 0], X[1, 1]]])
# Projects a function (unknown) from one function space to another function space
# It solves a linear system of equation inorder for a function to be projected.
def project(function, space):
p = ufl.TrialFunction(space)
q = ufl.TestFunction(space)
dz = ufl.Measure('dx', domain = domain, 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()
# Interpolate quadrature function evaluates the ufl conditions in the gauss points.
def interpolate_quadrature(ufl_expr, function):
expr_expr = fem.Expression(ufl_expr, quadrature_points)
expr_eval = expr_expr.eval(domain, cells)
function.x.array[:] = expr_eval.flatten()[:]
ppos = lambda x: ufl.max_value(x, 0)
# Constitutive law for the update of strain in case of plastic behaviour
# Predictor-Corrector Algorithm:
def constitutive_update(eps, old_p, eps_p):
sig_elas = elpl_behavior(eps,as_3D_tensor(eps_p))
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))
dx = ufl.Measure(
"dx",
domain=domain,
metadata={"quadrature_degree": deg_quad, "quadrature_scheme": "default"},
)
# Weak form Residual for the Mechanical Equation
Residual = ufl.inner(eps(u_),as_3D_tensor(sig)) * dx #- ufl.inner( loading * n, u_ ) * ds
# K_t tangent form
tangent_form = ufl.inner(eps(v), sigma_tang(eps(u_))) * dx
# Exports the cell type and calculate the number of guass points and weight for the specified degree '4'
basix_celltype = getattr(basix.CellType, domain.topology.cell_type.name)
quadrature_points, weights = basix.make_quadrature(basix_celltype, deg_quad)
# Maps the entire domain and obtains the total number of gauss points that are present in the entire domain
map_c = domain.topology.index_map(domain.topology.dim)
num_cells = map_c.size_local + map_c.num_ghosts
cells = np.arange(0, num_cells, dtype=np.int32)
ngauss = num_cells * len(weights)
class CustomLinearProblem(fem.petsc.LinearProblem):
def assemble_rhs(self):
# Assemble rhs
with self._b.localForm() as b_loc:
b_loc.set(0)
fem.petsc.assemble_vector(self._b, self._L)
# Apply boundary conditions to the rhs
fem.petsc.apply_lifting(self._b, [self._a], bcs=[self.bcs])
self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
fem.petsc.set_bc(self._b, self.bcs)
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):
# Solve linear system and update ghost values in the solution
self._solver.solve(self._b, self._x)
self.u.x.scatter_forward()
#---------------------------------------------------------------------------------------------------
# problem solves for correction displacement du
tangent_problem = CustomLinearProblem(
tangent_form,
-Residual,
u=du,
bcs=bcs,
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
},
)
Nitermax, tol = 15, 1e-6 # parameters of the Newton-Raphson procedure
load = 0.01 # greater than 0 otherwise nRes / nRes0 tends to infinity
delta_u = 0.005 # load step
total_u = 0.17 # total displacement
iteration_blocks = []
while load <=1:
disp.value = load*total_u
print(f"load = {load*total_u}")
tangent_problem.assemble_rhs()
nRes0 = tangent_problem._b.norm()
nRes = nRes0
Du.x.array[:] = 0
niter = 0
current_niters = []
current_nRes = []
print("nRes0", nRes0)
while nRes / nRes0 > tol and niter < Nitermax:
disp.value = 0.0
tangent_problem.assemble_lhs()
tangent_problem.solve_system()
Du.vector.axpy(1, du.vector)
Du.x.scatter_forward()
eps_t = eps(Du)
sig_, n_elas_, beta_, dp_ , delta_eps_p_ = constitutive_update(eps_t, p, eps_p_old)
# interpolate the new stresses and internal state variables
interpolate_quadrature(sig_, sig)
interpolate_quadrature(n_elas_, n_elas)
interpolate_quadrature(beta_, beta)
# compute the new residual
tangent_problem.assemble_rhs()
nRes = tangent_problem._b.norm()
print(f"nRes: {nRes}")
current_niters.append(niter)
current_nRes.append(nRes)
niter += 1
# Update the previous plastic strain
interpolate_quadrature(dp_, dp)
p.vector.axpy(1, dp.vector)
interpolate_quadrature(delta_eps_p_, delta_eps_p)
eps_p.vector.axpy(1, delta_eps_p.vector)
eps_p_old.x.array[:] = eps_p.x.array[:]
iteration_blocks.append((current_niters, current_nRes))
load += delta_u