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


## 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):
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

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)

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