How to set the constitutive equation for the elasto-plastic study of 3d structures

General
1

Dear community,
I am trying to simulate the elasto-plastic mechanics of 3d structures by some basic constitutive equation. For example bilinear model(as fig. 1 and fig. 2) and Hollomon model(fig. 3).

fig1
fig1292×150 11.2 KB

fig. 1
fig2
fig2279×213 5.18 KB

fig. 2
fig3
fig. 3
Refer to the content in this page https://jsdokken.com/dolfinx-tutorial/chapter2/nonlinpoisson_code.html#newtons-method. I wrote the following code for bilinear model:

elastic parameters
E=110
ν=0.3
lambda_ = E*ν/(1+ν)/(1-2*ν)
mu=E/(2+2*ν)
#plastic parameters
E_p=50
ν_p=0.3
lambda_p = E_p*ν_p/(1+ν_p)/(1-2*ν_p)
mu_p=E_p/(2+2*ν_p)
Sy=dfx.fem.Constant(domain, scalar(0.865))#yield stress

def epsilon(u):
    return ufl.sym(ufl.grad(u))  

def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)
def sigma_p(u):
    return lambda_p * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu_p * epsilon(u)

def rJ2(A):
    """Square root of J2 invariant of tensor A"""
    J2 = 1 / 2 * ufl.inner(A, A)
    rJ2 = ufl.sqrt(J2)
    return ufl.conditional(rJ2 < 1.0e-12, 0.0, rJ2)

def von_Mises(u):
    return ufl.sqrt(3) * rJ2(ufl.dev(sigma(u)))

def plasticity_model(strain,u):
    v_M_criterion=ufl.gt(von_Mises(u), Sy)
    stress_p= E_p * strain 
    stress_e=E * strain
    return ufl.conditional(v_M_criterion,stress_p,stress_e)

# weak form
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
strain = ufl.sym(ufl.grad(u))
stress = plasticity_model(strain,u)
F = ufl.inner(stress, ufl.grad(v)) * ufl.dx 
uh = fem.Function(V)
dis=[-0.01*i for i in range(1,5)]
xx_force=[]
von_Mises_force=[]
for i in dis:
    u_D2 = np.array([i, 0, 0], dtype=default_scalar_type)
    nodeids2 = fem.locate_dofs_topological(V, fdim, boundary_facets)
    bc2 = fem.dirichletbc(u_D2, nodeids2, V)
    problem = NonlinearProblem(F, uh, bcs=[bc, bc2])
    #bc is the bottom boundary of the domain, which is fixed. bc2 is the upper surface and applied displacement.

    solver = NewtonSolver(MPI.COMM_WORLD, problem)
    solver.convergence_criterion = "incremental"
    solver.rtol = 1e-6
    solver.report = True
    ksp = solver.krylov_solver
    opts = PETSc.Options()
    option_prefix = ksp.getOptionsPrefix()
    opts[f"{option_prefix}ksp_type"] = "gmres"
    opts[f"{option_prefix}ksp_rtol"] = 1.0e-8
    opts[f"{option_prefix}pc_type"] = "hypre"
    opts[f"{option_prefix}pc_hypre_type"] = "boomeramg"
    opts[f"{option_prefix}pc_hypre_boomeramg_max_iter"] = 1
    opts[f"{option_prefix}pc_hypre_boomeramg_cycle_type"] = "v"
    ksp.setFromOptions()

the error messages is shown in fig. 4

fig 4
fig 4774×366 89.9 KB

I’m a rookie for Fenicsx, so could you please solve my problems and tell me how to implement elastoplastic simulation I need.

2

See here Elasto-plastic analysis of a 2D von Mises material — Computational Mechanics Numerical Tours with FEniCSx