{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Linear Buckling Analysis of a 3D solid\n", "\n", "This demo has been written in collaboration with [Eder Medina (Harvard University)](mailto:e_medina@g.harvard.edu).\n", "\n", "In this numerical tour, we will demonstrate how to compute the buckling modes of a three-dimensional elastic solid under a conservative loading. \n", "The critical buckling loads and buckling modes are computed by solving a generalized non-hermitian eigenvalue problem using the `SLEPcEigensolver`. This example is closely related to the [Modal analysis of an elastic structure](https://comet-fenics.readthedocs.io/en/latest/demo/modal_analysis_dynamics/cantilever_modal.py.html) and the [Eulerian buckling of a beam](../beam_buckling/beam_buckling.ipynb) demos. For more details on the theoretical formulation, the reader can refer to [[NGU00]](#References).\n", "\n", "$$\\newcommand{\\bq}{\\boldsymbol{q}}\n", "\\newcommand{\\bu}{\\boldsymbol{u}}\n", "\\newcommand{\\bv}{\\boldsymbol{v}}\n", "\\newcommand{\\bt}{\\boldsymbol{t}}\n", "\\newcommand{\\be}{\\boldsymbol{e}}\n", "\\newcommand{\\beps}{\\boldsymbol{\\varepsilon}}\n", "\\newcommand{\\bsig}{\\boldsymbol{\\sigma}}\n", "\\newcommand{\\T}{^{\\text{T}}}$$" ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ ".. figure:: buckling_modes.png\n", " :scale: 50%\n", " :align: center" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Geometrically Nonlinear Elasticity\n", "\n", "To formulate the buckling problem, we must consider a general geometrically nonlinear framework. More precisely, we will consider the case of small-strain/large displacements that are expected in typical buckling analysis. \n", "\n", "We consider the Green-Lagrange strain tensor $\\be$ defined as:\n", "\n", "$$e_{ij} = \\frac{1}{2} (u_{i,j} + u_{j,i}) + \\frac{1}{2} u_{k,i}u_{k,j}$$\n", "\n", "Note that this can be rewritten in tensor format as\n", "\n", "$$ \\be = \\beps(\\bu) + \\dfrac{1}{2}\\bq(\\bu, \\bu) $$\n", "\n", "where we isolate the linear $\\beps(\\bu)$ and nonlinear quadratic part $\\bq(\\bu, \\bu)=\\nabla\\bu\\T\\nabla\\bu$.\n", "\n", "Since we restrict to a small strain assumption, we adopt a standard linear Hookean behaviour for the material material model which can be described by the following quadratic elastic energy density:\n", "\n", "$$ \\psi(\\be) = \\dfrac{1}{2}\\be:\\mathbb{C}: \\be $$\n", "\n", "where $\\mathbb{C}$ is the elastic moduli tensor.\n", "\n", "We seek to find the underlying equilibrium displacement field $\\bu$ as a function of a scalar loading parameter $\\lambda$ under a set of prescribed external surface forces $\\bt(\\lambda)$. \n", "\n", "## Equilibrium and stability conditions\n", "Let us consider the total potential energy functional:\n", "\n", "$$ \\mathcal{E}_\\text{pot}(\\bu;\\lambda) = \\int_\\Omega \\psi(\\be(\\bu)) d\\Omega - \\int_{\\Gamma_N} \\bt(\\lambda) \\cdot \\bu dS $$\n", "\n", "The equilibrium equations are obtained from the first variation which must vanish: \n", "\n", "$$ \\partial_u\\mathcal{E}_\\text{pot}[\\delta\\bu] = \\int_\\Omega \\be(\\bu):\\mathbb{C}:\\delta\\be[\\delta\\bu] d\\Omega - \\int_{\\Gamma_N} \\bt(\\lambda) \\cdot \\delta\\bu dS $$\n", "\n", "where $\\delta\\bu$ is any kinematically admissible perturbation direction and where:\n", "\n", "$$\n", "\\delta\\be[\\delta\\bu] = \\beps(\\delta\\bu) + \\nabla\\bu\\T\\nabla\\delta\\bu = \\beps(\\delta\\bu) + \\bq(\\bu,\\delta\\bu)\n", "$$\n", "\n", "The stability conditions of an equilibrium are obtained from the Lejeune-Dirichlet theorem stating that the energy second variation must be positive. Here the second variation bilinear form is given by:\n", "\n", "$$ \\partial_{uu}\\mathcal{E}_\\text{pot}[\\delta\\bu,\\delta \\bv] = \\int_\\Omega (\\delta \\be[\\delta\\bv]:\\mathbb{C}:\\delta\\be[\\delta\\bu] + \\be:\\mathbb{C}:\\delta^2\\be[\\delta\\bu,\\delta\\bv] )d\\Omega $$\n", "\n", "where:\n", "$$\n", "\\delta^2\\be[\\delta\\bu,\\delta\\bv] = \\delta\\bq(\\bu,\\delta \\bu)[\\delta\\bv] = \\nabla \\delta\\bv\\T\\nabla\\delta\\bu = \\bq(\\delta\\bv,\\delta\\bu)\n", "$$\n", "\n", "## Bifurcation point\n", "\n", "A given equilibrium solution $\\lambda^c, \\bu^c=\\bu(\\lambda^c)$ is a bifurcation point if the second variation vanishes for some $\\delta\\bu$. $\\lambda^c$ defines the bifurcation load and $\\delta\\bu$ the bifurcation mode. Hence, we have:\n", "$$\n", "\\int_\\Omega ((\\beps(\\delta\\bv) + \\bq(\\bu^c,\\delta\\bv)):\\mathbb{C}:(\\beps(\\delta\\bu) + \\bq(\\bu^c,\\delta\\bu)) + \\be(\\bu^c):\\mathbb{C}:\\bq(\\delta\\bv,\\delta\\bu))d\\Omega = 0\n", "$$\n", "\n", "## Linear buckling analysis\n", "\n", "When performing a linear buckling analysis, we look for bifurcation points on the fundamental equilibrium branch obtained by a small displacement assumption. As a result, in the above bifurcation equation the quadratic contribution in the first term can be neglected and the Green-Lagrange strain $\\be(\\bu^c)$ can be replaced by the linearized strain $\\beps(\\bu^c)$. Besides, in the small-displacement assumption, if the loading depends linearly on the load parameter $\\lambda$, the small displacement solution also depends linearly on it. We therefore have $\\bu^c= \\lambda^c\\bu_0$ and the bifurcation condition becomes:\n", "$$\n", "\\int_\\Omega (\\beps(\\delta\\bv):\\mathbb{C}:\\beps(\\delta\\bu) + \\lambda^c\\beps(\\bu_0):\\mathbb{C}:\\bq(\\delta\\bv,\\delta\\bu))d\\Omega = 0\n", "$$\n", "\n", "After, introducing the inital pre-stress $\\bsig_0 = \\mathbb{C}:\\beps(\\bu_0)$, the last term can be rewritten as $\\bsig_0:\\bq(\\delta\\bv,\\delta\\bu) = (\\nabla \\delta \\bv) \\bsig_0 (\\nabla \\delta\\bu)\\T$ so that :\n", "$$\n", "\\int_\\Omega (\\beps(\\delta\\bv):\\mathbb{C}:\\beps(\\delta\\bu) + \\lambda^c(\\nabla \\delta \\bv) \\bsig_0 (\\nabla \\delta\\bu)\\T)d\\Omega = 0\n", "$$\n", "\n", "We recognize here a linear eigenvalue problem where the first term corresponds to the classical linear elasticity bilinear form and the second term is an additional contribution depending on the pre-stressed state $\\bu_0$. Transforming this variational problem into matrix notations yields:\n", "$$\n", "[\\mathbf{K}]+\\lambda[\\mathbf{K}_G(\\bu_0)]=0\n", "$$\n", "\n", "where $[\\mathbf{K}]$ is the classical linear elastic stiffness matrix and $[\\mathbf{K}_G(\\bu_0)]$ the so-called *geometrical stiffness* matrix." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## FEniCS implementation\n", "\n", "### Defining the Geometry\n", "\n", "Here we will consider a fully three dimensional beam-like structure with length $L = 1$ and of rectangular cross-section with height $h = 0.01$, and width $b = 0.03$. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from dolfin import *\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "L, b, h = 1, 0.01, 0.03\n", "Nx, Ny, Nz = 51, 5, 5\n", "\n", "mesh = BoxMesh(Point(0, 0, 0), Point(L, b, h), Nx, Ny, Nz)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving for the pre-stressed state $\\bu_0$\n", "\n", "We first need to compute the pre-stressed linear state $\\bu_0$. It is obtained by solving a simple linear elasticity problem. The material is assumed here isotropic and we consider a pre-stressed state obtained from the application of a unit compression applied at the beam right end in the $X$ direction while the $Y$ and $Z$ displacement are fixed, mimicking a simple support condition. The beam is fully clamped on its left end." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "E , nu = 1e3, 0.\n", "\n", "mu = Constant(E/2/(1+nu))\n", "lmbda = Constant(E*nu/(1+nu)/(1-2*nu))\n", "\n", "def eps(v):\n", " return sym(grad(v))\n", "def sigma(v):\n", " return lmbda*tr(eps(v))*Identity(v.geometric_dimension())+2*mu*eps(v)\n", "\n", "# Compute the linearized unit preload \n", "N0 = 1\n", "T = Constant((-N0, 0, 0))\n", "V = VectorFunctionSpace(mesh, 'Lagrange', degree = 2)\n", "v = TestFunction(V)\n", "du = TrialFunction(V)\n", "\n", "# Two different ways to define how to impose boundary conditions\n", "def left(x,on_boundary):\n", " return near(x[0], 0.)\n", "\n", "def right(x, on_boundary):\n", " return near(x[0], L)\n", "\n", "boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)\n", "ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)\n", "AutoSubDomain(right).mark(boundary_markers, 1)\n", "\n", "# Clamped Boundary and Simply Supported Boundary Conditions\n", "bcs = [DirichletBC(V, Constant((0, 0, 0)), left),\n", " DirichletBC(V.sub(1), Constant(0), right),\n", " DirichletBC(V.sub(2), Constant(0), right)]\n", "\n", "# Linear Elasticity Bilinear Form\n", "a = inner(sigma(du), eps(v))*dx\n", "\n", "# External loading on the right end\n", "l = dot(T, v)*ds(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will reuse our stiffness matrix in the linear buckling computation and hence we supply a `PETScMatrix` and `PETScVector` to the system assembler. We then use the `solve` function on these objects and specify the linear solver method.\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "K = PETScMatrix()\n", "f = PETScVector()\n", "assemble_system(a, l, bcs, A_tensor=K, b_tensor=f)\n", "\n", "u = Function(V,name = \"Displacement\")\n", "solve(K, u.vector(), f, \"mumps\")\n", "\n", "# Output the trivial solution \n", "ffile = XDMFFile(\"output/solution.xdmf\")\n", "ffile.parameters[\"functions_share_mesh\"] = True\n", "ffile.parameters[\"flush_output\"] = True\n", "ffile.write(u, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Forming the geometric stiffness matrix and the buckling eigenvalue problem\n", "\n", "From the previous solution, the prestressed state $\\bsig_0$ entering the geometric stiffness matrix expression is simply given by `sigma(u)`. After having formed the geometric stiffness matrix, we will call the `SLEPcEigenSolver` for solving a generalized eigenvalue problem $\\mathbf{Ax}=\\lambda \\mathbf{Bx}$ with here $\\mathbf{A}=\\mathbf{K}$ and $\\mathbf{B}=-\\mathbf{K}_G$. We therefore include directly the negative sign in the definition of the geometric stiffness form:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "kgform = -inner(sigma(u), grad(du).T*grad(v))*dx\n", "\n", "KG = PETScMatrix()\n", "assemble(kgform, KG)\n", "\n", "# Zero out the rows to enforce boundary condition\n", "for bc in bcs:\n", " bc.zero(KG)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The buckling points of interest are the smallest critical loads, i.e. the points corresponding to the smallest $\\lambda^c$. \n", "\n", "We now pass some parameters to the SLEPc solver. In particular, we peform a shift-invert transform as discussed in the [Eulerian buckling of an elastic beam](https://comet-fenics.readthedocs.io/en/latest/demo/beam_buckling/beam_buckling.html) tour." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", " slepc_eigenvalue_solver | type value range access change\n", " --------------------------------------------------------------------\n", " maximum_iterations | int Not set 0 0\n", " problem_type | string Not set 0 0\n", " solver | string Not set 0 0\n", " spectral_shift | double Not set 0 0\n", " spectral_transform | string Not set 0 0\n", " spectrum | string Not set 0 0\n", " tolerance | double Not set 0 0\n", " verbose | bool Not set 0 0\n" ] } ], "source": [ "# What Solver Configurations exist?\n", "eigensolver = SLEPcEigenSolver(K, KG)\n", "print(eigensolver.parameters.str(True))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Computing 3 first eigenvalues...\n", "Number of converged eigenvalues: 3\n", "Critical load factor 1: 0.16821 FE | 0.16826 Beam\n", "Critical load factor 2: 0.49691 FE | 0.49405 Beam\n", "Critical load factor 3: 0.98918 FE | 0.99084 Beam\n" ] } ], "source": [ "eigensolver.parameters['problem_type'] = 'gen_hermitian'\n", "eigensolver.parameters['solver'] = \"krylov-schur\"\n", "eigensolver.parameters['tolerance'] = 1.e-12\n", "eigensolver.parameters['spectral_transform'] = 'shift-and-invert'\n", "eigensolver.parameters['spectral_shift'] = 1e-2\n", "\n", "# Request the smallest 3 eigenvalues\n", "N_eig = 3 \n", "print(\"Computing {} first eigenvalues...\".format(N_eig))\n", "eigensolver.solve(N_eig)\n", "print(\"Number of converged eigenvalues:\", eigensolver.get_number_converged())\n", "\n", "eigenmode = Function(V)\n", "\n", "# Analytical beam theory solution\n", "# F_cr,n = alpha_n^2*EI/L^2\n", "# where alpha_n are solutions to tan(alpha) = alpha\n", "alpha = np.array([1.4303, 2.4509, 3.4709] + [(n+1/2) for n in range(4, 10)])*pi\n", "I = b**3*h/12\n", "S = b*h\n", "\n", "for i in range(eigensolver.get_number_converged()):\n", " # Extract eigenpair\n", " r, c, rx, cx = eigensolver.get_eigenpair(i)\n", " \n", " critical_load_an = alpha[i]**2*float(E*I/N0/S)/L**2\n", " print(\"Critical load factor {}: {:.5f} FE | {:.5f} Beam\".format(i+1, r, critical_load_an))\n", "\n", " # Initialize function and assign eigenvector\n", " eigenmode.vector()[:] = rx\n", " eigenmode.rename(\"Eigenmode {}\".format(i+1), \"\")\n", "\n", " ffile.write(eigenmode, 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> Note that the analysis above is not limited to problems subject to Neumann boundary conditions. \n", "One can equivalenty compute the prestress resulting from a displacement control simulation and the critical buckling load will correspond to the critical buckling displacement of the structure." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "[NGU00]: Nguyen, Q. S. (2000). *Stability and nonlinear solid mechanics*. Wiley." ] } ], "metadata": { "celltoolbar": "Format de la Cellule Texte Brut", "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.8.10" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }