plasticity_mfront.py.rst 19 KB
Newer Older
1

2 3 4 5 6 7 8 9 10 11
.. raw:: html

 <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><p align="center"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/88x31.png"/></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">Creative Commons Attribution-ShareAlike 4.0 International License</a></p>

.. _PlasticityMFront:

======================================================================
Elasto-plastic analysis implemented using the `MFront` code generator
======================================================================

12 13
This numerical tour has been written in collaboration with **Thomas Helfer** (thomas.helfer@cea.fr), MFront's main developper.

14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
-------------
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 <https://bitbucket.org/fenics-apps/fenics-solid-mechanics>`_
is discussed.

The considered example is exactly the same as the pure FEniCS tutorial :ref:`vonMisesPlasticity` 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 
:download:`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`
38 39 40 41 42
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
43 44 45 46

-------------
Prerequisites
-------------
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
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 :download:`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.
64

65 66 67 68
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
69 70 71 72 73 74

-----------------
Problem position
-----------------

The present implementation will heavily rely on ``Quadrature`` elements to represent
75
previous and current stress states, internal state variables (cumulated
76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367
equivalent plastic strain :math:`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 :math:`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 the ``load`` 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 the ``Behaviour`` class and the
  number of integration points [#]_. The ``MaterialDataManager`` contains
  various interesting members:
  
  - ``s0``: data structure of the ``MaterialStateManager`` type which stands for the material state at the beginning of the time step.
  - ``s1``: data structure of the ``MaterialStateManager`` type which stands for the material state at the end of the time step.
  - ``K``: a ``numpy`` 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 the
      ``gradients_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 the ``thermodynamic_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 the
      ``internal_state_variables_stride`` member.
      
* the ``setMaterialProperty`` and ``setExternalStateVariable`` functions can
  be used to set the value a material property or a state variable
  respectively.
* the ``update`` function updates an instance of the
  ``MaterialStateManager`` by copying the state ``s1`` at the end of the
  time step in the state ``s0`` at the beginning of the time step.
* the ``revert`` function reverts an instance of the
  ``MaterialStateManager`` by copying the state ``s0`` at the beginning of
  the time step in the state ``s1`` 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 a ``MaterialStateManager``, 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 :download:`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 (:math:`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.IntegrationWithConsitentTangentOperator
         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)
368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389
     
     # 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 :math:`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
-------------------------
390 391

We can verify that the convergence of the Newton-Raphson procedure is extremely similar
392
between the `MFront`-based implementation and the pure FEniCS one, the same number of 
393 394 395 396
iterations per increment is obtained along with close values of the residual.

**Total computing time** took approximately:

397
* 5.9s for the present `MFront` implementation against
398 399
* 6.8s for the previous FEniCS-only implementation
  
400 401 402 403 404 405 406 407 408 409 410
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
   
411 412 413 414 415 416 417 418 419 420 421 422 423 424
------------
 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>.

.. [#] 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.