{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Finite-strain elastoplasticity within the logarithmic strain framework\n",
"\n",
"This demo is dedicated to the resolution of a finite-strain elastoplastic problem using the logarithmic strain framework proposed in (Miehe et al., 2002). \n",
"\n",
"## Logarithmic strains \n",
"\n",
"This framework expresses constitutive relations between the Hencky strain measure $\\boldsymbol{H} = \\dfrac{1}{2}\\log (\\boldsymbol{F}^T\\cdot\\boldsymbol{F})$ and its dual stress measure $\\boldsymbol{T}$. This approach makes it possible to extend classical small strain constitutive relations to a finite-strain setting. In particular, the total (Hencky) strain can be split **additively** into many contributions (elastic, plastic, thermal, swelling, etc.). Its trace is also linked with the volume change $J=\\exp(\\operatorname{tr}(\\boldsymbol{H}))$. As a result, the deformation gradient $\\boldsymbol{F}$ is used for expressing the Hencky strain $\\boldsymbol{H}$, a small-strain constitutive law is then written for the $(\\boldsymbol{H},\\boldsymbol{T})$-pair and the dual stress $\\boldsymbol{T}$ is then post-processed to an appropriate stress measure such as the Cauchy stress $\\boldsymbol{\\sigma}$ or Piola-Kirchhoff stresses.\n",
"\n",
"## MFront implementation\n",
"\n",
" TODO:\n",
"* Write the gallery example comments and refer to it?\n",
"\n",
"* or use standard brick with comments here ?\n",
" \n",
"\n",
"\n",
"## FEniCS implementation\n",
"\n",
"We define a box mesh representing half of a beam oriented along the $x$-direction. The beam will be fully clamped on its left side and symmetry conditions will be imposed on its right extremity. The loading consists of a uniform self-weight."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"from dolfin import *\n",
"import mfront_wrapper as mf\n",
"import numpy as np\n",
"import ufl\n",
"\n",
"length, width, height = 1., 0.04, 0.1\n",
"nx, ny, nz = 30, 5, 10\n",
"mesh = BoxMesh(Point(0, -width/2, -height/2.), Point(length, width/2, height/2.), nx, ny, nz)\n",
"\n",
"V = VectorFunctionSpace(mesh, \"CG\", 2)\n",
"u = Function(V, name=\"Displacement\")\n",
"\n",
"def left(x, on_boundary):\n",
" return near(x[0], 0) and on_boundary\n",
"def right(x, on_boundary):\n",
" return near(x[0], length) and on_boundary\n",
"\n",
"bc = [DirichletBC(V, Constant((0.,)*3), left),\n",
" DirichletBC(V.sub(0), Constant(0.), right)]\n",
"\n",
"selfweight = Expression((\"0\", \"0\", \"-t*qmax\"), t=0., qmax = 50e6, degree=0)\n",
"\n",
"file_results = XDMFFile(\"results/finite_strain_plasticity.xdmf\")\n",
"file_results.parameters[\"flush_output\"] = True\n",
"file_results.parameters[\"functions_share_mesh\"] = True"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `MFrontNonlinearMaterial` instance is loaded from the MFront `LogarithmicStrainPlasticity` behaviour. This behaviour is a finite-strain behaviour (`material.finite_strain=True`) which relies on a kinematic description using the total deformation gradient $\\boldsymbol{F}$. By default, a MFront behaviour always returns the Cauchy stress as the stress measure after integration. However, the stress variable dual to the deformation gradient is the first Piloa-Kirchhoff (PK1) stress. An internal option of the MGIS interface is therefore used in the finite-strain context to return the PK1 stress as the \"flux\" associated to the \"gradient\" $\\boldsymbol{F}$. Both quantities are non-symmetric tensors, aranged as a 9-dimensional vector in 3D following [MFront conventions on tensors](http://tfel.sourceforge.net/tensors.html)."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"StandardFiniteStrainBehaviour\n",
"F_CAUCHY\n",
"['DeformationGradient'] [9]\n",
"['FirstPiolaKirchhoffStress'] [9]\n"
]
}
],
"source": [
"material = mf.MFrontNonlinearMaterial(\"../materials/src/libBehaviour.so\",\n",
" \"LogarithmicStrainPlasticity\")\n",
"print(material.behaviour.getBehaviourType())\n",
"print(material.behaviour.getKinematic())\n",
"print(material.get_gradient_names(), material.get_gradient_sizes())\n",
"print(material.get_flux_names(), material.get_flux_sizes())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `MFrontNonlinearProblem` instance must therefore register the deformation gradient as `Identity(3)+grad(u)`. This again done automatically since `\"DeformationGradient\"` is a predefined gradient. The following message will be shown upon calling `solve`:\n",
"```\n",
"Automatic registration of 'DeformationGradient' as I + (grad(Displacement)).\n",
"```\n",
"The loading is then defined and, as for the [small-strain elastoplasticity example](small_strain_elastoplasticity.ipynb), state variables include the `ElasticStrain` and `EquivalentPlasticStrain` since the same behaviour is used as in the small-strain case with the only difference that the total strain is now given by the Hencky strain measure. In particular, the `ElasticStrain` is still a symmetric tensor (vector of dimension 6). Note that it has not been explicitly defined as a state variable in the MFront behaviour file since this is done automatically when using the `IsotropicPlasticMisesFlow` parser.\n",
"\n",
"Finally, we setup a few parameters of the Newton non-linear solver."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(6,)\n"
]
}
],
"source": [
"problem = mf.MFrontNonlinearProblem(u, material, bcs=bc)\n",
"problem.set_loading(dot(selfweight, u)*dx)\n",
"\n",
"p = problem.get_state_variable(\"EquivalentPlasticStrain\")\n",
"epsel = problem.get_state_variable(\"ElasticStrain\")\n",
"print(ufl.shape(epsel))\n",
"\n",
"prm = problem.solver.parameters\n",
"prm[\"absolute_tolerance\"] = 1e-6\n",
"prm[\"relative_tolerance\"] = 1e-6\n",
"prm[\"linear_solver\"] = \"mumps\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"During the load incrementation, we monitor the evolution of the vertical downwards displacement at the middle of the right extremity.\n",
"\n",
"This simulation is a bit heavy to run so we suggest running it in parallel:\n",
"```bash\n",
"mpirun -np 16 python3\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Increment 1\n",
"Automatic registration of 'DeformationGradient' as I + (grad(Displacement)).\n",
"\n",
"Automatic registration of 'Temperature' as a Constant value = 293.15.\n",
"\n",
"Increment 2\n",
"Increment 3\n",
"Increment 4\n"
]
}
],
"source": [
"P0 = FunctionSpace(mesh, \"DG\", 0)\n",
"p_avg = Function(P0, name=\"Plastic strain\")\n",
"\n",
"Nincr = 30\n",
"load_steps = np.linspace(0., 1., Nincr+1)\n",
"results = np.zeros((Nincr+1, 3))\n",
"for (i, t) in enumerate(load_steps[1:]):\n",
" selfweight.t = t\n",
" print(\"Increment \", i+1)\n",
" problem.solve(u.vector())\n",
"\n",
" results[i+1, 0] = -u(length, 0, 0)[2]\n",
" results[i+1, 1] = t\n",
"\n",
" file_results.write(u, t)\n",
" p_avg.assign(project(p, P0))\n",
" file_results.write(p_avg, t)\n",
"\n",
"plt.figure()\n",
"plt.plot(results[:, 0], results[:, 1], \"-o\")\n",
"plt.xlabel(r\"Displacement\")\n",
"plt.ylabel(\"Load\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The load-displacement curve exhibits a classical elastoplastic behaviour rapidly followed by a stiffening behaviour due to membrane catenary effects. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.9"
}
},
"nbformat": 4,
"nbformat_minor": 4
}