This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License
Elasto-plastic analysis implemented using the MFront code generator¶
This numerical tour has been written in collaboration with Thomas Helfer (thomas.helfer@cea.fr), MFront’s main developper.
Introduction¶
This example is concerned with the incremental analysis of an elasto-plastic von Mises material. The behaviour is implemented using the code generator tool MFront [HEL2015] (see http://tfel.sourceforge.net). The behaviour integration is handled by the MFrontGenericInterfaceSupport project (see https://github.com/thelfer/MFrontGenericInterfaceSupport) which allows to load behaviours generated by MFront and handle the behaviour integration. Other information concerning FEniCS binding in MFront can be found at https://thelfer.github.io/mgis/web/FEniCSBindings.html, in particular binding through the C++ interface and inspired by the Fenics Solid Mechanics project is discussed.
The considered example is exactly the same as the pure FEniCS tutorial Elasto-plastic analysis of a 2D von Mises material so
that both implementations can be compared. Many implementation steps are
common to both demos and will not be discussed again here, the reader can refer
to the previous tutorial for more details. The sources can be downloaded from
plasticity_mfront.zip
.
Let us point out that a pure FEniCS implementation was possible only for this specific constitutive law, namely von Mises plasticity with isotropic linear hardening, since the return mapping step is analytical in this case and can be expressed using simple UFL operators. The use of MFront can therefore make it possible to consider more complex material laws. Indeed, one only has to change the names of the behaviour and library generated by MFront to change the material behaviours, which can be arbitraly complex (although limited to small strains). This will be illustrated in forthcoming tutorials. Interested users may have a look at the examples of the MFront gallery to have a small overview of MFront abilities: http://tfel.sourceforge.net/gallery.html
Prerequisites¶
In order to run this numerical tour, you must first install the TFEL library on which MFront is built as well as MGIS (MFrontGenericInterfaceSupport). Please note that development versions are used for now.
To proceed with the installation, we provide the following installation shell script install.sh
.
Before running it, please install the necessary prerequisites mentioned in the script, e.g. on Ubuntu run:
> sudo apt-get install cmake libboost-all-dev g++ gfortran
> sudo apt-get install git libqt5svg5-dev qtwebengine5-dev
> sudo apt-get install python3-matplotlib
Once the installation script finished installing TFEL and MGIS, you can use MGIS by running the following command:
> source <PREFIX>/codes/mgis/master/install/env.sh
in which <PREFIX>
is the installation directory you have chosen.
As recalled later in the tutorial, the MFront behaviour file must first be compiled before running this Python script. The compilation command is:
> mfront --obuild --interface=generic IsotropicLinearHardeningPlasticity.mfront
Problem position¶
The present implementation will heavily rely on Quadrature
elements to represent
previous and current stress states, internal state variables (cumulated
equivalent plastic strain \(p\) in the present case) as well as components
of the tangent stiffness matrix. The use of quadrature
representation is therefore
still needed to circumvent the known issue with uflacs
:
from __future__ import print_function
from dolfin import *
import numpy as np
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
In order to bind MFront behaviours with FEniCS Python interface we will first load
the MFrontGenericInterfaceSupport Python package named mgis
and more precisely
the behaviour
subpackage:
import mgis.behaviour as mgis_bv
As before, we load a Gmsh
generated mesh representing a portion of a thick cylinder:
Re, Ri = 1.3, 1. # external/internal radius
mesh = Mesh("thick_cylinder.xml")
facets = MeshFunction("size_t", mesh, "thick_cylinder_facet_region.xml")
ds = Measure('ds')[facets]
We now define appropriate function spaces. Standard CG-space of degree 2 will still be used for the displacement whereas various Quadrature spaces are considered for:
- stress/strain-like variables (represented here as 4-dimensional vector since the \(zz\) component must be considered in the plane strain plastic behaviour)
- scalar variables for the cumulated plastic strain
- the consistent tangent matrix represented here as a tensor of shape 4x4
As in the previous tutorial a degree=2
quadrature rule (i.e. 3 Gauss points)
will be used. In the end, the total number of Gauss points in the mesh is retrieved
as it will be required to instantiate MFront objects (note that it can be obtained
from the dimension of the Quadrature function spaces or computed from the number of
mesh cells and the chosen quadrature degree):
deg_u = 2
deg_stress = 2
stress_strain_dim = 4
V = VectorFunctionSpace(mesh, "CG", deg_u)
# Quadrature space for sigma
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress,
dim=stress_strain_dim, quad_scheme='default')
W = FunctionSpace(mesh, We)
# Quadrature space for p
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
# Quadrature space for tangent matrix
Wce = TensorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress,
shape=(stress_strain_dim, stress_strain_dim),
quad_scheme='default')
WC = FunctionSpace(mesh, Wce)
# get total number of gauss points
ngauss = W0.dim()
Various functions are defined to keep track of the current internal state (stresses, current strain estimate, cumulative plastic strain and cpnsistent tangent matrix) as well as the previous displacement, current displacement estimate and current iteration correction:
# Define functions based on Quadrature spaces
sig = Function(W, name="Current stresses")
sig_old = Function(W)
Eps1 = Function(W, name="Current strain increment estimate at the end of the end step")
Ct = Function(WC, name="Consistent tangent operator")
p = Function(W0, name="Cumulative plastic strain")
p_old = Function(W0)
u = Function(V, name="Displacement at the beginning of the time step")
u1 = Function(V, name="Current displacement estimate at the end of the end step")
du = Function(V, name="Iteration correction")
v = TrialFunction(V)
u_ = TestFunction(V)
Material constitutive law definition using MFront¶
We now define the material. First let us make a rapid description of some classes and functions introduced by the MFrontGenericInterface project that will be helpful for this tutorial:
the
Behaviour
class handles all the information about a specific MFront behaviour. It is created by theload
function which takes the path to a library, the name of a behaviour and a modelling hypothesis.the
MaterialDataManager
class handles a bunch of integration points. It is instantiated using an instance of theBehaviour
class and the number of integration points [1]. TheMaterialDataManager
contains various interesting members:s0
: data structure of theMaterialStateManager
type which stands for the material state at the beginning of the time step.s1
: data structure of theMaterialStateManager
type which stands for the material state at the end of the time step.K
: anumpy
array containing the consistent tangent operator at each integration points.
the
MaterialStateManager
class describe the state of a material. The following members will be useful in the following:gradients
: a numpy array containing the value of the gradients at each integration points. The number of components of the gradients at each integration points is given by thegradients_stride
member.thermodynamic_forces
: a numpy array containing the value of the thermodynamic forces at each integration points. The number of components of the thermodynamic forces at each integration points is given by thethermodynamic_forces_stride
member.internal_state_variables
: a numpy array containing the value of the internal state variables at each integration points. The number of internal state variables at each integration points is given by theinternal_state_variables_stride
member.
the
setMaterialProperty
andsetExternalStateVariable
functions can be used to set the value a material property or a state variable respectively.the
update
function updates an instance of theMaterialStateManager
by copying the states1
at the end of the time step in the states0
at the beginning of the time step.the
revert
function reverts an instance of theMaterialStateManager
by copying the states0
at the beginning of the time step in the states1
at the end of the time step.the
integrate
function triggers the behaviour integration at each integration points. Various overloads of this function exist. We will use a version taking as argument aMaterialStateManager
, the time increment and a range of integration points.
In the present case, we compute a plane strain von Mises plasticity with isotropic
linear hardening. The material behaviour is implemented in the IsotropicLinearHardeningPlasticity.mfront
file
which must first be compiled to generate the appropriate librairies as follows:
> mfront --obuild --interface=generic IsotropicLinearHardeningPlasticity.mfront
We can then setup the MaterialDataManager
:
# Defining the modelling hypothesis
h = mgis_bv.Hypothesis.PlaneStrain
# Loading the behaviour
b = mgis_bv.load('src/libBehaviour.so','IsotropicLinearHardeningPlasticity',h)
# Setting the material data manager
m = mgis_bv.MaterialDataManager(b, ngauss)
# elastic parameters
E = 70e3
nu = 0.3
# yield strength
sig0 = 250.
Et = E/100.
# hardening slope
H = E*Et/(E-Et)
for s in [m.s0, m.s1]:
mgis_bv.setMaterialProperty(s, "YoungModulus", E)
mgis_bv.setMaterialProperty(s, "PoissonRatio", nu)
mgis_bv.setMaterialProperty(s, "HardeningSlope", H)
mgis_bv.setMaterialProperty(s, "YieldStrength", sig0)
mgis_bv.setExternalStateVariable(s, "Temperature", 293.15)
Boundary conditions and external loading are defined as before along with the analytical limit load solution:
bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]
n = FacetNormal(mesh)
q_lim = float(2/sqrt(3)*ln(Re/Ri)*sig0)
loading = Expression("-q*t", q=q_lim, t=0, degree=2)
def F_ext(v):
return loading*dot(n, v)*ds(4)
Global problem and Newton-Raphson procedure¶
Before writing the global Newton system, we first define the strain measure
function eps_MFront
consistent with the MFront conventions (see
http://tfel.sourceforge.net/tensors.html for details). We also define the custom
quadrature measure and the projection function onto Quadrature spaces:
def eps_MFront(v):
e = sym(grad(v))
return as_vector([e[0, 0], e[1, 1], 0, sqrt(2)*e[0, 1]])
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
def local_project(v, V, u=None):
"""
projects v on V with custom quadrature scheme dedicated to
FunctionSpaces V of `Quadrature` type
if u is provided, result is appended to u
"""
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
The bilinear form of the global problem is obtained using the consistent tangent
matrix Ct
and the MFront strain measure, whereas the right-hand side consists of
the residual between the internal forces associated with the current
stress state sig
and the external force vector.
a_Newton = inner(eps_MFront(v), dot(Ct, eps_MFront(u_)))*dxm
res = -inner(eps_MFront(u_), sig)*dxm + F_ext(u_)
Before defining the Newton-Raphson loop, we set up the output file and appropriate FunctionSpace (here piecewise constant) and Function for output of the equivalent plastic strain since XDMF output does not handle Quadrature elements:
file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
The tangent stiffness is also initialized with the elasticity matrix:
it = mgis_bv.IntegrationType.PredictionWithElasticOperator
mgis_bv.integrate(m, it, 0, 0, m.n);
Ct.vector().set_local(m.K.flatten())
Ct.vector().apply("insert")
The main difference with respect to the pure FEniCS implementation of the previous
tutorial is that MFront computes the current stress state and stiffness matrix
(integrate
method) based on the value of the total strain which is computed
from the total displacement estimate u1
. The associated strain is projected
onto the appropriate Quadrature function space so that its array of values at all
Gauss points is given to MFront via the m.s1.gradients
attribute. The flattened
array of stress and tangent stiffness values are then used to update the current
stress and tangent stiffness variables. The cumulated plastic strain is also
retrieved from the internal_state_variables
attribute (\(p\) being the last
column in the present case). At the end of the iteration loop, the material
behaviour and the previous displacement variable are updated:
Nitermax, tol = 200, 1e-8 # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):
loading.t = t
A, Res = assemble_system(a_Newton, res, bc)
nRes0 = Res.norm("l2")
nRes = nRes0
u1.assign(u)
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "mumps")
# update the current estimate of the displacement at the end of the time step
u1.assign(u1+du)
# compute the current estimate of the strain at the end of the
# time step using `MFront` conventions
local_project(eps_MFront(u1), W, Eps1)
# copy the strain values to `MGIS`
m.s1.gradients[:, :] = Eps1.vector().get_local().reshape((m.n, stress_strain_dim))
# integrate the behaviour
it = mgis_bv.IntegrationType.IntegrationWithConsistentTangentOperator
mgis_bv.integrate(m, it, 0, 0, m.n);
# getting the stress and consistent tangent operator back to
# the FEniCS world.
sig.vector().set_local(m.s1.thermodynamic_forces.flatten())
sig.vector().apply("insert")
Ct.vector().set_local(m.K.flatten())
Ct.vector().apply("insert")
# retrieve cumulated plastic strain values
p.vector().set_local(m.s1.internal_state_variables[:, -1])
p.vector().apply("insert")
# solve Newton system
A, Res = assemble_system(a_Newton, res, bc)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
# update the displacement for the next increment
u.assign(u1)
# update the material
mgis_bv.update(m)
# postprocessing results
file_results.write(u, t)
p_avg.assign(project(p, P0))
file_results.write(p_avg, t)
results[i+1, :] = (u(Ri, 0)[0], t)
import matplotlib.pyplot as plt
plt.plot(results[:, 0], results[:, 1], "-o")
plt.xlabel("Displacement of inner boundary")
plt.ylabel(r"Applied pressure $q/q_{lim}$")
plt.show()
Note
Note that we defined the cumulative plastic strain variable \(p\) in FEniCS only for post-processing purposes. In fact, FEniCS deals only with the global equilibrium whereas MFront manages the history of internal state variables, so that this variable would not have been needed if we were not interested in post-processing it.
Results and future works¶
We can verify that the convergence of the Newton-Raphson procedure is extremely similar between the MFront-based implementation and the pure FEniCS one, the same number of iterations per increment is obtained along with close values of the residual.
Total computing time took approximately:
- 5.9s for the present MFront implementation against
- 6.8s for the previous FEniCS-only implementation
Several points need to be mentioned regarding this implementation efficiency:
- MGIS can handle parallel integration of the constitutive law which has not been used for the present computation
- the present approach can be improved by letting MGIS reuse the memory already allocated by FEniCS which will reduce information transfer times and memory consumption
- extension to large strains is a work in progress
- this FEniCS/MGIS coupling will make it possible, in a near future, to test in a rapid manner generalized constitutive laws (higher-order theories, phase-field) and/or multiphysics couplings
References¶
[HEL2015] | Helfer, Thomas, Bruno Michel, Jean-Michel Proix, Maxime Salvo, Jérôme Sercombe, and Michel Casella. 2015. Introducing the Open-Source Mfront Code Generator: Application to Mechanical Behaviours and Material Knowledge Management Within the PLEIADES Fuel Element Modelling Platform. Computers & Mathematics with Applications. <https://doi.org/10.1016/j.camwa.2015.06.027>. |
[1] | Note that an instance of MaterialDataManager keeps a reference to the behaviour |
which has been used for its initialization: the user must ensure that this behaviour outlives the instance of the MaterialDataManager, otherwise memory corruption may occur.