Thermal Stress in a Steel Column (Thermoelasticity u + T)

Overview

Every engineer has seen the consequences of thermal expansion: railway rails that buckle in a heat wave, concrete slabs that crack because nobody left an expansion joint, engine parts that seize when they get hot. The physics is simple — a material expands by ε = α ΔT — but the stresses it produces depend entirely on how the part is held: a free part grows quietly, while a constrained one develops stresses that can reach the yield strength.

This tutorial runs the two canonical cases on a real part — a 1 m length of 50 × 50 mm steel square bar, the section of a window mullion or a welded frame member. ThermoElasticity() couples the displacement field :u and the temperature field :T; the material supplies E, ν, k and the thermal-expansion coefficient α, and a temperature field induces stresses through the coupling β = (3λ + 2μ) α.

The workflow is the canonical one:

shape → boundary sets → mesh → boundary conditions → solve

The physics

The thermal strain is ε_th = α ΔT in every direction. If the part is free to expand, it simply grows: a bar of length L heated uniformly by ΔT stretches to δ = α ΔT L — and carries no stress at all. If it is fully constrained, it cannot grow, and the supports must push back with the force that would be needed to compress the free expansion back to zero — a uniform stress

\[ \sigma_{th} = E\, \alpha\, \Delta T \]

that depends only on the material and the temperature swing, never on the cross-section. The thermo-elastic weak form adds the temperature gradient as a body-like load −β T ∇·v on top of the elastic K u = f, which is why the system is block-triangular (one-way): the temperature solves first, then feeds the mechanical problem.

Drawing the member

The member is a 1 m length of 50 × 50 mm steel square bar — the kind of section used for window mullions, railings, and welded frames — drawn as a single box, in metres:

L, s = 1.0, 0.050                       # 1 m long, 50 × 50 mm section (m)

bar = box(L, s, s)
b = brep(bar)

BRep{Vertex{Float64}, Edge{LineCurve{Float64}}, Loop, Face{PlaneSurface{Float64}}, Shell, Solid}(Vertex{Float64}[Vertex{Float64}([-0.5, -0.025, -0.025]), Vertex{Float64}([0.5, -0.025, -0.025]), Vertex{Float64}([0.5, 0.025, -0.025]), Vertex{Float64}([-0.5, 0.025, -0.025]), Vertex{Float64}([-0.5, -0.025, 0.025]), Vertex{Float64}([0.5, -0.025, 0.025]), Vertex{Float64}([0.5, 0.025, 0.025]), Vertex{Float64}([-0.5, 0.025, 0.025])], Edge{LineCurve{Float64}}[Edge{LineCurve{Float64}}(1, 2, LineCurve{Float64}([-0.5, -0.025, -0.025], [0.5, -0.025, -0.025])), Edge{LineCurve{Float64}}(2, 3, LineCurve{Float64}([0.5, -0.025, -0.025], [0.5, 0.025, -0.025])), Edge{LineCurve{Float64}}(3, 4, LineCurve{Float64}([0.5, 0.025, -0.025], [-0.5, 0.025, -0.025])), Edge{LineCurve{Float64}}(4, 1, LineCurve{Float64}([-0.5, 0.025, -0.025], [-0.5, -0.025, -0.025])), Edge{LineCurve{Float64}}(5, 6, LineCurve{Float64}([-0.5, -0.025, 0.025], [0.5, -0.025, 0.025])), Edge{LineCurve{Float64}}(6, 7, LineCurve{Float64}([0.5, -0.025, 0.025], [0.5, 0.025, 0.025])), Edge{LineCurve{Float64}}(7, 8, LineCurve{Float64}([0.5, 0.025, 0.025], [-0.5, 0.025, 0.025])), Edge{LineCurve{Float64}}(8, 5, LineCurve{Float64}([-0.5, 0.025, 0.025], [-0.5, -0.025, 0.025])), Edge{LineCurve{Float64}}(1, 5, LineCurve{Float64}([-0.5, -0.025, -0.025], [-0.5, -0.025, 0.025])), Edge{LineCurve{Float64}}(2, 6, LineCurve{Float64}([0.5, -0.025, -0.025], [0.5, -0.025, 0.025])), Edge{LineCurve{Float64}}(3, 7, LineCurve{Float64}([0.5, 0.025, -0.025], [0.5, 0.025, 0.025])), Edge{LineCurve{Float64}}(4, 8, LineCurve{Float64}([-0.5, 0.025, -0.025], [-0.5, 0.025, 0.025]))], Loop[Loop([-1, -4, -3, -2]), Loop([5, 6, 7, 8]), Loop([1, 10, -5, -9]), Loop([3, 12, -7, -11]), Loop([9, -8, -12, 4]), Loop([2, 11, -6, -10])], Face{PlaneSurface{Float64}}[Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.5, -0.025, -0.025], [1.0, 0.0, 0.0], [0.0, 0.05, 0.0]), [1]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.5, -0.025, 0.025], [1.0, 0.0, 0.0], [0.0, 0.05, 0.0]), [2]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.5, -0.025, -0.025], [1.0, 0.0, 0.0], [0.0, 0.0, 0.05]), [3]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.5, 0.025, -0.025], [1.0, 0.0, 0.0], [0.0, 0.0, 0.05]), [4]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([-0.5, -0.025, -0.025], [0.0, 0.05, 0.0], [0.0, 0.0, 0.05]), [5]), Face{PlaneSurface{Float64}}(PlaneSurface{Float64}([0.5, -0.025, -0.025], [0.0, 0.05, 0.0], [0.0, 0.0, 0.05]), [6])], Shell[Shell([1, 2, 3, 4, 5, 6])], Solid[Solid([1])])

Because the part is an axis-aligned box, the canonical workflow applies end to end: the exact BRep carries its face sets onto the mesh, so the ends are selected before meshing.

vtk(b)

Zoom Opacity

Boundary sets and meshing

The bar is welded into a frame at its left end; the right end is free (or later, clamped for the constrained case). Both are the square cross-sections at x = ±L/2:

mesh = tetrahedralize(b, 10)

left  = faces_in_box(b, [-0.501, -0.026, -0.026], [-0.499, 0.026, 0.026])
right = faces_in_box(b, [ 0.499, -0.026, -0.026], [ 0.501, 0.026, 0.026])

println("nodes: ", size(mesh.nodes, 2), ", elements: ", size(mesh.elements, 2))
println("left face: ", length(left.ids), " BRep faces,  right face: ", length(right.ids), " BRep faces")

nodes: 44, elements: 60
left face: 1 BRep faces,  right face: 1 BRep faces

Case A — free thermal expansion

Clamp the left end, heat the whole bar uniformly by ΔT = 100 K, and measure the stretch. The right end should move by exactly α ΔT L:

ΔT = 100.0
α  = thermal_expansion(steel)

u_free = solve(ThermoElasticity(), mesh, steel;
               fixed_u = left => [0.0, 0.0, 0.0],
               fixed_T = [left => ΔT, right => ΔT])

# the gauge check: the expansion per metre, measured mid-bar away from the
# welded end (the clamp restrains the lateral Poisson contraction there)
i1 = argmin([abs(mesh.nodes[1, i] - 0.1) + abs(mesh.nodes[2, i]) + abs(mesh.nodes[3, i]) for i in axes(mesh.nodes, 2)])
i2 = argmin([abs(mesh.nodes[1, i] - 0.4) + abs(mesh.nodes[2, i]) + abs(mesh.nodes[3, i]) for i in axes(mesh.nodes, 2)])
uu = ustrip(u_free.u)

println("FEM strain:   ", round((uu[3i2 - 2] - uu[3i1 - 2]) / 0.3; digits = 8))
println("hand, α ΔT:   ", round(α * ΔT; digits = 8))

FEM strain:   0.00117
hand, α ΔT:   0.00117

The mid-bar strain matches α ΔT to numerical precision — and because nothing resists the expansion, the stress is zero everywhere. This is the quiet case.

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

Zoom Opacity

Case B — the constrained member

Now weld both ends into the frame and heat it again. The solver returns only the small lateral bulge of the unclamped sides — the end-to-end stretch is gone — and that does not mean the member is unstressed: it means the supports are holding it in compression. The stress they must provide is σ = E α ΔT, exactly the stress that would cancel the free expansion:

u_c = solve(ThermoElasticity(), mesh, steel;
            fixed_u = [left => [0.0, 0.0, 0.0], right => [0.0, 0.0, 0.0]],
            fixed_T = [left => ΔT, right => ΔT])

println("max |u| with both ends welded: ", round(maximum(abs, ustrip.(u_c.u)); digits = 6), " m  (just the lateral bulge)")

E = youngs_modulus(steel)
σth = E * α * ΔT
println("thermal stress σ = E α ΔT:        ", round(σth / 1e6; digits = 1), " MPa")
println("yield strength of the steel:       ", round(yield_strength(steel) / 1e6; digits = 0), " MPa")

max |u| with both ends welded: 9.7e-5 m  (just the lateral bulge)
thermal stress σ = E α ΔT:        234.0 MPa
yield strength of the steel:       250.0 MPa

At a 100 K swing the constrained stress is ≈ 240 MPa — within a few percent of the steel’s 250 MPa yield. This is the rail story in one number: a welded track that cannot expand on a 100 K day carries nearly its yield strength in compression, which is why rails get expansion joints and why they buckle when those joints jam.

Case C — the service case (pull + temperature gradient)

Real parts see a gradient, not a uniform swing: the bar end at the furnace side runs 100 K hotter than the far end, while the frame pulls it with a 1 MPa traction. The one-way coupling solves the temperature first and then the mechanics:

res = solve(ThermoElasticity(), mesh, steel;
            fixed_u = left => [0.0, 0.0, 0.0],
            loads_u = right => [1e6, 0.0, 0.0],
            fixed_T = [left => ΔT, right => 0.0])
println("max |u| = ", round(maximum(abs, ustrip.(res.u)); digits = 5),
        ", T range = ", round(minimum(ustrip.(res.T)); digits = 2), " … ",
        round(maximum(ustrip.(res.T)); digits = 2), " K")

max |u| = 0.00066, T range = 0.0 … 100.0 K

vtk(mesh; field = res.T, fieldname = "temperature", colormap = :inferno)

Zoom Opacity

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

Zoom Opacity

Design (topology optimization)

ThermoElasticTopOpt() minimizes the thermo-mechanical compliance with the density method — the density scales the coupled block (elasticity + thermal coupling), solved with LU:

φ, hist = solve(ThermoElasticTopOpt(), mesh, steel;
                fixed_u = left => [0.0, 0.0, 0.0], loads_u = right => [1e6, 0.0, 0.0],
                fixed_T = [left => ΔT, right => 0.0],
                volfrac = 0.5, max_iter = 40)

println("final thermo-mechanical compliance = ", round(hist[end]; digits = 4))
println("final volume fraction              = ", round(sum(φ) / length(φ); digits = 4))

final thermo-mechanical compliance = 0.0313
final volume fraction              = 0.4991

vtk(mesh; field = φ, fieldname = "thermoelastic density")

Zoom Opacity

Notes

The solution is returned as a NamedTuple (; u, T) with the fields in the same order as field_names(ThermoElasticity()) == (:u, :T).
The coupled system is block-triangular (one-way thermal→mechanical), so the default solver is LU.
fixed_T = face => value prescribes the temperature on a face set; a uniform temperature on both ends gives a uniform field, a hot/cold pair gives a linear gradient.

Next steps

Linear elasticity — the K u = f core that the thermal coupling feeds.
Thermoelectricity — the same Laplace temperature field, now driving the Seebeck voltage.
Heat conduction — the scalar heat problem behind the temperature half of this tutorial.