Incompressible elasticity
This example is also available as a Jupyter notebook: incompressible_elasticity.ipynb
.
Introduction
Mixed elements can be used to overcome locking when the material becomes incompressible. However, for an element to be stable, it needs to fulfill the LBB condition. In this example we will consider two different element formulations
- linear displacement with linear pressure approximation (does not fulfill LBB)
- quadratic displacement with linear pressure approximation (does fulfill LBB)
The quadratic/linear element is also known as the Taylor-Hood element. We will consider Cook's Membrane with an applied traction on the right hand side.
Commented program
What follows is a program spliced with comments. The full program, without comments, can be found in the next section.
using Ferrite, Tensors
using BlockArrays, SparseArrays, LinearAlgebra
First we generate a simple grid, specifying the 4 corners of Cooks membrane.
function create_cook_grid(nx, ny)
corners = [
Vec{2}((0.0, 0.0)),
Vec{2}((48.0, 44.0)),
Vec{2}((48.0, 60.0)),
Vec{2}((0.0, 44.0)),
]
grid = generate_grid(Triangle, (nx, ny), corners)
# facesets for boundary conditions
addfacetset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0)
addfacetset!(grid, "traction", x -> norm(x[1]) ≈ 48.0)
return grid
end;
Next we define a function to set up our cell- and FacetValues.
function create_values(interpolation_u, interpolation_p)
# quadrature rules
qr = QuadratureRule{RefTriangle}(3)
facet_qr = FacetQuadratureRule{RefTriangle}(3)
# cell and FacetValues for u
cellvalues_u = CellValues(qr, interpolation_u)
facetvalues_u = FacetValues(facet_qr, interpolation_u)
# cellvalues for p
cellvalues_p = CellValues(qr, interpolation_p)
return cellvalues_u, cellvalues_p, facetvalues_u
end;
We create a DofHandler, with two fields, :u
and :p
, with possibly different interpolations
function create_dofhandler(grid, ipu, ipp)
dh = DofHandler(grid)
add!(dh, :u, ipu) # displacement
add!(dh, :p, ipp) # pressure
close!(dh)
return dh
end;
We also need to add Dirichlet boundary conditions on the "clamped"
facetset. We specify a homogeneous Dirichlet bc on the displacement field, :u
.
function create_bc(dh)
dbc = ConstraintHandler(dh)
add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "clamped"), x -> zero(x), [1, 2]))
close!(dbc)
return dbc
end;
The material is linear elastic, which is here specified by the shear and bulk moduli
struct LinearElasticity{T}
G::T
K::T
end
Now to the assembling of the stiffness matrix. This mixed formulation leads to a blocked element matrix. Since Ferrite does not force us to use any particular matrix type we will use a BlockedArray
from BlockArrays.jl
.
function doassemble(
cellvalues_u::CellValues,
cellvalues_p::CellValues,
facetvalues_u::FacetValues,
K::SparseMatrixCSC, grid::Grid, dh::DofHandler, mp::LinearElasticity
)
f = zeros(ndofs(dh))
assembler = start_assemble(K, f)
nu = getnbasefunctions(cellvalues_u)
np = getnbasefunctions(cellvalues_p)
fe = BlockedArray(zeros(nu + np), [nu, np]) # local force vector
ke = BlockedArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix
# traction vector
t = Vec{2}((0.0, 1 / 16))
# cache ɛdev outside the element routine to avoid some unnecessary allocations
ɛdev = [zero(SymmetricTensor{2, 2}) for i in 1:getnbasefunctions(cellvalues_u)]
for cell in CellIterator(dh)
fill!(ke, 0)
fill!(fe, 0)
assemble_up!(ke, fe, cell, cellvalues_u, cellvalues_p, facetvalues_u, grid, mp, ɛdev, t)
assemble!(assembler, celldofs(cell), ke, fe)
end
return K, f
end;
The element routine integrates the local stiffness and force vector for all elements. Since the problem results in a symmetric matrix we choose to only assemble the lower part, and then symmetrize it after the loop over the quadrature points.
function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facetvalues_u, grid, mp, ɛdev, t)
n_basefuncs_u = getnbasefunctions(cellvalues_u)
n_basefuncs_p = getnbasefunctions(cellvalues_p)
u▄, p▄ = 1, 2
reinit!(cellvalues_u, cell)
reinit!(cellvalues_p, cell)
# We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
for q_point in 1:getnquadpoints(cellvalues_u)
for i in 1:n_basefuncs_u
ɛdev[i] = dev(symmetric(shape_gradient(cellvalues_u, q_point, i)))
end
dΩ = getdetJdV(cellvalues_u, q_point)
for i in 1:n_basefuncs_u
divδu = shape_divergence(cellvalues_u, q_point, i)
δu = shape_value(cellvalues_u, q_point, i)
for j in 1:i
Ke[BlockIndex((u▄, u▄), (i, j))] += 2 * mp.G * ɛdev[i] ⊡ ɛdev[j] * dΩ
end
end
for i in 1:n_basefuncs_p
δp = shape_value(cellvalues_p, q_point, i)
for j in 1:n_basefuncs_u
divδu = shape_divergence(cellvalues_u, q_point, j)
Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ
end
for j in 1:i
p = shape_value(cellvalues_p, q_point, j)
Ke[BlockIndex((p▄, p▄), (i, j))] += - 1 / mp.K * δp * p * dΩ
end
end
end
symmetrize_lower!(Ke)
# We integrate the Neumann boundary using the FacetValues.
# We loop over all the facets in the cell, then check if the facet
# is in our `"traction"` facetset.
for facet in 1:nfacets(cell)
if (cellid(cell), facet) ∈ getfacetset(grid, "traction")
reinit!(facetvalues_u, cell, facet)
for q_point in 1:getnquadpoints(facetvalues_u)
dΓ = getdetJdV(facetvalues_u, q_point)
for i in 1:n_basefuncs_u
δu = shape_value(facetvalues_u, q_point, i)
fe[i] += (δu ⋅ t) * dΓ
end
end
end
end
return
end
function symmetrize_lower!(Ke)
for i in 1:size(Ke, 1)
for j in (i + 1):size(Ke, 1)
Ke[i, j] = Ke[j, i]
end
end
return
end;
To evaluate the stresses after solving the problem we once again loop over the cells in the grid. Stresses are evaluated in the quadrature points, however, for export/visualization you typically want values in the nodes of the mesh, or as single data points per cell. For the former you can project the quadrature point data to a finite element space (see the example with the L2Projector
in Post processing and visualization). In this example we choose to compute the mean value of the stress within each cell, and thus end up with one data point per cell. The mean value is computed as
\[\bar{\boldsymbol{\sigma}}_i = \frac{1}{ |\Omega_i|} \int_{\Omega_i} \boldsymbol{\sigma}\, \mathrm{d}\Omega, \quad |\Omega_i| = \int_{\Omega_i} 1\, \mathrm{d}\Omega\]
where $\Omega_i$ is the domain occupied by cell number $i$, and $|\Omega_i|$ the volume (area) of the cell. The integrals are evaluated using numerical quadrature with the help of cellvalues for u and p, just like in the assembly procedure.
Note that even though all strain components in the out-of-plane direction are zero (plane strain) the stress components are not. Specifically, $\sigma_{33}$ will be non-zero in this formulation. Therefore we expand the strain to a 3D tensor, and then compute the (3D) stress tensor.
function compute_stresses(
cellvalues_u::CellValues, cellvalues_p::CellValues,
dh::DofHandler, mp::LinearElasticity, a::Vector
)
ae = zeros(ndofs_per_cell(dh)) # local solution vector
u_range = dof_range(dh, :u) # local range of dofs corresponding to u
p_range = dof_range(dh, :p) # local range of dofs corresponding to p
# Allocate storage for the stresses
σ = zeros(SymmetricTensor{2, 3}, getncells(dh.grid))
# Loop over the cells and compute the cell-average stress
for cc in CellIterator(dh)
# Update cellvalues
reinit!(cellvalues_u, cc)
reinit!(cellvalues_p, cc)
# Extract the cell local part of the solution
for (i, I) in pairs(celldofs(cc))
ae[i] = a[I]
end
# Loop over the quadrature points
σΩi = zero(SymmetricTensor{2, 3}) # stress integrated over the cell
Ωi = 0.0 # cell volume (area)
for qp in 1:getnquadpoints(cellvalues_u)
dΩ = getdetJdV(cellvalues_u, qp)
# Evaluate the strain and the pressure
ε = function_symmetric_gradient(cellvalues_u, qp, ae, u_range)
p = function_value(cellvalues_p, qp, ae, p_range)
# Expand strain to 3D
ε3D = SymmetricTensor{2, 3}((i, j) -> i < 3 && j < 3 ? ε[i, j] : 0.0)
# Compute the stress in this quadrature point
σqp = 2 * mp.G * dev(ε3D) - one(ε3D) * p
σΩi += σqp * dΩ
Ωi += dΩ
end
# Store the value
σ[cellid(cc)] = σΩi / Ωi
end
return σ
end;
Now we have constructed all the necessary components, we just need a function to put it all together.
function solve(ν, interpolation_u, interpolation_p)
# material
Emod = 1.0
Gmod = Emod / 2(1 + ν)
Kmod = Emod * ν / ((1 + ν) * (1 - 2ν))
mp = LinearElasticity(Gmod, Kmod)
# Grid, dofhandler, boundary condition
n = 50
grid = create_cook_grid(n, n)
dh = create_dofhandler(grid, interpolation_u, interpolation_p)
dbc = create_bc(dh)
# CellValues
cellvalues_u, cellvalues_p, facetvalues_u = create_values(interpolation_u, interpolation_p)
# Assembly and solve
K = allocate_matrix(dh)
K, f = doassemble(cellvalues_u, cellvalues_p, facetvalues_u, K, grid, dh, mp)
apply!(K, f, dbc)
u = K \ f
# Compute the stress
σ = compute_stresses(cellvalues_u, cellvalues_p, dh, mp, u)
σvM = map(x -> √(3 / 2 * dev(x) ⊡ dev(x)), σ) # von Mise effective stress
# Export the solution and the stress
filename = "cook_" *
(interpolation_u == Lagrange{RefTriangle, 1}()^2 ? "linear" : "quadratic") *
"_linear"
VTKGridFile(filename, grid) do vtk
write_solution(vtk, dh, u)
for i in 1:3, j in 1:3
σij = [x[i, j] for x in σ]
write_cell_data(vtk, σij, "sigma_$(i)$(j)")
end
write_cell_data(vtk, σvM, "sigma von Mises")
end
return u
end
We now define the interpolation for displacement and pressure. We use (scalar) Lagrange interpolation as a basis for both, and for the displacement, which is a vector, we vectorize it to 2 dimensions such that we obtain vector shape functions (and 2nd order tensors for the gradients).
linear_p = Lagrange{RefTriangle, 1}()
linear_u = Lagrange{RefTriangle, 1}()^2
quadratic_u = Lagrange{RefTriangle, 2}()^2
All that is left is to solve the problem. We choose a value of Poissons ratio that results in incompressibility ($ν = 0.5$) and thus expect the linear/linear approximation to return garbage, and the quadratic/linear approximation to be stable.
u1 = solve(0.5, linear_u, linear_p);
u2 = solve(0.5, quadratic_u, linear_p);
Plain program
Here follows a version of the program without any comments. The file is also available here: incompressible_elasticity.jl
.
using Ferrite, Tensors
using BlockArrays, SparseArrays, LinearAlgebra
function create_cook_grid(nx, ny)
corners = [
Vec{2}((0.0, 0.0)),
Vec{2}((48.0, 44.0)),
Vec{2}((48.0, 60.0)),
Vec{2}((0.0, 44.0)),
]
grid = generate_grid(Triangle, (nx, ny), corners)
# facesets for boundary conditions
addfacetset!(grid, "clamped", x -> norm(x[1]) ≈ 0.0)
addfacetset!(grid, "traction", x -> norm(x[1]) ≈ 48.0)
return grid
end;
function create_values(interpolation_u, interpolation_p)
# quadrature rules
qr = QuadratureRule{RefTriangle}(3)
facet_qr = FacetQuadratureRule{RefTriangle}(3)
# cell and FacetValues for u
cellvalues_u = CellValues(qr, interpolation_u)
facetvalues_u = FacetValues(facet_qr, interpolation_u)
# cellvalues for p
cellvalues_p = CellValues(qr, interpolation_p)
return cellvalues_u, cellvalues_p, facetvalues_u
end;
function create_dofhandler(grid, ipu, ipp)
dh = DofHandler(grid)
add!(dh, :u, ipu) # displacement
add!(dh, :p, ipp) # pressure
close!(dh)
return dh
end;
function create_bc(dh)
dbc = ConstraintHandler(dh)
add!(dbc, Dirichlet(:u, getfacetset(dh.grid, "clamped"), x -> zero(x), [1, 2]))
close!(dbc)
return dbc
end;
struct LinearElasticity{T}
G::T
K::T
end
function doassemble(
cellvalues_u::CellValues,
cellvalues_p::CellValues,
facetvalues_u::FacetValues,
K::SparseMatrixCSC, grid::Grid, dh::DofHandler, mp::LinearElasticity
)
f = zeros(ndofs(dh))
assembler = start_assemble(K, f)
nu = getnbasefunctions(cellvalues_u)
np = getnbasefunctions(cellvalues_p)
fe = BlockedArray(zeros(nu + np), [nu, np]) # local force vector
ke = BlockedArray(zeros(nu + np, nu + np), [nu, np], [nu, np]) # local stiffness matrix
# traction vector
t = Vec{2}((0.0, 1 / 16))
# cache ɛdev outside the element routine to avoid some unnecessary allocations
ɛdev = [zero(SymmetricTensor{2, 2}) for i in 1:getnbasefunctions(cellvalues_u)]
for cell in CellIterator(dh)
fill!(ke, 0)
fill!(fe, 0)
assemble_up!(ke, fe, cell, cellvalues_u, cellvalues_p, facetvalues_u, grid, mp, ɛdev, t)
assemble!(assembler, celldofs(cell), ke, fe)
end
return K, f
end;
function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facetvalues_u, grid, mp, ɛdev, t)
n_basefuncs_u = getnbasefunctions(cellvalues_u)
n_basefuncs_p = getnbasefunctions(cellvalues_p)
u▄, p▄ = 1, 2
reinit!(cellvalues_u, cell)
reinit!(cellvalues_p, cell)
# We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
for q_point in 1:getnquadpoints(cellvalues_u)
for i in 1:n_basefuncs_u
ɛdev[i] = dev(symmetric(shape_gradient(cellvalues_u, q_point, i)))
end
dΩ = getdetJdV(cellvalues_u, q_point)
for i in 1:n_basefuncs_u
divδu = shape_divergence(cellvalues_u, q_point, i)
δu = shape_value(cellvalues_u, q_point, i)
for j in 1:i
Ke[BlockIndex((u▄, u▄), (i, j))] += 2 * mp.G * ɛdev[i] ⊡ ɛdev[j] * dΩ
end
end
for i in 1:n_basefuncs_p
δp = shape_value(cellvalues_p, q_point, i)
for j in 1:n_basefuncs_u
divδu = shape_divergence(cellvalues_u, q_point, j)
Ke[BlockIndex((p▄, u▄), (i, j))] += -δp * divδu * dΩ
end
for j in 1:i
p = shape_value(cellvalues_p, q_point, j)
Ke[BlockIndex((p▄, p▄), (i, j))] += - 1 / mp.K * δp * p * dΩ
end
end
end
symmetrize_lower!(Ke)
# We integrate the Neumann boundary using the FacetValues.
# We loop over all the facets in the cell, then check if the facet
# is in our `"traction"` facetset.
for facet in 1:nfacets(cell)
if (cellid(cell), facet) ∈ getfacetset(grid, "traction")
reinit!(facetvalues_u, cell, facet)
for q_point in 1:getnquadpoints(facetvalues_u)
dΓ = getdetJdV(facetvalues_u, q_point)
for i in 1:n_basefuncs_u
δu = shape_value(facetvalues_u, q_point, i)
fe[i] += (δu ⋅ t) * dΓ
end
end
end
end
return
end
function symmetrize_lower!(Ke)
for i in 1:size(Ke, 1)
for j in (i + 1):size(Ke, 1)
Ke[i, j] = Ke[j, i]
end
end
return
end;
function compute_stresses(
cellvalues_u::CellValues, cellvalues_p::CellValues,
dh::DofHandler, mp::LinearElasticity, a::Vector
)
ae = zeros(ndofs_per_cell(dh)) # local solution vector
u_range = dof_range(dh, :u) # local range of dofs corresponding to u
p_range = dof_range(dh, :p) # local range of dofs corresponding to p
# Allocate storage for the stresses
σ = zeros(SymmetricTensor{2, 3}, getncells(dh.grid))
# Loop over the cells and compute the cell-average stress
for cc in CellIterator(dh)
# Update cellvalues
reinit!(cellvalues_u, cc)
reinit!(cellvalues_p, cc)
# Extract the cell local part of the solution
for (i, I) in pairs(celldofs(cc))
ae[i] = a[I]
end
# Loop over the quadrature points
σΩi = zero(SymmetricTensor{2, 3}) # stress integrated over the cell
Ωi = 0.0 # cell volume (area)
for qp in 1:getnquadpoints(cellvalues_u)
dΩ = getdetJdV(cellvalues_u, qp)
# Evaluate the strain and the pressure
ε = function_symmetric_gradient(cellvalues_u, qp, ae, u_range)
p = function_value(cellvalues_p, qp, ae, p_range)
# Expand strain to 3D
ε3D = SymmetricTensor{2, 3}((i, j) -> i < 3 && j < 3 ? ε[i, j] : 0.0)
# Compute the stress in this quadrature point
σqp = 2 * mp.G * dev(ε3D) - one(ε3D) * p
σΩi += σqp * dΩ
Ωi += dΩ
end
# Store the value
σ[cellid(cc)] = σΩi / Ωi
end
return σ
end;
function solve(ν, interpolation_u, interpolation_p)
# material
Emod = 1.0
Gmod = Emod / 2(1 + ν)
Kmod = Emod * ν / ((1 + ν) * (1 - 2ν))
mp = LinearElasticity(Gmod, Kmod)
# Grid, dofhandler, boundary condition
n = 50
grid = create_cook_grid(n, n)
dh = create_dofhandler(grid, interpolation_u, interpolation_p)
dbc = create_bc(dh)
# CellValues
cellvalues_u, cellvalues_p, facetvalues_u = create_values(interpolation_u, interpolation_p)
# Assembly and solve
K = allocate_matrix(dh)
K, f = doassemble(cellvalues_u, cellvalues_p, facetvalues_u, K, grid, dh, mp)
apply!(K, f, dbc)
u = K \ f
# Compute the stress
σ = compute_stresses(cellvalues_u, cellvalues_p, dh, mp, u)
σvM = map(x -> √(3 / 2 * dev(x) ⊡ dev(x)), σ) # von Mise effective stress
# Export the solution and the stress
filename = "cook_" *
(interpolation_u == Lagrange{RefTriangle, 1}()^2 ? "linear" : "quadratic") *
"_linear"
VTKGridFile(filename, grid) do vtk
write_solution(vtk, dh, u)
for i in 1:3, j in 1:3
σij = [x[i, j] for x in σ]
write_cell_data(vtk, σij, "sigma_$(i)$(j)")
end
write_cell_data(vtk, σvM, "sigma von Mises")
end
return u
end
linear_p = Lagrange{RefTriangle, 1}()
linear_u = Lagrange{RefTriangle, 1}()^2
quadratic_u = Lagrange{RefTriangle, 2}()^2
u1 = solve(0.5, linear_u, linear_p);
u2 = solve(0.5, quadratic_u, linear_p);
This page was generated using Literate.jl.