Solving with physics

Overview

The assemblers in The finite element method give you every building block, but a real simulation is always the same sequence: pick a finite-element space, write the weak form, assemble, apply boundary conditions, solve. Rasmah packages that whole sequence behind a single call:

solve(Physics(), mesh, material; boundary_conditions...)

You name the physics, hand over a mesh and a material, and describe the boundary conditions. Rasmah picks the element order, the weak form, and the solver, and returns the solved field (or fields). This chapter covers the mechanics of that one line.

Naming the physics

A physics is just a type that names a problem. Passing it to solve selects the governing equation, the boundary conditions, and the defaults.

AbstractPhysics

julia> Elasticity() isa AbstractPhysics
true

Each marker has its own solve(physics, mesh, material; …) method. The signature is the same everywhere — only the returned fields and the required material properties differ.

Linear elasticity

Elasticity solves for the displacement field. The material must supply $E$ and $\nu$:

Elasticity

julia> m = tetrahedralize_box(1.0, 1.0, 1.0, 4, 4, 4);

julia> left  = nodes_in_box(m, [0.0, 0.0, 0.0], [0.0, 1.0, 1.0]);

julia> right = nodes_in_box(m, [1.0, 0.0, 0.0], [1.0, 1.0, 1.0]);

julia> steel = Material(youngs_modulus = 200e9, poisson_ratio = 0.3);

julia> u = solve(Elasticity(), m, steel; fixed = left, loads = right => (1e8, 0.0, 0.0));

julia> u isa QuantityVector && length(u) == 3 * size(m.nodes, 2)
true

Boundary conditions use the region => value syntax: a PointSet/FaceSet region paired with a prescribed value. fixed is a Dirichlet support (a fixed displacement), loads a Neumann load (a force or traction).

Heat conduction

HeatConduction solves for the temperature field from the material's conductivity $k$:

HeatConduction

julia> m = tetrahedralize_box(1.0, 1.0, 1.0, 4, 4, 4);

julia> left  = nodes_in_box(m, [0.0, 0.0, 0.0], [0.0, 1.0, 1.0]);

julia> right = nodes_in_box(m, [1.0, 0.0, 0.0], [1.0, 1.0, 1.0]);

julia> copper = Material(youngs_modulus = 120e9, poisson_ratio = 0.34,
                         thermal_conductivity = 400.0);

julia> T = solve(HeatConduction(), m, copper; fixed = left => 0.0, loads = right => 1.0);

julia> T isa QuantityVector && length(T) == size(m.nodes, 2)
true

The same region => value syntax applies — fixed = left => 0.0 pins the temperature on the left face to zero, loads = right => q imposes a heat flux.

Coupled physics

Multi-field problems return a NamedTuple of fields. ThermoElasticity couples the temperature and displacement fields one-way (thermal strain drives deformation):

ThermoElasticity

julia> m = tetrahedralize_box(1.0, 1.0, 1.0, 4, 4, 4);

julia> left = nodes_in_box(m, [0.0, 0.0, 0.0], [0.0, 1.0, 1.0]);

julia> steel = Material(youngs_modulus = 200e9, poisson_ratio = 0.3,
                        thermal_conductivity = 50.0, thermal_expansion = 1.2e-5);

julia> r = solve(ThermoElasticity(), m, steel; fixed_u = left, fixed_T = left => 0.0);

julia> keys(r) == (:u, :T) && r.u isa QuantityVector && r.T isa QuantityVector
true

m = tetrahedralize_box(1.0, 1.0, 1.0, 4, 4, 4)
left  = nodes_in_box(m, [0.0, 0.0, 0.0], [0.0, 1.0, 1.0])
right = nodes_in_box(m, [1.0, 0.0, 0.0], [1.0, 1.0, 1.0])

steel = Material(youngs_modulus = 200e9, poisson_ratio = 0.3,
                 thermal_conductivity = 50.0, thermal_expansion = 1.2e-5)

r = solve(ThermoElasticity(), m, steel;
          fixed_u = left,
          fixed_T = (left => 0.0, right => 100.0))

(u = [0.0 m, 0.0 m, 0.0 m, -5.742605809004645e-5 m, -0.00015733452373501415 m, -0.00015733452373501504 m, 9.989056318185969e-6 m, -0.00031213410666607375 m, -0.00031213410666607614 m, 0.0001824390299229943 m, …, 0.00012621475477644732 m, 7.095747117668785e-5 m, 0.0002670143439750902 m, 0.0002670143439750886 m, 0.0002696389058521776 m, 0.0004011323133777439 m, 0.00040113231337774076 m, 0.0005511286158403335 m, 0.0005034453675248777 m, 0.0005034453675248737 m], T = [0.0 K, 25.00000000000005 K, 50.00000000000008 K, 75.00000000000006 K, 100.0 K, 0.0 K, 25.000000000000064 K, 50.000000000000114 K, 75.00000000000009 K, 100.0 K, …, 0.0 K, 25.00000000000004 K, 50.000000000000064 K, 75.00000000000004 K, 100.0 K, 0.0 K, 25.000000000000032 K, 50.00000000000006 K, 75.00000000000007 K, 100.0 K])

round(maximum(abs.(r.T)); digits = 1)

100.0 K

The $100$ K temperature rise produces a thermal strain, and therefore a displacement field, even though no mechanical load was applied:

round(maximum(abs.(r.u)); digits = 10)

0.0007144484 m

The hot end ($+100$ K) expands against the cold, clamped end, so the free end moves in the $+x$ direction. The ideal free-expansion estimate $\alpha\,\Delta T\,L = 1.2 \times 10^{-5} \cdot 100 \cdot 1 = 1.2 \times 10^{-3}$ m is an upper bound: clamping the whole $x = 0$ face also restrains the lateral Poisson contraction, so the axial stretch comes out a little smaller ($\approx 7.1 \times 10^{-4}$ m).

The physics catalogue

Elasticity, heat, and thermo-elasticity are three of a larger family. The same marker pattern — solve(Physics(), mesh, material; …) — dispatches across thirteen built-in physics, each a small struct naming a governing equation, its solved fields, and the material properties it consumes:

A single-field physics returns its field; a coupled physics returns a NamedTuple ((; u, p), (; u, T), …). The flow solvers take a FluidMaterial and read its viscosity and density; the solid solvers take a Material. Every marker is documented in the reference and uses the same region => value boundary-condition syntax shown above, so switching physics is just changing the marker (and supplying the right material).

Materials

The third argument to solve is any AbstractMaterial. Rasmah ships a library of ready-made materials (steel, aluminum, copper, …) and a Material(...) constructor for custom properties. The material decides which coefficients the physics uses — Elasticity reads $E, \nu$, HeatConduction reads $k$, and ThermoElasticity reads all four. The Materials chapter covers them in detail.

A worked example: temperature through a wall

The one-line solve makes it easy to check a result against theory. Steady heat conduction through a unit slab with the two faces held at $0$ and $1$ has the exact linear solution $T(x) = x$. Solve it and confirm the mid-plane temperature is $0.5$:

m = tetrahedralize_box(1.0, 1.0, 1.0, 4, 4, 4)
left  = nodes_in_box(m, [0.0, 0.0, 0.0], [0.0, 1.0, 1.0])
right = nodes_in_box(m, [1.0, 0.0, 0.0], [1.0, 1.0, 1.0])

copper = Material(youngs_modulus = 120e9, poisson_ratio = 0.34,
                  thermal_conductivity = 400.0)

T = solve(HeatConduction(), m, copper; fixed = (left => 0.0, right => 1.0))

mid = nodes_in_box(m, [0.5, 0.0, 0.0], [0.5, 1.0, 1.0])

PointSet([3, 8, 13, 18, 23, 28, 33, 38, 43, 48  …  78, 83, 88, 93, 98, 103, 108, 113, 118, 123])

all(isapprox(T[i], 0.5u"K"; atol = 1e-12) for i in mid.ids)

true

Colour the slab by temperature — a smooth linear ramp from the cold to the hot face:

vtk(m; field = T, fieldname = "temperature", edges = false)

Next steps

With solve in hand you can run any built-in physics. The natural next chapters are the Topology optimization that uses these solves as its inner loop, and Materials, which documents the material library that feeds them.