# Plasticity in fenics

Dear all,

I am trying to implement plasticity in FEniCS. I find an extremely valuable tutorial here (Courtesy @bleyerj ). Thanks a ton @bleyerj . I tried to match the result of the FEniCS code with the Abaqus so that I can build upon that. But unfortunately, the displacement at the end of the simulation from FEniCS and Abaqus are not matching.

. As you can see, at the end displacement are very different. Though, the onset of plasticity is approximately at the same displacement. Am I missing something in Abaqus? For the reference, I have attached the input file of Abaqus here (for the suggestions from any of folks here who are familiar with Abaqus as well might help).

Also, are there anymore sources of plasticity in FEniCS so that I can explore further? Any pointers in that direction will be useful.

OS: Ubuntu 16.04LTS
FEniCS version: 2017.2.0

Peter

Hi,
thanks for the comment on the plasticity example. I haven’t compared to Abaqus but I get identical results between the FeniCS implementation and another completely independent implementation of FE plasticity.

I am not an experienced Abaqus user but your time goes up to 1., what is your imposed loading ? it seems from your input file that your final load is 75.751 right ? What does it correspond to? In the example, I go up to q_{lim}*\sqrt{1.1} which is around 79.4347

Hi @bleyerj,

Thanks for your reply. My bad. I did not notice that q_{lim} * \sqrt{1.1}. Now, it matches with the answer provided by you. Extremely sorry for the mess up.

Hi @bleyerj,

Firstly I want to thanks for providing this elasto-plastic tutorial. I have one simple question related to this model.
It is force-controlled model in this tutorial, could you please help me for implementing the same model with displacement-controlled? I had tried commenting F_ext and adding DirichletBC with certain value (> 0) into bc. However, the manual Newton-Raphson method cannot coverage.

Thank you.

1 Like

Hello,

using displacement controlled requires more care. First, from how I have written the demo, the Newton Raphson system is expressed on the displacement increment, so you need to add a DirichletBC with imposed value equal to the displacement increment over the time step and not the total displacement. Second, only the first iteration should have a non-zero DirichletBC all other iterations must be computed with zero DirichletBC, otherwise you will accumulate imposed displacements over all iterations and will not converge (see for instance: Custom Newton Solver problem with Dirichlet conditions)

1 Like

Hello,

Were you able to implement a solution to this problem?

Thank you.

Yes, as the link mentioned in @bleyerj 's reply, I got reasonable results by adding for bc in bcs:bc.homogenize() after system_assemble() of first increment.

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):
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.

Hi, it seems OK to me.

Hi all,

At the end of these plasticity models, how to plot the equivalent or von mises stress - sqrt(3/2.*inner(s, s)) ?

I tried something like this at the end of code, but it is not working. Can you please point out where I’m going wrong ?
Thank you !

sig_n1 = as_3D_tensor(sig)
s1 = dev(sig_n1)
sig_eq1 = sqrt(3/2.*inner(s1, s1))
# sig_avg.assign(project(sig_eq1, P0))

Vstr = VectorFunctionSpace(mesh, "DG", 1, dim=1)
von_mises = project(sig_eq1, Vstr)

# For Plotting
import matplotlib.pyplot as plt
p = plot(von_mises, mode='color')
plt.colorbar(p)
plt.show()


What do you mean when you are saying it is not working? Please supply an error message.
In general, you can also use VTKPLOTTER (New vtkplotter module for visualization) or Paraview for visualization.

@dokken

In this plasticity example, since sig was a 4 dimensional, projecting on a Vector Function Space of dimension 4 helped in visualizing each stress component.

Vstr = VectorFunctionSpace(mesh, "DG", 1, dim=4)
strDG = local_project(sig, Vstr)

import matplotlib.pyplot as plt
p = plot(strDG[0], mode='color')
plt.colorbar(p)
plt.title(r"$\sigma_{xx}$",fontsize=26)
plt.show()



But when I attmept to do a similar this for the equivalent stress, it shows the following.
How do I project the sig_eq1 ? Thank you.

RuntimeError: Don't know how to plot given object:

and projection failed:


sig_n1 = as_3D_tensor(sig)
s1 = dev(sig_n1)
sig_eq1 = sqrt(3/2.*inner(s1, s1))

# For Plotting
import matplotlib.pyplot as plt
p = plot(sig_eq1, mode='color')
plt.colorbar(p)
plt.title(r"Von Mises",fontsize=26)
plt.show()


As far as I can see \sigma^{eq}_{elas} in the plastic criterion a scalar. So projecting on an appropriate FunctionSpace and plotting should work.

Vs = FunctionSpace(mesh, "DG", 1)
sig_eq1 = project(sqrt(3/2.*inner(s1, s1)), Vs)

import matplotlib.pyplot as plt
p = plot(sig_eq1, mode='color')
plt.colorbar(p)
plt.title(r"Von Mises",fontsize=26)
plt.show()



should work unless I am missing something

Thanks for the reply. Like you said, sig_eq is a scalar, but it still shows

AssertionError: Mismatch of quadrature points!


You are dealing with functions belonging to the Quadrature-elements. Take a look at this for instance. The following code works fine for me


sig_n1 = as_3D_tensor(sig)
s1 = dev(sig_n1)
Vs = FunctionSpace(mesh, "DG", 1)
sig_eq1 = local_project(sqrt(3/2.*inner(s1, s1)), Vs)

# For Plotting
import matplotlib.pyplot as plt
p = plot(sig_eq1, mode='color')
plt.colorbar(p)
plt.title(r"Von Mises",fontsize=26)
plt.show()


1 Like

@bhaveshshrimali
Thanks a lot !
It works now. I had used project instead of local_project previously, which caused mismatch of quadrature points.

I am currently trying to solve a nonlinear return mapping equation. I try to use fsolve from scipy but keep getting the error “Result from function call is not a proper array of floats.” Is there some kind of special treatment that should be considered when solving for the equivalent plastic strain?

Hi,
for your information, we have now included a FEniCS interface for the constitutive law code generator MFront in the MFrontGenericInterfaceSupport project. Documentation and commented demos are available here: Overview of the mgis.fenics module
This enables to solve problems involving complex non-linear constitutive equations inside FEniCS including plasticity of course.

1 Like

Can someone please help me on how to plot the plastic strain at the end of the calculations ?
@bhaveshshrimali @bleyerj @dokken

I tried doing the following

import matplotlib.pyplot as plt
p1 = plot(p_avg, mode=‘color’)
plt.colorbar(p1)
plt.title(r"Plastic Strain",fontsize=26)
plt.show()

It is displays the following error.

AttributeError: ‘PolyCollection’ object has no property ‘mode’

Hi,
try this for plotting via matplotlib

import matplotlib.pyplot as plt
from matplotlib import cm
p = plot(p_avg, cmap=cm.jet)
plt.colorbar(p)
plt.show()


this works for me

1 Like