# 2023-11-06 Solid Mechanics#

## Last time#

• Stabilized methods for transport

• VMS and SUPG

• FE interfaces

## Today#

• Mixed finite elements

• Intro to solid mechanics, Ratel

using Plots
default(linewidth=3)
using LinearAlgebra
using SparseArrays
import NLsolve: nlsolve

function my_spy(A)
cmax = norm(vec(A), Inf)
s = max(1, ceil(120 / size(A, 1)))
spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

function vander_legendre_deriv(x, k=nothing)
if isnothing(k)
k = length(x) # Square by default
end
m = length(x)
Q = ones(m, k)
dQ = zeros(m, k)
Q[:, 2] = x
dQ[:, 2] .= 1
for n in 1:k-2
Q[:, n+2] = ((2*n + 1) * x .* Q[:, n+1] - n * Q[:, n]) / (n + 1)
dQ[:, n+2] = (2*n + 1) * Q[:,n+1] + dQ[:,n]
end
Q, dQ
end

x, _ = gausslobatto(P)
Pk, _ = vander_legendre_deriv(x)
Bp, Dp = vander_legendre_deriv(q, P)
B = Bp / Pk
D = Dp / Pk
x, q, w, B, D
end

function fe1_mesh(P, nelem)
x = LinRange(-1, 1, nelem+1)
rows = Int[]
cols = Int[]
for i in 1:nelem
append!(rows, (i-1)*P+1:i*P)
append!(cols, (i-1)*(P-1)+1:i*(P-1)+1)
end
x, sparse(cols, rows, ones(nelem*P))'
end

function xnodal(x, P)
xn = Float64[]
xref, _ = gausslobatto(P)
for i in 1:length(x)-1
xL, xR = x[i:i+1]
append!(xn, (xL+xR)/2 .+ (xR-xL)/2 * xref[1+(i>1):end])
end
xn
end

struct FESpace
P::Int
Q::Int
nelem::Int
x::Vector
xn::Vector
Et::SparseMatrixCSC{Float64, Int64}
q::Vector
w::Vector
B::Matrix
D::Matrix
x, E = fe1_mesh(P, nelem)
xn = xnodal(x, P)
_, q, w, B, D = febasis(P, Q, quadrature)
new(P, Q, nelem, x, xn, E', q, w, B, D)
end
end

# Extract out what we need for element e
function fe_element(fe, e)
xL, xR = fe.x[e:e+1]
q = (xL+xR)/2 .+ (xR-xL)/2*fe.q
w = (xR - xL)/2 * fe.w
E = fe.Et[:, (e-1)*fe.P+1:e*fe.P]'
dXdx = ones(fe.Q) * 2 / (xR - xL)
q, w, E, dXdx
end

function fe_residual(u_in, fe, fq; bci=[1], bcv=[1.])
u = copy(u_in); v = zero(u)
u[bci] = bcv
for e in 1:fe.nelem
q, w, E, dXdx = fe_element(fe, e)
B, D = fe.B, fe.D
ue = E * u
uq = B * ue
Duq = dXdx .* (D * ue)
f0, f1 = fq(q, uq, Duq)
ve = B' * (w .* f0) + D' * (dXdx .* w .* f1)
v += E' * ve
end
v[bci] = u_in[bci] - u[bci]
#println("residual")
v
end

function fe_jacobian(u_in, fe, dfq; bci=[1], bcv=[1.])
u = copy(u_in); u[bci] = bcv
rows, cols, vals = Int[], Int[], Float64[]
for e in 1:fe.nelem
q, w, E, dXdx = fe_element(fe, e)
B, D, P = fe.B, fe.D, fe.P
ue = E * u
uq = B * ue; Duq = dXdx .* (D * ue)
K = zeros(P, P)
for j in 1:fe.P
du = B[:,j]
Ddu = dXdx .* D[:,j]
df0, df1 = dfq(q, uq, du, Duq, Ddu)
K[:,j] = B' * (w .* df0) + D' * (dXdx .* w .* df1)
end
inds = rowvals(E')
append!(rows, kron(ones(P), inds))
append!(cols, kron(inds, ones(P)))
append!(vals, vec(K))
end
A = sparse(rows, cols, vals)
A[bci, :] .= 0; A[:, bci] .= 0
A[bci,bci] = diagm(ones(length(bci)))
A
end

fe_jacobian (generic function with 1 method)


LFAToolkit.jl

# Finite element interfaces: Deal.II#

$\begin{gather*} v^T F(u) \sim \int_\Omega v \cdot \color{olive}{f_0(u, \nabla u)} + \nabla v \!:\! \color{olive}{f_1(u, \nabla u)} \quad v^T J w \sim \int_\Omega \begin{bmatrix} v \\ \nabla v \end{bmatrix}^T \color{teal}{\begin{bmatrix} f_{0,0} & f_{0,1} \\ f_{1,0} & f_{1,1} \end{bmatrix}} \begin{bmatrix} w \\ \nabla w \end{bmatrix} \\ J w = \sum_e \mathcal E_e^T \underbrace{\begin{bmatrix} B_I \\ B_{\nabla} \end{bmatrix}^T \begin{bmatrix} I & \\ & \left( \frac{\partial X}{\partial x}\right)^T \end{bmatrix}}_{\texttt{fe\_values}} W_q \color{teal}{\begin{bmatrix} f_{0,0} & f_{0,1} \\ f_{1,0} & f_{1,1} \end{bmatrix}} \underbrace{\begin{bmatrix} I & \\ & \left( \frac{\partial X}{\partial x}\right) \end{bmatrix} \begin{bmatrix} B_I \\ B_{\nabla} \end{bmatrix}}_{\texttt{fe\_values}} \mathcal E_e w_L \end{gather*}$
for e in elems:
fe_values.reinit()
for q in q_points:
for i in test_functions:
for j in trial_functions
K_e[i,j] += ...
f_e[i] += ...
for f in e.faces:
if f.at_boundary():
fe_face_values.reinit()
for q in q_points:
...


# Finite element interfaces: MOOSE#

## Materials#

• Can be written without knowledge of finite elements

• Registration allows libraries of materials (some in MOOSE, others packaged separatle)

• Example: crystal plasticity

• Code is C++, so can do dirty things

• table lookups, proprietary code

• implicit materials (Newton solve at each quadrature point)

## Composition in configuration files#

• Select materials from libraries

• Multiphysics composition

• Multiscale coupling

# Periodic table of finite elements#

## Careful choice of mixed elements allows exactly satisfying discrete identities#

$\nabla\cdot\nabla\times u = 0$
$\nabla\times\nabla u = 0$
• Generalized concept: Finite Element Exterior Calculus (FEEC)

• Improved stability or numerical properties

• Higher order of accuracy for quantity of interest despite non-smooth problem

## $$H^1$$ Poisson#

Find $$p$$ such that

$\int \nabla q \cdot \kappa \nabla p = \int q f, \quad \forall q$

## Mixed Poisson#

Find $$\mathbf u, p$$ such that

(55)#\begin{align} \int \mathbf v \cdot \kappa^{-1} \mathbf u - p \nabla\cdot \mathbf v &= 0 & \forall \mathbf v \\ \int -q \nabla\cdot \mathbf u &= -\int q f & \forall q \end{align}

Dirichlet (essential) and Neumann (natural) boundary conditions are swapped, and we get a “saddle point” linear problem \begin{pmatrix} M & B^T \ B & 0 \end{pmatrix} \begin{pmatrix} \mathbf u \ p \end{pmatrix} =

(56)#$\begin{pmatrix} 0 \\ -f \end{pmatrix}$

$$$With appropriate choice of spaces,$$p,q$$can live in a piecewise constant space while attaining second order accuracy in the fluxes$$\mathbf u$.

## Stokes: (slow) incompressible flow#

Find velocity $$\mathbf u$$ and pressure $$p$$ such that

(57)#\begin{align} -\nabla\cdot(\eta \nabla^s \mathbf u - pI) &= 0 \\ \nabla\cdot \mathbf u &= 0 \end{align}

where $$\nabla^s \mathbf u$$ is the symmetric part of the $$3\times 3$$ gradient $$\nabla \mathbf u$$.

Weak form: find $$(\mathbf u, p)$$ such that

$\int \nabla^s \mathbf v \!:\! \eta \!:\! \nabla^s \mathbf u - p \nabla\cdot \mathbf v - q \nabla\cdot \mathbf u = 0, \quad \forall (\mathbf v, q)$

## Inf-sup stability#

For this problem to be well posed, it is necessary that the divergence of velocity spans the pressure space “nicely”. This is quantified by the “inf-sup constant” (aka. Ladyzhenskaya–Babuška–Brezzi constant)

$\beta = \inf_{\lVert p \rVert =1} \sup_{\mathbf u \ne 0} \frac{\int p \nabla\cdot \mathbf u}{\int \nabla^s \mathbf u \!:\! \eta \!:\! \nabla^s \mathbf u}$

The method loses accuracy if $$\beta$$ decays: you want it to stay of order 1 uniformly. Stability can be quantified numerically by solving an eigenvalue problem (1993). This short paper has many figures quantifying stable and unstable elements.

# Solid mechanics: Ratel Theory Guide#

## Material coordinates#

The current configuration $$x$$ is a function of the initial configuration $$X$$. We typically solve for displacement $$u = x - X$$, and define the deformation gradient

$F = \frac{\partial x}{\partial X} = I + \underbrace{\frac{\partial u}{\partial X}}_H$
where $$H$$ is known as displacement gradient.

# Momentum balance in initial configuration#

$\nabla\cdot \mathbf F \mathbf S = 0$

where $$\mathbf F = I + H$$ and $$\mathbf S$$ is the symmetric stress tensor (Second Piola-Kirchhoff tensor).

# Strain measures#

Stress $$\mathbf S$$ must be defined as a function of displacement. A common choice is to define the right Cauchy-Green tensor

$C = F^T F$

This has value $$I$$ for zero strain. A better formulation uses

$E = \frac 1 2 (C - I) = \underbrace{\frac 1 2 (H + H^T + H^T H)}_{\text{stable}}$

## Neo-Hookean model#

Strain energy density \begin{aligned} \psi \left(\mathbf{E} \right) &= \frac{\lambda}{2} \left( \log J \right)^2 - \mu \log J + \frac \mu 2 \left( \operatorname{trace} \mathbf{C} - 3 \right) \ &= \frac{\lambda}{2} \left( \log J \right)^2 - \mu \log J + \mu \operatorname{trace} \mathbf{E}, \end{aligned}

$\mathbf S = \frac{\partial \psi}{\partial \mathbf E}$