Hi @bleyerj, one more question of this plastic model. Since I will apply this model for both 2D and 3D problem, I modified you model as below:
dim = 2 # or 3
def epsilon(u):
return fe.sym(fe.grad(u))
def sigma(eps_u):
return lambda*fe.tr(eps_u)*fe.Identity(dim) + 2 * mu * eps_u
def vector2tensor(X):
if dim is 2:
return fe.as_tensor([[X[0], X[3]], [X[3], X[1]]])
else:
return fe.as_tensor([[X[0], X[3], X[5]], [X[3], X[1], X[4]], [X[5], X[4], X[2]]])
def tensor2vector(X):
if dim is 2:
return fe.as_vector([X[0, 0], X[1, 1], 0, X[0, 1]])
else:
return fe.as_vector([X[0, 0], X[1, 1], X[2, 2], X[0, 1], X[1, 2], X[0, 2]])
def proj_sig(deps, old_sig, old_p):
sig_elas = vector2tensor(old_sig) + sigma(deps)
s = fe.dev(sig_elas) # deviatoric elastic stress
# isotropic elasto-plastic von Mises yield condition
sig_eq = fe.sqrt(3/2.*fe.inner(s, s))
f_elas = sig_eq - sig0 - H*old_p # plasticity criterion
f_elas_pp = (f_elas + abs(f_elas))/2 # positive part of x
dp = f_elas_pp/(3*mu+H) # plastic strain increment
n_elas = s/sig_eq*f_elas_pp/f_elas # normal vector to the final yield surface
beta = 3*mu*dp/sig_eq
new_sig = sig_elas-beta*s # final stress state
return tensor2vector(new_sig), tensor2vector(n_elas), beta, dp
def sigma_tang(eps):
N_elas = vector2tensor(n_elas)
return sigma(eps) - 3*mu*(3*mu/(3*mu+H)-beta)*fe.inner(N_elas, eps)*N_elas-2*mu*beta*fe.dev(eps)
Is it correct?
I saw you mentioned the strain tensor will be represented in a 3D fashion by appending zeros on the out-of-plane components since, even if the problem is 2D, the plastic constitutive relation will involve out-of-plane plastic strains.
Does the way I deal with 2D and 3D tensors work?
Thank you.
