This demo is the third installment in a series of two previous demos considering the implementation of von Mises elastoplasticity in FEniCS published on Numerical tours of continuum mechanics using FEniCS:

Elasto-plastic analysis of a 2D von Mises material: this first installment presents a pure FEniCS implementation of von Mises plasticity which uses UFL operators for implementing the return mapping procedure. It is therefore not general since it relies heavily on the fact that the return mapping can be expressed in closed-form analytic expression in the case of von Mises plasticity with linear isotropic hardening. It nevertheless shows how to integrate the constitutive relation using

`Quadrature`

elements, including the use of consistent tangent operators, inside a global Newton-Raphson method.Elasto-plastic analysis implemented using the MFront code generator: this second installment is in fact the original demo which led to the development of the proposed FEniCS/MFront interface. It shows how to rely on MFront for the constitutive relation update using the MGIS interface. The way how information is exchanged between FEniCS and MFront in this demo is extremely similar to how the

`MFrontNonlinearProblem`

class is implemented.

The present demo therefore directly follows up on this second installment by offering to the use a much more compact syntax, especially avoiding the cumbersome definition of `Quadrature`

spaces and the implementation of a Newton method by hand. The mechanical problem is exactly the same as in these two demos, namely the expansion of a hollow cylinder under internal pressure. The main difference in the present demo is that the same problem will also be solved in axisymmetric conditions, referring therefore to the expansion of a hollow sphere under internal pressure.

Source files:

- Jupyter notebook: mgis_fenics_small_strain_elastoplasticity.ipynb
- Python file: mgis_fenics_small_strain_elastoplasticity.py
- MFront behaviour file: IsotropicLinearHardeningPlasticity.mfront
- MFront behaviour file: IsotropicLinearHardeningPlasticityBrick.mfront

In this document, we consider a \(J_2\) isotropic plastic model with linear isotropic hardening. As the hardening is linear and the von Mises stress is quadratic, the behaviour integration is explicit, as described here.

The implementation used here is modified to use material properties instead of parameters. The difference between material properties and parameters is that material properties are allowed to be distinct on each integration points and may evolve over time, whereas parameters are global values.

`StandardElastoViscoPlasticity`

brickAnother approach, less efficient in this case, is to use the `StandardElastoViscoPlasticity`

brick which allows for a very compact and declarative implementation but relies on a general implicit scheme.

The advantage of using the `StandardElastoViscoPlasticity`

brick is that the user may rapidely adapt this document to try more advanced plastic behaviours.

```
@DSL Implicit;
@Behaviour IsotropicLinearHardeningPlasticityBrick;
@Author Thomas Helfer/Jérémy Bleyer;
@Date 07 / 04 / 2020;
@Algorithm NewtonRaphson;
@Epsilon 1.e-14;
@Theta 1.;
@MaterialProperty stress s0;
"YieldStress");
s0.setGlossaryName(@MaterialProperty stress H0;
"HardeningSlope");
H0.setEntryName(
@Brick StandardElastoViscoPlasticity{
"Hooke",
stress_potential : "Plastic" {
inelastic_flow : "Mises",
criterion : "Linear" {H : "H0", R0 : "s0"}
isotropic_hardening :
} };
```

Since we did not explicitely specify the values of the elastic properties, those are automatically declared as material properties. To use material properties to define yield stress and the hardening slope, one uses here a small trick, which consists in defining the material properties before the brick declaration and declare the plastic coefficients `H`

and `R0`

as functions of those material properties.

`FEniCS`

implementationWe first load the mesh and define values for the material properties. We will see later how to pass them to MFront. We will use a standard \(P_2\) Lagrange interpolation for the displacement field and roller conditions are defined on the horizontal (tagged `1`

) and vertical boundaries (tagged `3`

) of the domain. The inner surface is tagged as `4`

in the `facets`

MeshFunction.

```
%matplotlib notebook
import matplotlib.pyplot as plt
from dolfin import *
import mgis.fenics as mf
import numpy as np
import ufl
= 1.3, 1. # external/internal radius
Re, Ri = 70e3, 0.3 # elastic parameters
E, nu = 250. # yield strength
sig0 = E/100. # hardening slope
Et = E*Et/(E-Et)
H
= Mesh("meshes/thick_cylinder.xml")
mesh
plt.figure()=0.5)
plot(mesh, linewidth
plt.show()
= MeshFunction("size_t", mesh, "meshes/thick_cylinder_facet_region.xml")
facets = Measure('ds', subdomain_data=facets)
ds
= VectorFunctionSpace(mesh, "CG", 2)
V = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]
bc = FacetNormal(mesh) n
```

`<IPython.core.display.Javascript object>`

We now define the loading expression which consists of a uniform internal pressure \(q\) which will be progressively increased from 0 to a value slightly larger than \(q_{lim}\), the analytical limit load in the case of perfect plasticity, which is known to be \(q_{lim} = \dfrac{2}{\sqrt{3}}\ln\left(\dfrac{R_e}{R_i}\right)\sigma_0\) for the hollow cylinder (plane strain conditions) and \(q_{lim} = 2\ln\left(\dfrac{R_e}{R_i}\right)\sigma_0\) for the hollow sphere (axisymmetric conditions) and in which \(R_i\) (resp. \(R_e\)) denotes the internal (resp. external) radius.

The `MFrontNonlinearMaterial`

instance is loaded from the `MFront`

`IsotropicLinearHardeningPlasticity`

behaviour. Material properties can be defined by passing a dictionary to the `material_properties`

keyword argument. Note that the keys must match the entry names defined in the `MFront`

behaviour. Temperature is not involved in this behaviour so we will let it be automatically registered as a constant default value.

The `MFrontNonlinearProblem`

instance is then defined by specifying the quadrature degree used for the computation of the constitutive relation. The present computation is a standard quasi-static mechanical problem involving only the stress and strains as the pair of dual flux/gradient variables. As in the heat equation example, we should therefore register the gradient variable named `Strain`

as `sym(grad(u))`

in plane strain conditions. However, `Strain`

is again a predefined gradient variable so that we can rely on automatic registration. This is especially efficient in the present case where the registration is different depending on the hypothesis. For the `"plane_strain"`

hypothesis, we get the following message:

`Automatic registration of 'Strain' as sym(grad(Displacement)).`

whereas for the `"axisymmetric"`

hypothesis (see the `FEniCS`

tour on axisymmetric computations), we get:

```
Automatic registration of 'Strain' as sym([
[(grad(Displacement[0]))[0], (grad(Displacement[0]))[1], 0],
[(grad(Displacement[1]))[0], (grad(Displacement[1]))[1], 0],
[0, 0, Displacement[0] / x[0]]
]).
```

in which the first two components refer to the \(r\) and \(z\) directions whereas the third component refers to the out-of-plane \(\theta\) direction.

The linear form corresponding to the work of external forces can be specified through the `set_loading`

method, namely in the present case:

\[ W_{ext}(u) = \int_{r=R_i} -q(t)\boldsymbol{n}\cdot\boldsymbol{u} \text{ dS} \]

Note that we pass directly the action of the form on the solution `u`

rather than the linear form itself, in order to avoid the need of defining a `TestFunction`

. Note also the use of the additional integration measure \(2\pi r\) in the axisymmetric case.

Finally, load-stepping is done very easily by just updating the necessary expressions involved in the linear form definition (in the present case) or in the boundary conditions. At each step, the `solve`

method is called on the displacement field, no updates are necessary from the user side.

State variables can also be retrieved using the `get_state_variable`

method, again with names matching the entry names defined in the `MFront`

behaviour. In the present case, state variables include the elastic strain (represented by a 4-dimensional vector) and the equivalent plastic strain \(p\) (a scalar). The `get_state_variable`

method returns a `dolfin.Function`

defined on the `Quadrature`

space corresponding to the problem `quadrature_degree`

. Such functions must therefore be projected on function spaces suitable for visualization or other post-processing purposes, such as (DG, 0) in the present case for the equivalent plastic strain using the `project_on`

optional keyword of the `get_state_variable`

method.

```
for hypothesis in ["plane_strain", "axisymmetric"]:
= Function(V, name="Displacement")
u if hypothesis == "plane_strain":
= "Expansion of a thick cylinder in plane strain conditions"
title = float(2/sqrt(3)*ln(Re/Ri)*sig0)
q_lim = 1
measure elif hypothesis == "axisymmetric":
= "Expansion of a thick sphere in axisymmetric conditions"
title = SpatialCoordinate(mesh)
x = float(2*ln(Re/Ri)*sig0)
q_lim = 2*pi*abs(x[0])
measure print(title + "\n"+"-"*56)
= Expression("-q*t", q=q_lim, t=0, degree=2)
loading
= {"YoungModulus": E,
mat_prop "PoissonRatio": nu,
"HardeningSlope": H,
"YieldStrength": sig0}
= mf.MFrontNonlinearMaterial("./src/libBehaviour.so",
material "IsotropicLinearHardeningPlasticity",
=hypothesis,
hypothesis=mat_prop)
material_properties= mf.MFrontNonlinearProblem(u, material, quadrature_degree=2, bcs=bc)
problem *dot(n, u)*measure*ds(4))
problem.set_loading(loading
= problem.get_state_variable("EquivalentPlasticStrain")
p assert (ufl.shape(p) == ())
assert (p.ufl_element().family() == "Quadrature")
= problem.get_state_variable("ElasticStrain")
epsel assert (ufl.shape(epsel)==(4, ))
= XDMFFile("results/plasticity_{}_results.xdmf".format(hypothesis))
file_results "flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.parameters[
= 20
Nincr = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
load_steps = np.zeros((Nincr+1, 2))
results for (i, t) in enumerate(load_steps):
= t
loading.t
problem.solve(u.vector())= problem.get_state_variable("EquivalentPlasticStrain", project_on=("DG", 0))
p0
file_results.write(u, t)
file_results.write(p0, t)
+1, :] = (u(Ri, 0)[0], t)
results[i
plt.figure()
plt.title(title)0], results[:, 1], "-o")
plt.plot(results[:, "Displacement of inner boundary")
plt.xlabel(r"Applied pressure $q/q_{lim}$")
plt.ylabel( plt.show()
```

```
Expansion of a thick cylinder in plane strain conditions
--------------------------------------------------------
Behaviour 'IsotropicLinearHardeningPlasticity' has not been found in './src/libBehaviour.so'.
Attempting to compile 'IsotropicLinearHardeningPlasticity.mfront' in './'...
Automatic registration of 'Strain' as sym(grad(Displacement)).
Automatic registration of 'Temperature' as a Constant value = 293.15.
<IPython.core.display.Javascript object>
```

```
Expansion of a thick sphere in axisymmetric conditions
--------------------------------------------------------
Automatic registration of 'Strain' as sym([
[(grad(Displacement[0]))[0], (grad(Displacement[0]))[1], 0],
[(grad(Displacement[1]))[0], (grad(Displacement[1]))[1], 0],
[0, 0, Displacement[0] / x[0]]
]).
Automatic registration of 'Temperature' as a Constant value = 293.15.
<IPython.core.display.Javascript object>
```

We can check that the obtained results match those of the first two installments, both in terms of load-displacement curves, number of iterations and residual values for the plane strain case. In the axisymmetric case, we obtain similar results with plastic flow occuring at a critical pressure very close to the analytical limit load \(q_{lim}\) in the perfectly plastic case.