Implementing dynamic analysis for von-Mises plasticity

General
1

Hello,

I’ve implemented a dynamic analysis for the von-Mises plasticity based on vonMises_Plasticity and Elastodynamics tutorials. The simulation works well in the elastic region but when it enters the plastic region, the residual becomes very large and does not converge in the first increment in the plastic region.

Here are the main functions in this approach. The other necessary functions (integration scheme, etc.) are not changed. I am not sure if I have to consider anything additional for the plasticity.

def eps(v):
    e = sym(grad(v))
    return as_tensor([[e[0, 0], e[0, 1], 0],
                      [e[0, 1], e[1, 1], 0],
                      [0, 0, 0]])

def sigma(eps_el):
    return lmbda*tr(eps_el)*Identity(3) + 2*mu*eps_el

def as_3D_tensor(X):
    return as_tensor([[X[0], X[3], 0],
                      [X[3], X[1], 0],
                      [0, 0, X[2]]])

ppos = lambda x: (x+abs(x))/2.
def proj_sig(deps, old_sig, old_p):
    sig_n = as_3D_tensor(old_sig)
    sig_elas = sig_n + sigma(deps)
    s = dev(sig_elas)
    sig_eq = sqrt(3/2.*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
    return as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \
           as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \
           beta, dp

def sigma_tang(e):
    N_elas = as_3D_tensor(n_elas)
    return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)

def m(u, u_):
    return rho*inner(u, u_)*dxm

def k(u, u_):
    return inner(eps(u_), as_3D_tensor(sig))*dxm + inner(eps(u), sigma_tang(eps(u_)))*dxm

def c(u, u_):
    return eta_m*m(u, u_) + eta_k*k(u, u_)

def Wext(u_):
    return loading*dot(n, u_)*ds(4)

res = m(avg(a_old, a_new, alpha_m), u_) + c(avg(v_old, v_new, alpha_f), u_) \
       + k(avg(u_old, v, alpha_f), u_) - Wext(u_)

Below is the load-displacement graph from the vonMises tutorial. The geometry enters plastic region after 0.003
image
Below is the residual when the displacement is around 0.006 (plastic).

image
image431×530 9.92 KB

I have tried implementing the same problem using MFront library instead of directly implementing the material law. I still get the same behaviour. That’s why I tried with direct implementation of material law to see where the problem is.

Any suggestion would be highly appreciated. Thanks in advance!

1 Like
2

Why do you have both sig and sig_tang in the stiffness term ? It should only be sig.

3

I have modified the stiffness term and residual as below

def k(u_):
    return inner(eps(u_), as_3D_tensor(sig))*dxm

res = m(avg(a_old, a_new, alpha_m), u_) + c(avg(v_old, v_new, alpha_f), u_) \
       + k(u_) - Wext(u_)
a_form = inner(eps(v), sigma_tang(eps(u_)))*dxm
L_form = -res
...
for (i, t) in enumerate(np.diff(time)):
    t = time[i+1]
    loading.t = t-float(alpha_f*dt)
    A, Res = assemble_system(a_form, L_form, bc)
    ...

But I get the following error

image
image1094×391 21.7 KB

This is due to using TrialFunction in nonlinear problems right? I changed the TrialFunction to Function and the system cannot be assembled then

image
image1079×450 15 KB

Here is the list of function spaces

deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)

sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)

u_old = Function(V)
v_old = Function(V)
a_old = Function(V)

Should I also update the displacement (avg(u_old, du, alpha_f)) in the sig_tang?

4

Hi,

Bleyer wrote a package that can be useful for plasticity. See his 2019 and 2020 papers. I tried the legacy version of the package, but had to switch to its dolfinx version due to incompatibility.

If you don’t mind posting a minimal working code I can try running the bugs myself.