RasmahRasmah

Linear Elasticity with the Finite Element Method

using Rasmah
using LinearAlgebra

Overview

Every load-bearing part — a bracket, a strap, a weldment — starts with the same design question: given the loads, will it hold, and will it hold without bending too much? This is the canonical Rasmah simulation workflow. We build a part, select its boundary faces, mesh it, apply the service loads, and solve for the displacement field — then run the two checks an engineer actually cares about: stress against the material’s yield (a factor of safety) and displacement against an allowable deflection.

The part is a tension strap — a 100 × 30 × 5 mm steel plate, the kind that takes a lifting load in a bracket assembly — pulled by 10 kN at its ends. We work in millimetres and megapascals, the units a mechanical engineer thinks in.

The workflow follows a fixed order of operations, which every finite-element and topology-optimization tutorial in Rasmah repeats:

shape → boundary sets → mesh → boundary conditions → solve

A quick reminder of the physics

Linear elasticity is the simplest solid-mechanics model. Under a load, a material point at x moves to x + u(x); the symmetric strain is

\[ \varepsilon(u) = \frac{1}{2}\big(\nabla u + \nabla u^{\mathsf T}\big), \]

and for an isotropic material the stress is linear in the strain (Hooke’s law),

\[ \sigma = 2\mu\,\varepsilon + \lambda\,\operatorname{tr}(\varepsilon)\,I, \]

where μ and λ are the Lamé parameters (equivalently the Young’s modulus E and Poisson ratio ν, which steel supplies). The strong form is ∇·σ + f = 0 with u = g on the fixed boundary; the finite-element method turns this into the sparse linear system K u = f, which solve assembles and solves for you.

Geometry, face selection, and meshing

# 1. the geometry — a 100 × 30 × 5 mm steel strap, drawn in millimetres
part = box(100u"mm", 30u"mm", 5u"mm")

# 2. exact boundary representation + face selection (before meshing)
b = brep(part)
left  = faces_in_box(b, [-55u"mm", -20u"mm", -5u"mm"], [-45u"mm", 20u"mm", 5u"mm"])
right = faces_in_box(b, [ 45u"mm", -20u"mm", -5u"mm"], [ 55u"mm", 20u"mm", 5u"mm"])

# 3. mesh the geometry — the face sets carry over
mesh = tetrahedralize(b, 10)

println("nodes:    ", size(mesh.nodes, 2))
println("elements: ", size(mesh.elements, 2))
nodes:    88
elements: 180

Visualize the geometry:

vtk(b)

Solving for displacement

Clamp the left end and pull the right end with the service load — a total of 10 kN, applied as a uniform traction of 10 kN / (30 × 5 mm²) = 66.7 MPa. Elasticity() is a vector field :u; steel supplies E and ν, and the region => value pair syntax resolves the BRep face sets onto the mesh:

u = solve(Elasticity(), mesh, steel;
          fixed = left => [0.0, 0.0, 0.0],
          loads = right => [66.6667u"MPa", 0.0u"MPa", 0.0u"MPa"])

println("maximum |displacement| = ", round(uconvert(u"mm", maximum(abs, u)); digits = 4), " mm")
maximum |displacement| = 0.0334 mm

Colour the mesh by the per-node displacement magnitude:

umag = [norm(u[(3i-2):(3i)]) for i in 1:size(mesh.nodes, 2)]
vtk(mesh; field = umag, fieldname = "|displacement|")

What the API did

  • faces_in_box(b, …) selected the two faces on the BRep, before meshing; tetrahedralize(b, 10) carried those face sets onto the mesh.
  • fixed = left => [0, 0, 0] prescribed the Dirichlet support (one vector repeated per node); loads = right => [66.7, 0, 0] MPa applied a surface traction, entered in megapascals and converted to SI inside the solver.
  • steel supplied E and ν; the default element order (1) and direct solver were used. Pass order, method = :amg / :gmg / :gmg_mf, or levels to change them.

The same problem solved with an algebraic-multigrid solver agrees:

u_amg = solve(Elasticity(), mesh, steel;
              fixed = left => [0.0, 0.0, 0.0], loads = right => [66.6667u"MPa", 0.0u"MPa", 0.0u"MPa"],
              method = :amg)
println("amg agrees: ", isapprox(ustrip(u), ustrip(u_amg); atol = 1e-6))
amg agrees: true

Design checks

The point of the simulation is a pass/fail decision. Two checks matter for a strap like this.

Displacement against the hand calculation

A bar under pure tension stretches by δ = FL/(EA). With F = 10 kN, L = 100 mm, E = 200 GPa and A = 150 mm²:

δ_hand = 10000.0 * 0.1 / (200.0e9 * 150.0e-6) * 1e3   # mm
println("hand calculation:  ", round(δ_hand; digits = 4), " mm")
println("finite elements:  ", round(uconvert(u"mm", maximum(abs, u)); digits = 4), " mm")
hand calculation:  0.0333 mm
finite elements:  0.0334 mm

The answers agree to well under a percent: the finite-element result is a fraction larger because clamping the whole x = 0 face also restrains the lateral Poisson contraction near the support.

Stress against the yield — the factor of safety

von_mises_stress evaluates the stress element-by-element from the displacement. The von Mises (equivalent) stress is what a ductile material yields under:

σv = von_mises_stress(mesh, ustrip(u), stiffness(steel))

# the stress at mid-span, where the field is uniform
mid   = [sum(mesh.nodes[1, mesh.elements[:, e]]) / 4 for e in axes(mesh.elements, 2)]
σmid  = maximum(σv[i] for i in eachindex(mid) if 0.04 < mid[i] < 0.06)

println("mid-span von Mises: ", round(σmid / 1e6; digits = 1), " MPa")
println("hand calculation:   ", round(10000.0 / 150.0; digits = 1), " MPa")
println("factor of safety:   ", round(yield_strength(steel) / σmid; digits = 2), "  (yield: ", round(yield_strength(steel) / 1e6; digits = 0), " MPa)")
mid-span von Mises: 66.7 MPa
hand calculation:   66.7 MPa
factor of safety:   3.75  (yield: 250.0 MPa)

The mid-span stress equals F/A = 66.7 MPa — the mesh recovers the analytic value exactly. A36 structural steel yields at 250 MPa, so the factor of safety is about 3.7: the strap passes with margin.

Notice the deliberate choice of where we read the stress. The maximum over the whole part sits at the corners of the clamped face:

println("maximum von Mises: ", round(maximum(σv) / 1e6; digits = 1), " MPa (at the clamped corners)")
maximum von Mises: 70.2 MPa (at the clamped corners)
vtk(mesh; field = σv, fieldname = "von Mises stress", colormap = :inferno, scalarbar = true)

That corner peak is a stress concentration, not a load-bearing concern: fixing every degree of freedom on a face restrains the corner in a way a real bolt or weld never could, and the peak grows as the mesh is refined instead of converging. The engineering rule is to evaluate stress in the uniform-stress field away from supports and load-application points.

The deflection limit

The strap is a positioning member: the application allows 0.05 mm of stretch.

δ_max = uconvert(u"mm", maximum(abs, u))
println("stretch = ", round(δ_max; digits = 4), " mm;  allowed = 0.05 mm;  passes: ", δ_max < 0.05)
stretch = 0.0334 mm;  allowed = 0.05 mm;  passes: true

Did I mesh finely enough?

A mesh-convergence check is part of any credible FEA: rerun the solve at finer resolutions and confirm the quantity you report stops changing. For the strap, both the stretch and the mid-span stress are already converged at the coarsest mesh:

for res in (10, 20, 40)
    m = tetrahedralize(b, res)
    u_ = solve(Elasticity(), m, steel;
               fixed = left => [0.0, 0.0, 0.0], loads = right => [66.6667u"MPa", 0.0u"MPa", 0.0u"MPa"])
    σv_ = von_mises_stress(m, ustrip(u_), stiffness(steel))
    mid_ = [sum(m.nodes[1, m.elements[:, e]]) / 4 for e in axes(m.elements, 2)]
    σmid_ = maximum(σv_[i] for i in eachindex(mid_) if 0.04 < mid_[i] < 0.06)
    println("res = ", res, ":  δ = ", round(uconvert(u"mm", maximum(abs, u_)); digits = 4),
            " mm,  mid-span σ = ", round(σmid_ / 1e6; digits = 2), " MPa")
end
res = 10:  δ = 0.0334 mm,  mid-span σ = 66.67 MPa
res = 20:  δ = 0.0333 mm,  mid-span σ = 66.67 MPa
res = 40:  δ = 0.0333 mm,  mid-span σ = 66.67 MPa

The stretch settles at 0.0333 mm and the stress at 66.67 MPa — the answers at resolution 10 were already accurate, which is typical for a part in pure tension. (A part with stress gradients — a hole, a fillet, a crack — needs local refinement instead; see adaptive refinement for that case.)

Compliance

Compliance c = fᵀu measures how much work the load does — a stiffer structure has lower compliance. It is the standard “make it stiff” objective:

loads = face_load(mesh, faces_in_box(mesh, [45e-3, -20e-3, -5e-3], [55e-3, 20e-3, 5e-3]), [66.6667e6, 0.0, 0.0])
c = dot(loads, ustrip(u))
println("compliance = ", round(c; digits = 4), " J")
println("applied load total = ", round(sum(loads[1:3:end]) / 1e3; digits = 2), " kN")
compliance = 0.3289 J
applied load total = 10.0 kN

The last line is the engineer’s bookkeeping check: the forces spread over the pulled face must sum to the 10 kN we applied — they do, exactly. (The equal and opposite 10 kN is what the fixed face reacts with; every supported simulation should balance forces this way.)

Next steps