{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "$\\newcommand{\\sym}{\\operatorname{sym}}\n", "\\newcommand{\\skew}{\\operatorname{skew}}\n", "\\renewcommand{\\tr}{\\operatorname{tr}}\n", "\\newcommand{\\T}{^\\text{T}}\n", "\\newcommand{\\utilde}[1]{\\tilde{#1}}\n", "\\newcommand{\\uttilde}[1]{\\tilde{\\tilde{#1}}}\n", "\\newcommand{\\tgrad}{\\utilde{\\nabla}}\n", "\\newcommand{\\bX}{\\boldsymbol{X}}\n", "\\newcommand{\\bsig}{\\boldsymbol{\\sigma}}\n", "\\newcommand{\\btheta}{\\boldsymbol{\\theta}}\n", "\\newcommand{\\bvarepsilon}{\\boldsymbol{\\varepsilon}}\n", "\\newcommand{\\bepsilon}{\\boldsymbol{\\epsilon}}\n", "\\newcommand{\\bgamma}{\\boldsymbol{\\gamma}}\n", "\\newcommand{\\bchi}{\\boldsymbol{\\kappa}}\n", "\\newcommand{\\bX}{\\boldsymbol{X}}\n", "\\newcommand{\\bxi}{\\boldsymbol{\\xi}}\n", "\\newcommand{\\ba}{\\boldsymbol{a}}\n", "\\newcommand{\\be}{\\boldsymbol{e}}\n", "\\newcommand{\\bt}{\\boldsymbol{t}}\n", "\\newcommand{\\bu}{\\boldsymbol{u}}\n", "\\newcommand{\\bv}{\\boldsymbol{v}}\n", "\\newcommand{\\bx}{\\boldsymbol{x}}\n", "\\newcommand{\\bF}{\\boldsymbol{F}}\n", "\\newcommand{\\bI}{\\boldsymbol{I}}\n", "\\newcommand{\\bM}{\\boldsymbol{M}}\n", "\\newcommand{\\bN}{\\boldsymbol{N}}\n", "\\newcommand{\\bQ}{\\boldsymbol{Q}}\n", "\\newcommand{\\bR}{\\boldsymbol{R}}\n", "\\newcommand{\\bV}{\\boldsymbol{V}}\n", "\\newcommand{\\be}{\\boldsymbol{e}}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear shell model\n", "\n", "This tour is dedicated to the implementation of small-strain elastic shell model. In the [fenics-shells project](https://fenics-shells.readthedocs.io), the shell surface is described by a mapping between a parametric 2D domain and the shell initial stress-free configuration. This approach makes it impossible to tackle non-manifold surfaces or closed surfaces for which the mapping would exhibit singularities.\n", "\n", "In the present implementation, we will instead work directly with the shell surface (either manifold or non-manifold) described as an assembly of planar facets. In this respect, the approach will bear important similarities with the [Elastic 3D beam structures](https://comet-fenics.readthedocs.io/en/latest/demo/beams_3D/beams_3D.html). In particular, we will take advantage of `UFL` ability to compute gradients over a manifold.\n", "\n", "Regarding the interpolation spaces, we will use a `P2/CR` mixed space for displacements and rotations, i.e. the same choice as in the [Locking-free Reissner-Mindlin plate with Crouzeix-Raviart interpolation](https://comet-fenics.readthedocs.io/en/latest/demo/reissner_mindlin_crouzeix_raviart/reissner_mindlin_CR.html) tour. The present shell model will therefore be a small-displacement version of the nonlinear shell model developed in [[CAM03]](#References).\n", "\n", "> The implementation of the fully nonlinear shell model of [[CAM03]](#References) is developed in a separate project for research purposes. You can directly contact the author for more details on this matter.\n", "\n", "For the mesh generation, we refer to the [Generating a shell model for a shallow arch I-beam with the Gmsh Python API](I_beam_gmsh.ipynb) tour." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Shell kinematics and 3D linearized strain\n", "\n", "The shell initial reference configuration consists of piecewise flat portions and $\\bxi$ in $\\mathbb{R}^3$ will denote an initial point on this surface. $(\\be_1,\\be_2,\\be_3)$ will be a reference local frame with $\\be_3$ the normal to the shell surface, this frame being orthonormal and piecewise constant. A material point $\\bX$ in the shell reference configuration will then be given by:\n", "\\begin{equation}\n", "\\bX(\\bxi,\\zeta) = \\bxi + \\zeta \\be_3 \\quad \\text{with } \\zeta \\in [-h/2;h/2]\n", "\\end{equation}\n", "where $h$ is the shell thickness\n", "\n", "The shell kinematics will be described by its mid-surface displacement $\\bu$ and the infinitesimal rotation vector $\\btheta$ of its director. The new normal director is $\\ba_3 = \\bR(\\btheta)\\be_3 = \\be_3+\\btheta\\times\\be_3$ with $\\bR$ the infinitesimal rotation matrix associated with $\\btheta$. Neglecting any thickness change in the shell kinematics, the material point $\\bx$ in the deformed configuration associated with $\\bX$ will then be given by:\n", "\n", "\\begin{equation}\n", "\\bx = \\bxi + \\bu(\\bxi) + \\zeta\\ba_3 = \\bxi + \\bu(\\bxi) + \\zeta(\\be_3+\\btheta(\\bxi)\\times\\be_3)\n", "\\end{equation}\n", "Differentiating with respect to $\\bX$, we get:\n", "\n", "\\begin{align}\n", "d\\bx &= d\\bxi + \\nabla\\bu\\cdot d\\bxi + \\zeta(\\nabla\\btheta\\cdot d\\bxi)\\times\\be_3 + d\\zeta(\\be_3+\\btheta\\times\\be_3) \\\\\n", "&= d\\bX + \\nabla\\bu\\cdot d\\bxi - \\zeta (\\be_3\\times \\nabla\\btheta )\\cdot d\\bxi - d\\zeta (\\be_3\\times\\btheta) + \\text{ h.o.t.} \\notag\n", "\\end{align}\n", "\n", "where we retained only up to first order terms in $\\bu,\\btheta$. The 3D deformation gradient is then given by:\n", "\\begin{equation}\n", "\\bF = \\dfrac{d\\bx}{d\\bX} = \\bI + \\tgrad\\bu - \\zeta \\be_3\\times \\tgrad\\btheta - (\\be_3\\times\\btheta)\\otimes\\be_3 \n", "\\end{equation} \n", "\n", "where we introduced the in-plane gradient (i.e. the gradient with respect to the shell local tangent plane $(\\be_1,\\be_2)$) as follows $\\tgrad\\bv = \\partial_1\\bv \\otimes \\be_1+ \\partial_2\\bv \\otimes \\be_2$. More generally, we will use the following notation $\\utilde{\\bv} = v_1\\be_1+v_2\\be_2$ to denote the in-plane part of $\\bv$.\n", "\n", "The linearized strain tensor is then given by:\n", "\\begin{equation}\n", "\\bvarepsilon = \\sym\\left(\\tgrad\\bu - \\zeta \\be_3\\times \\tgrad\\btheta - (\\be_3\\times\\btheta)\\otimes\\be_3 \\right) + \\epsilon(\\zeta) \\be_3\\otimes\\be_3\n", "\\end{equation} \n", "where $\\sym$ denotes the symmetrized part of a tensor and where we added an incompatible out-of-plane strain $\\epsilon(\\zeta)$ which will enable to enforce a plane-stress constraint. The 3D strain can be split into its in-plane components $\\utilde{\\utilde{\\bvarepsilon}}$, out-of-plane shear components $\\bgamma = 2\\utilde{\\bvarepsilon}_3$ and out-of-plane transverse components $\\varepsilon_{33}$:\n", "\\begin{align}\n", "\\uttilde{\\bvarepsilon} &= \\sym\\left(\\tgrad\\utilde{\\bu} - \\zeta \\be_3\\times \\tgrad\\btheta \\right)= \\boldsymbol{\\epsilon} - \\zeta \\bchi \\\\\n", "\\bgamma &= \\tgrad u_3 - \\be_3\\times\\btheta \\tag{1}\\\\\n", "\\varepsilon_{33} &= \\epsilon(\\zeta)\n", "\\end{align}\n", "where $\\bepsilon=\\sym(\\tgrad\\utilde{\\bu})$ is the membrane strain and $\\bchi = \\sym(\\be_3\\times \\tgrad\\btheta)$ the bending curvature. We see that we recover the definition of the curvature and shear strain of the Reissner-Mindlin plate model." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Elastic energy density\n", "\n", "\n", "The internal work of deformation density per shell unit surface is then given by:\n", "\n", "\\begin{align}\n", "w_\\text{def} &= \\int_{-h/2}^{h/2} \\bsig:\\bvarepsilon d\\zeta \\\\\n", "&= \\left(\\int_{-h/2}^{h/2} \\uttilde{\\bsig} d\\zeta\\right):\\bepsilon + \\left(\\int_{-h/2}^{h/2} (-\\zeta)\\uttilde{\\bsig} d\\zeta\\right):\\bchi \\notag\\\\\n", "&\\phantom{=}+ \\left(\\int_{-h/2}^{h/2} \\utilde{\\bsig}_3 d\\zeta\\right)\\cdot\\bgamma + \\int_{-h/2}^{h/2} \\sigma_{33}\\epsilon(\\zeta)d\\zeta \\notag \\\\\n", "&= \\bN:\\bepsilon + \\bM:\\bchi + \\bQ\\cdot\\bgamma + \\int_{-h/2}^{h/2} \\sigma_{33}\\epsilon(\\zeta)d\\zeta \\notag\n", "\\end{align}\n", "\n", "where $\\bN$ is the shell membrane tensor, $\\bM$ the bending moment tensor and $\\bQ$ the shear force vector appearing in duality with $\\bepsilon,\\bchi,\\bgamma$ respectively. The out-of-plane stress $\\sigma_{33}$ appears in duality with the out-of-plane strain $\\epsilon(\\zeta)$. The latter is in fact a purely local variable which can be locally eliminated to enforce the plane stress condition $\\sigma_{33}=0$.\n", "\n", "Indeed, the shell free energy density can be defined as:\n", "\n", "\\begin{align}\n", "\\Psi &= \\int_{-h/2}^{h/2} \\psi(\\bvarepsilon) d\\zeta \\\\\n", "&= \\int_{-h/2}^{h/2} \\frac{1}{2}(\\lambda (\\tr\\bvarepsilon)^2 + 2\\mu \\bvarepsilon:\\bvarepsilon) d\\zeta\n", "\\end{align}\n", "\n", "for isotropic linear elasticity.\n", "\n", "When minimizing this energy over the shell degrees of freedom, the minimization over $\\epsilon$ is in fact local and we can replace the above free-energy by its plane-stress counterpart:\n", "\n", "\n", "\\begin{equation}\n", "\\begin{array}{rl}\\Psi_\\text{ps} &= \\displaystyle{\\int_{-h/2}^{h/2} \\inf_{\\epsilon(z)}\\psi(\\bvarepsilon) d\\zeta}\\\\\n", "& =\\displaystyle{\\int_{-h/2}^{h/2}\\left(\\psi_\\text{ps}(\\uttilde{\\bvarepsilon}) + \\dfrac{1}{2}\\mu \\bgamma^2\\right) d\\zeta }\n", "\\end{array}\n", "\\end{equation}\n", "\n", "where $\\psi_\\text{ps}$ is the elastic strain energy density for a 2D plane-stress material:\n", "\n", "\\begin{equation}\n", "\\psi_\\text{ps}(\\uttilde{\\bvarepsilon}) = \\frac{1}{2}\\left(\\lambda_\\text{ps} (\\tr\\uttilde{\\bvarepsilon})^2 + 2\\mu \\uttilde{\\bvarepsilon}:\\uttilde{\\bvarepsilon}\\right) \n", "\\end{equation}\n", "\n", "with the plane-stress pseudo Lamé coefficient (see the [2D linear elasticity](https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/2D_elasticity.py.html) tour):\n", "\n", "\\begin{equation}\n", "\\lambda_\\text{ps} = \\dfrac{2\\lambda\\mu}{\\lambda+2\\mu}\n", "\\end{equation}\n", "\n", "Taking the derivative of the strain energy with respect to the 3D strain components gives the local 3D stress-strain constitutive relation:\n", "\n", "\\begin{align}\n", "\\uttilde{\\bsig} &= \\lambda_\\text{ps} \\tr(\\uttilde{\\bvarepsilon})\\mathbf{1} + 2\\mu \\bvarepsilon = \\mathbb{C}_\\text{ps}:\\uttilde{\\bvarepsilon} \\\\\n", "\\utilde{\\bsig}_{3} &= \\mu \\bgamma \\\\\n", "\\sigma_{33} &= 0\n", "\\end{align}\n", "\n", "where $\\mathbb{C}_\\text{ps}$ is the plane-stress elasticity tensor." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stress-resultant constitutive equations\n", "\n", "Plugging relations (1) into the previous 3D constitutive relation yields the following stress-resultant equations, assuming a homogeneous shell across the thickness:\n", "\n", "\\begin{align}\n", "\\bN &= \\int_{-h/2}^{h/2} \\uttilde{\\bsig} d\\zeta = h \\int_{-h/2}^{h/2} \\mathbb{C}_\\text{ps}:(\\boldsymbol{\\epsilon} - \\zeta \\bchi) d\\zeta = h\\mathbb{C}_\\text{ps}:\\boldsymbol{\\epsilon} \\\\\n", "\\bM &= \\int_{-h/2}^{h/2} (-\\zeta)\\uttilde{\\bsig} d\\zeta = h \\int_{-h/2}^{h/2} (-\\zeta)\\mathbb{C}_\\text{ps}:(\\boldsymbol{\\epsilon} - \\zeta \\bchi) d\\zeta = \\dfrac{h^3}{12}\\mathbb{C}_\\text{ps}:\\boldsymbol{\\bchi}\\\\\n", "\\bQ &= \\int_{-h/2}^{h/2} \\utilde{\\bsig}_3 d \\zeta = \\mu h \\bgamma\n", "\\end{align}\n", "\n", "> Note that with the present purely kinematic approach, no shear correction factor appears for the shear force constitutive relation.\n", "\n", "## FEniCS implementation\n", "\n", "### Import statements and physical parameters\n", "\n", "We first load the relevant packages and some UFL functions. We then define the corressponding physical parameters. In particular, we consider here a uniformly distributed loading of unit intensity along the downward global $Z$ direction and of intensity 0.2 along the $Y$ direction perpendicular to the arch plane in order to induce a combined bending and twisting deflection mode." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "from dolfin import *\n", "from ufl import Jacobian, replace\n", "import numpy as np\n", "import meshio\n", "\n", "# material parameters\n", "thick = Constant(1e-3)\n", "E = Constant(210e9)\n", "nu = Constant(0.3)\n", "lmbda = E * nu / (1 + nu) / (1 - 2 * nu)\n", "mu = E / 2 / (1 + nu)\n", "lmbda_ps = 2 * lmbda * mu / (lmbda + 2 * mu)\n", "\n", "# loading (self-weight)\n", "f = Constant((0, 0.2, -1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Loading the mesh and computing a local tangent frame\n", "\n", "We now load the mesh as well as the facet MeshFunction from the XDMF files generated from a Gmsh script using Gmsh's Python API (see [this notebook](I_beam_gmsh.ipynb)).\n", "\n", "We then define the `local_frame` function which returns the UFL representation of the local frame $(\\be_1,\\be_2,\\be_3)$. We first use UFL's `Jacobian` function to retrieve the Jacobian of the mapping from reference cells to spatial coordinates. Contrary to the `fenics-shells` approach which uses a global analytical mapping between the reference domain and the shell surface, we use here the reference element mapping of the finite-element method. For a shell, the Jacobian is of size $3\\times 2$ with both columns $\\bt_1,\\bt_2$ being a vector of the local tangent plane. We therefore compute $\\be_3$, the normal to the mid-surface from $\\bt_1\\times \\bt_2$. We then have to choose a convention to compute $\\be_1,\\be_2$, knowing $\\be_3$. Our convention is that $\\be_1$ lies in the plane orthogonal to $\\be_Y$ and $\\be_3$. If $\\be_3$ happens to be colinear to $\\be_Y$, we set $\\be_1=\\be_Z$. $\\be_2$ then follows from $\\be_2=\\be_3\\times\\be_1$." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "filename = \"I_beam.xdmf\"\n", "\n", "mesh = Mesh()\n", "with XDMFFile(filename) as mesh_file:\n", " mesh_file.read(mesh)\n", "\n", "\n", "mvc = MeshValueCollection(\"size_t\", mesh, 1)\n", "with XDMFFile(filename.replace(\".\", \"_facet_region.\")) as infile:\n", " infile.read(mvc, \"name_to_read\")\n", "facets = MeshFunction(\"size_t\", mesh, mvc)\n", "\n", "def local_frame(mesh):\n", " t = Jacobian(mesh)\n", " if mesh.geometric_dimension() == 2:\n", " t1 = as_vector([t[0, 0], t[1, 0], 0])\n", " t2 = as_vector([t[0, 1], t[1, 1], 0])\n", " else:\n", " t1 = as_vector([t[0, 0], t[1, 0], t[2, 0]])\n", " t2 = as_vector([t[0, 1], t[1, 1], t[2, 1]])\n", " e3 = cross(t1, t2)\n", " e3 /= sqrt(dot(e3, e3))\n", " ey = as_vector([0, 1, 0])\n", " ez = as_vector([0, 0, 1])\n", " e1 = cross(ey, e3)\n", " norm_e1 = sqrt(dot(e1, e1))\n", " e1 = conditional(lt(norm_e1, 0.5), ez, e1 / norm_e1)\n", "\n", " e2 = cross(e3, e1)\n", " e2 /= sqrt(dot(e2, e2))\n", " return e1, e2, e3\n", "\n", "\n", "frame = local_frame(mesh)\n", "e1, e2, e3 = frame\n", "\n", "VT = VectorFunctionSpace(mesh, \"DG\", 0, dim=3)\n", "\n", "ffile = XDMFFile(\"output.xdmf\")\n", "ffile.parameters[\"functions_share_mesh\"] = True\n", "for (i, ei) in enumerate(frame):\n", " ei = Function(VT, name=\"e{}\".format(i + 1))\n", " ei.assign(project(frame[i], VT))\n", " ffile.write(ei, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The local frame is then stored in an output XDMFFile. Note that when it is represented via `DG0` Functions but Paraview does not seem to handle very well such fields on a surface shell. The fields seem to be displayed as `CG1` Functions, inducing erratic orientations near the boundary. Overall, the orientation is however consistent with what has been defined ($\\be_1$ in blue, $\\be_2$ in green, $\\be_3$ in red).\n", "\n", "