Linear Elasticity

Multi-material structures with finite elements

Run AllReset

This tutorial solves a linear elasticity problem with the finite element method: an elastic bar made of two different materials, clamped on the left and pulled on the right. It is the 1D analogue of the classic multi-material cantilever benchmark, and small enough that we can check every number against the exact solution.

Problem: A bar of length 1 with cross-section $A = 1$ is clamped at $x = 0$ and loaded with a tension force $F$ at $x = 1$ . The left half is stiff steel ( $E_1 = 200$ ), the right half soft aluminum ( $E_2 = 70$ ). Find the displacement $u(x)$ and the stress $\sigma(x)$ .

Governing Equations

Stress is proportional to strain (Hooke's law), and equilibrium without body forces means the stress is divergence-free:

$\sigma = E \, \varepsilon = E \frac{du}{dx}, \qquad \frac{d\sigma}{dx} = 0.$

So the stress must be constant along the bar ( $\sigma = F/A = F$ ), while the strain jumps at the material interface — the soft half stretches more. The exact displacement is piecewise linear:

$u(x) = \begin{cases} \dfrac{F}{E_1} x & x \le \tfrac12 \\[2mm] \dfrac{F}{2 E_1} + \dfrac{F}{E_2}\left(x - \tfrac12\right) & x > \tfrac12 \end{cases}$

Our FEM solution should reproduce this exactly (P1 elements are exact for piecewise-linear solutions).

Step 1: Mesh and Materials

We divide the bar into ne equal elements. Each element gets a Young's modulus depending on which half it lies in.

Code [1]Run

ne = 20                  # number of elements
nn = ne + 1              # number of nodes
h = 1.0 / ne             # element length
F = 100.0                # tension force at the right end

E1 = 200.0               # steel (left half)
E2 = 70.0                # aluminum (right half)

# Per-element material: E[e] for element e spanning [x_e, x_e + h]
E = zeros(ne)
for e in 1:ne {
  xmid = (e - 0.5) * h
  if xmid < 0.5 {
    E[e] = E1
  } else {
    E[e] = E2
  }
}

println("=== Bar ===")
println("Elements: ${ne}, nodes: ${nn}")
println("Material split at x = 0.5: E = ${E1} | ${E2}")

Step 2: Element Stiffness and Assembly

For a 1D P1 element of length $h$ and stiffness $E$ , the element stiffness matrix is

$k_e = \frac{E_e}{h} \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix}$

We loop over elements and add each 2×2 block into the global matrix as triplets. Shared nodes receive contributions from both neighboring elements — sparse sums duplicate entries, which is exactly FEM assembly.

Code [2]Run

# 4 triplets per element
rows = zeros(4 * ne)
cols = zeros(4 * ne)
vals = zeros(4 * ne)

idx = 1
for e in 1:ne {
  ke = E[e] / h        # scalar stiffness factor
  n1 = e               # left node of element e
  n2 = e + 1           # right node

  rows[idx] = n1
  cols[idx] = n1
  vals[idx] = ke
  idx = idx + 1

  rows[idx] = n1
  cols[idx] = n2
  vals[idx] = -ke
  idx = idx + 1

  rows[idx] = n2
  cols[idx] = n1
  vals[idx] = -ke
  idx = idx + 1

  rows[idx] = n2
  cols[idx] = n2
  vals[idx] = ke
  idx = idx + 1
}

K = sparse(rows, cols, vals, nn, nn)
println("Global stiffness assembled: ${nn} x ${nn}, ${nnz(K)} nonzeros")

Step 3: Boundary Conditions

Two kinds of boundary conditions, treated differently:

Dirichlet (essential): $u(0) = 0$ — eliminated from the system by dropping node 1.
Neumann (natural): traction $F$ at $x = 1$ — appears in the load vector at the last node.

Code [3]Run

# Reduced system: unknowns are nodes 2..nn (node 1 is clamped)
nred = nn - 1

rr = zeros(4 * ne)
cc = zeros(4 * ne)
vv = zeros(4 * ne)
cnt = 0
for t in 1:(idx - 1) {
  i = rows[t]
  j = cols[t]
  if i > 1 {
    if j > 1 {
      cnt = cnt + 1
      rr[cnt] = i - 1
      cc[cnt] = j - 1
      vv[cnt] = vals[t]
    }
  }
}

r3 = zeros(cnt)
c3 = zeros(cnt)
v3 = zeros(cnt)
for t in 1:cnt {
  r3[t] = rr[t]
  c3[t] = cc[t]
  v3[t] = vv[t]
}
Kred = sparse(r3, c3, v3, nred, nred)

# Load vector: point force F at the last node
b = zeros(nred)
b[nred] = F

println("Clamped node eliminated: system is ${nred} x ${nred}")
println("Traction F = ${F} applied at x = 1")

Step 4: Solve

The reduced stiffness matrix is symmetric positive definite — solve directly, then rebuild the full displacement vector including the clamped node.

Code [4]Run

ured = spsolve(Kred, b)

# Full displacement field (node 1 clamped at 0)
u = zeros(nn)
for k in 1:nred {
  u[k + 1] = ured[k]
}

println("=== Solution ===")
println("Tip displacement u(1) = ${u[nn]}")

# Exact: F/(2 E1) + F/(2 E2)
u_exact = F / (2.0 * E1) + F / (2.0 * E2)
println("Exact tip displacement = ${u_exact}")
println("Error: ${abs(u[nn] - u_exact)}")

Step 5: Stress Recovery

Stress is computed per element from the strain $\varepsilon_e = (u_{n+1} - u_n)/h$ . Equilibrium demands constant stress $\sigma = F$ across the whole bar — including across the material interface.

Code [5]Run

sigma = zeros(ne)
strain = zeros(ne)
for e in 1:ne {
  strain[e] = (u[e + 1] - u[e]) / h
  sigma[e] = E[e] * strain[e]
}

println("Stress in first element:  ${sigma[1]}   (steel)")
println("Stress in last element:   ${sigma[ne]}   (aluminum)")
println("Strain in first element:  ${strain[1]}")
println("Strain in last element:   ${strain[ne]}")

The stress is the same in both materials, but the soft aluminum strains $E_1/E_2 \approx 2.9\times$ more — that's the strain jump at the interface.

Step 6: Visualize Displacement and Stress

Code [6]Run

x = linspace(0.0, 1.0, nn)

plot(x, u)
title("Displacement u(x) — kink at the material interface")
xlabel("x")
ylabel("u")

Code [7]Run

# Element midpoints for the piecewise-constant stress
xm = zeros(ne)
for e in 1:ne {
  xm[e] = (e - 0.5) * h
}

plot(xm, sigma)
title("Stress σ(x) — constant, as equilibrium demands")
xlabel("x")
ylabel("σ")

Step 7: Convergence Sanity Check

P1 elements represent piecewise-linear fields exactly, so the FEM solution is exact for any number of elements — as long as the material interface falls on an element boundary. Verify with a coarse mesh:

Code [8]Run

# 2 elements: one per material
ke1 = E1 / 0.5
ke2 = E2 / 0.5

Kc = sparse([1.0, 1.0, 2.0, 2.0], [1.0, 2.0, 1.0, 2.0], [ke1 + ke2, -ke2, -ke2, ke2], 2, 2)
bc = [0.0, F]
uc = spsolve(Kc, bc)

println("Tip displacement with 2 elements:  ${uc[2]}")
println("Tip displacement with ${ne} elements: ${u[nn]}")
println("Both match the exact value ${u_exact}")

Summary

The elasticity workflow mirrors every FEM problem:

Constitutive law — Hooke's law links stress and strain via material stiffness $E$
Element matrices — $\frac{E_e}{h}\begin{pmatrix}1 & -1\\ -1 & 1\end{pmatrix}$ per element, materials vary element-by-element
Assembly — triplets with duplicate summing; shared nodes accumulate both neighbors
Boundary conditions — Dirichlet by elimination, Neumann via the load vector
Solve & recover — displacement from spsolve, stress per element from the strain
Verify — constant stress and the exact tip displacement confirm the physics

In 2D/3D the element matrices grow (vector displacements, the full strain tensor) but the pipeline — materials, assembly, BCs, solve, stress recovery — is identical.

Workbench

Clear

No variables in workbench

Next Steps

Finite Elements 1: Poisson Problem

FEM fundamentals: weak forms, assembly, sparse solves

Heat Conduction

Time-dependent problems with implicit stepping