Fluid flow

Overview

The fields solved so far — displacement, temperature, potential — were scalars or a single vector on a solid. A fluid is different in one crucial way: it is incompressible. The velocity field $\mathbf{u}$ and the pressure $p$ are coupled, because the pressure is exactly the force that keeps the fluid from compressing. This chapter covers the two built-in flow solvers — creeping Stokes flow and inertial Navier–Stokes flow.

The shared idea is a force balance with a constraint. The momentum equation says viscous forces balance the pressure gradient (plus inertia, for Navier–Stokes); the continuity equation says the velocity is divergence-free. The two together form a saddle-point problem, and Rasmah solves it with Taylor–Hood elements (quadratic velocity, linear pressure).

Creeping flow: Stokes

When a flow is slow — vanishing Reynolds number — inertia is negligible, and the momentum balance is purely viscous versus pressure:

\[\mu \nabla^2 \mathbf{u} = \nabla p, \qquad \nabla \cdot \mathbf{u} = 0 .\]

Rasmah.StokesFlow — Type

StokesFlow

Stokes flow: solves for the velocity :u and pressure :p of a creeping Newtonian fluid, $\mu \nabla^2 \mathbf{u} = \nabla p$ with $\nabla \cdot \mathbf{u} = 0$.

Theory

When a flow is slow (vanishing Reynolds number), inertia is negligible and the momentum balance reduces to the Stokes equations: the viscous term balances the pressure gradient, and the velocity field is divergence-free (incompressible). The two fields are coupled — pressure is the Lagrange multiplier that enforces $\nabla \cdot \mathbf{u} = 0$ — so the discrete problem is a saddle-point system solved with Taylor–Hood (P2 velocity / P1 pressure) elements.

Solving

solve(StokesFlow(), m, fluid; fixed_u, loads_u, fixed_p, …) returns (; u, p). The FluidMaterial supplies the dynamic viscosity $\mu$; fixed_u pins the velocity (no-slip walls), fixed_p pins the pressure (inlet/outlet).

source

Inertial flow: Navier–Stokes

Add the convective term $\rho(\mathbf{u}\cdot\nabla)\mathbf{u}$ and the problem becomes the full Navier–Stokes equations — nonlinear, so the solve iterates:

\[\rho (\mathbf{u} \cdot \nabla) \mathbf{u} = \mu \nabla^2 \mathbf{u} - \nabla p, \qquad \nabla \cdot \mathbf{u} = 0 .\]

Rasmah.NavierStokes — Type

NavierStokes

Steady incompressible Navier–Stokes: solves for the velocity :u and pressure :p, $\rho (\mathbf{u} \cdot \nabla) \mathbf{u} = \mu \nabla^2 \mathbf{u} - \nabla p$ with $\nabla \cdot \mathbf{u} = 0$.

Theory

Stokes flow ignores inertia; the Navier–Stokes equations keep it, adding the convective term $\rho (\mathbf{u} \cdot \nabla) \mathbf{u}$ to the momentum balance. That term is what makes the problem nonlinear, so the solve iterates (Newton by default) rather than solving one linear system. The density $\rho$ joins the viscosity $\mu$ as the second fluid property.

Solving

solve(NavierStokes(), m, fluid; fixed_u, fixed_p, …) returns (; u, p). The FluidMaterial supplies $\mu$ and $\rho$. As with StokesFlow, fixed_u pins the velocity and fixed_p the pressure.

source

A worked example: pressure-driven channel flow

The classic creeping-flow benchmark is Poiseuille flow through a channel: a pressure drop between inlet and outlet drives a parabolic velocity profile, with the fluid at rest on the walls (no-slip). Set it up on a $[0,1] \times [0,1]$ channel — no-slip velocity on the top and bottom walls, pressure $1$ at the inlet and $0$ at the outlet:

function channel(nx, ny)
    nv = (nx + 1) * (ny + 1)
    verts = Matrix{Float64}(undef, 3, nv)
    for j in 0:ny, i in 0:nx
        verts[:, i + j * (nx + 1) + 1] = [i / nx, j / ny, 0.0]
    end
    faces = Matrix{Int}(undef, 3, 2nx * ny)
    k = 0
    for j in 0:(ny - 1), i in 0:(nx - 1)
        a = i + j * (nx + 1) + 1; b = a + 1; c = a + nx + 1; d = c + 1
        k += 1; faces[:, k] = [a, b, c]
        k += 1; faces[:, k] = [b, d, c]
    end
    TriangleMesh(verts, faces)
end

m = channel(8, 8)
U = FESpace(m, ReferenceFE(2, 2); valuetype = VectorValue{2})
P = FESpace(m, ReferenceFE(2, 1))

ufixed = vcat(dofs(U, BoundingBox([0.0, 0.0], [1.0, 0.0])),
              dofs(U, BoundingBox([0.0, 1.0], [1.0, 1.0])))

inlet  = dofs(P, BoundingBox([0.0, 0.0], [0.0, 1.0]))
outlet = dofs(P, BoundingBox([1.0, 0.0], [1.0, 1.0]))
pfixed = vcat(inlet, outlet)
pg     = vcat(fill(1.0, length(inlet)), fill(0.0, length(outlet)))

res = solve(StokesFlow(), m, water;
            bc = BoundaryConditions(FixedVelocity(ufixed),
                                    FixedPressure(pfixed; value = pg)))

(u = [0.0 m/s, 0.0 m/s, 0.0 m/s, 0.0 m/s, 0.0 m/s, 0.0 m/s, 0.0 m/s, 0.0 m/s, 0.0 m/s, 0.0 m/s, …, -0.08056833688620116 m/s, 0.007993804906480734 m/s, 0.0 m/s, 0.0 m/s, -0.08295036117237477 m/s, -0.008563713245927633 m/s, -0.08204456944473051 m/s, -0.05171650097263466 m/s, 0.0 m/s, 0.0 m/s], p = [1.0 Pa, -0.4875251618478223 Pa, 0.2167956782270357 Pa, -0.06627047555794277 Pa, 0.019580387953753838 Pa, -0.005091302339313065 Pa, 0.0007888098322469244 Pa, -0.0014799638441740626 Pa, 0.0 Pa, 1.0 Pa, …, 0.0 Pa, 1.0 Pa, -0.11092804519383512 Pa, -0.016924589585886714 Pa, 0.0022801498830181226 Pa, -0.0038567052739747814 Pa, -0.001085575942762524 Pa, -0.0004730203824722099 Pa, -0.0013763454524562055 Pa, 0.0 Pa])

round(maximum(abs, res.u); digits = 1)

19.9 m/s

round(minimum(res.p); digits = 3), round(maximum(res.p); digits = 3)

(-0.488 Pa, 1.0 Pa)

The maximum velocity sits mid-channel where the parabola peaks, and the pressure falls from $1$ at the inlet to $0$ at the outlet. Colour the channel by pressure:

vtk(m; field = res.p, fieldname = "pressure", edges = false)

Zoom Opacity

Next steps

Flows couple to heat (natural convection), to fields (MHD), and to structure (fluid–structure interaction) — see the transport and coupled chapters. The solve pattern and Boundary conditions chapter cover the mechanics used here.