Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License

The corresponding files can be obtained from:

Linear Buckling Analysis of a 3D solid

This demo has been written in collaboration with Eder Medina (Harvard University).

In this numerical tour, we will demonstrate how to compute the buckling modes of a three-dimensional elastic solid under a conservative loading. 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 and the Eulerian buckling of a beam demos. For more details on the theoretical formulation, the reader can refer to [NGU00].

\[\newcommand{\bq}{\boldsymbol{q}} \newcommand{\bu}{\boldsymbol{u}} \newcommand{\bv}{\boldsymbol{v}} \newcommand{\bt}{\boldsymbol{t}} \newcommand{\be}{\boldsymbol{e}} \newcommand{\beps}{\boldsymbol{\varepsilon}} \newcommand{\bsig}{\boldsymbol{\sigma}} \newcommand{\T}{^{\text{T}}}\]
../../_images/buckling_modes.png

Geometrically Nonlinear Elasticity

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.

We consider the Green-Lagrange strain tensor \(\be\) defined as:

\[e_{ij} = \frac{1}{2} (u_{i,j} + u_{j,i}) + \frac{1}{2} u_{k,i}u_{k,j}\]

Note that this can be rewritten in tensor format as

\[\be = \beps(\bu) + \dfrac{1}{2}\bq(\bu, \bu)\]

where we isolate the linear \(\beps(\bu)\) and nonlinear quadratic part \(\bq(\bu, \bu)=\nabla\bu\T\nabla\bu\).

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:

\[\psi(\be) = \dfrac{1}{2}\be:\mathbb{C}: \be\]

where \(\mathbb{C}\) is the elastic moduli tensor.

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)\).

Equilibrium and stability conditions

Let us consider the total potential energy functional:

\[\mathcal{E}_\text{pot}(\bu;\lambda) = \int_\Omega \psi(\be(\bu)) d\Omega - \int_{\Gamma_N} \bt(\lambda) \cdot \bu dS\]

The equilibrium equations are obtained from the first variation which must vanish:

\[\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\]

where \(\delta\bu\) is any kinematically admissible perturbation direction and where:

\[\delta\be[\delta\bu] = \beps(\delta\bu) + \nabla\bu\T\nabla\delta\bu = \beps(\delta\bu) + \bq(\bu,\delta\bu)\]

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:

\[\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\]

where:

\[\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)\]

Bifurcation point

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:

\[\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\]

Linear buckling analysis

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:

\[\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\]

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 :

\[\int_\Omega (\beps(\delta\bv):\mathbb{C}:\beps(\delta\bu) + \lambda^c(\nabla \delta \bv) \bsig_0 (\nabla \delta\bu)\T)d\Omega = 0\]

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:

\[[\mathbf{K}]+\lambda[\mathbf{K}_G(\bu_0)]=0\]

where \([\mathbf{K}]\) is the classical linear elastic stiffness matrix and \([\mathbf{K}_G(\bu_0)]\) the so-called geometrical stiffness matrix.

FEniCS implementation

Defining the Geometry

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\).

[1]:

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

L, b, h = 1, 0.01, 0.03
Nx, Ny, Nz = 51, 5, 5

mesh = BoxMesh(Point(0, 0, 0), Point(L, b, h), Nx, Ny, Nz)

Solving for the pre-stressed state \(\bu_0\)

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.

[2]:

E , nu = 1e3, 0.

mu = Constant(E/2/(1+nu))
lmbda = Constant(E*nu/(1+nu)/(1-2*nu))

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(v.geometric_dimension())+2*mu*eps(v)

# Compute the linearized unit preload
N0 = 1
T = Constant((-N0, 0, 0))
V = VectorFunctionSpace(mesh, 'Lagrange', degree = 2)
v = TestFunction(V)
du = TrialFunction(V)

# Two different ways to define how to impose boundary conditions
def left(x,on_boundary):
    return near(x[0], 0.)

def right(x, on_boundary):
    return near(x[0], L)

boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
AutoSubDomain(right).mark(boundary_markers, 1)

# Clamped Boundary and Simply Supported Boundary Conditions
bcs = [DirichletBC(V, Constant((0, 0, 0)), left),
       DirichletBC(V.sub(1), Constant(0), right),
       DirichletBC(V.sub(2), Constant(0), right)]

# Linear Elasticity Bilinear Form
a = inner(sigma(du), eps(v))*dx

# External loading on the right end
l = dot(T, v)*ds(1)

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.

[3]:

K = PETScMatrix()
f = PETScVector()
assemble_system(a, l, bcs, A_tensor=K, b_tensor=f)

u = Function(V,name = "Displacement")
solve(K, u.vector(), f, "mumps")

# Output the trivial solution
ffile = XDMFFile("output/solution.xdmf")
ffile.parameters["functions_share_mesh"] = True
ffile.parameters["flush_output"] = True
ffile.write(u, 0)

Forming the geometric stiffness matrix and the buckling eigenvalue problem

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:

[4]:

kgform = -inner(sigma(u), grad(du).T*grad(v))*dx

KG = PETScMatrix()
assemble(kgform, KG)

# Zero out the rows to enforce boundary condition
for bc in bcs:
    bc.zero(KG)

The buckling points of interest are the smallest critical loads, i.e. the points corresponding to the smallest \(\lambda^c\).

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 tour.

[5]:

# What Solver Configurations exist?
eigensolver = SLEPcEigenSolver(K, KG)
print(eigensolver.parameters.str(True))

<Parameter set "slepc_eigenvalue_solver" containing 8 parameter(s) and parameter set(s)>

  slepc_eigenvalue_solver  |    type    value    range  access  change
  --------------------------------------------------------------------
  maximum_iterations       |     int  <unset>  Not set       0       0
  problem_type             |  string  <unset>  Not set       0       0
  solver                   |  string  <unset>  Not set       0       0
  spectral_shift           |  double  <unset>  Not set       0       0
  spectral_transform       |  string  <unset>  Not set       0       0
  spectrum                 |  string  <unset>  Not set       0       0
  tolerance                |  double  <unset>  Not set       0       0
  verbose                  |    bool  <unset>  Not set       0       0

[6]:

eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['solver'] = "krylov-schur"
eigensolver.parameters['tolerance'] = 1.e-12
eigensolver.parameters['spectral_transform']  = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 1e-2

# Request the smallest 3 eigenvalues
N_eig = 3
print("Computing {} first eigenvalues...".format(N_eig))
eigensolver.solve(N_eig)
print("Number of converged eigenvalues:", eigensolver.get_number_converged())

eigenmode = Function(V)

# Analytical beam theory solution
# F_cr,n = alpha_n^2*EI/L^2
# where alpha_n are solutions to tan(alpha) = alpha
alpha = np.array([1.4303, 2.4509, 3.4709] + [(n+1/2) for n in range(4, 10)])*pi
I = b**3*h/12
S = b*h

for i in range(eigensolver.get_number_converged()):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)

    critical_load_an = alpha[i]**2*float(E*I/N0/S)/L**2
    print("Critical load factor {}: {:.5f} FE | {:.5f} Beam".format(i+1, r, critical_load_an))

    # Initialize function and assign eigenvector
    eigenmode.vector()[:] = rx
    eigenmode.rename("Eigenmode {}".format(i+1), "")

    ffile.write(eigenmode, 0)

Computing 3 first eigenvalues...
Number of converged eigenvalues: 3
Critical load factor 1: 0.16821 FE | 0.16826 Beam
Critical load factor 2: 0.49691 FE | 0.49405 Beam
Critical load factor 3: 0.98918 FE | 0.99084 Beam
Note that the analysis above is not limited to problems subject to Neumann boundary conditions. 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.

References

[NGU00]: Nguyen, Q. S. (2000). Stability and nonlinear solid mechanics. Wiley.