{ "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 }