RasmahRasmah

Magnetohydrodynamics: Hartmann Flow

Overview

When an electrically conducting fluid flows through a magnetic field, the field induces currents that produce a Lorentz force opposing the flow. In the inductionless (low magnetic Reynolds number) limit this force on a fluid in a uniform field B₀ is

\[F = \sigma\,(u \times B_0) \times B_0\]

For a transverse field B₀ = B₀ ẑ this is a damping −σB₀²u that flattens the channel profile. The exact solution of this Hartmann flow (Hartmann 1937; Davidson 2001) between plates at y = ±h under a pressure gradient G is

\[u(y) = \frac{G}{\sigma B_0^2}\Bigl[1 - \frac{\cosh(Ha\, y/h)}{\cosh(Ha)}\Bigr], \qquad Ha = B_0 h \sqrt{\frac{\sigma}{\mu}}\]

where the Hartmann number Ha measures the magnetic damping relative to viscosity.

using Rasmah
using LinearAlgebra

The channel and the Lorentz term

Build a 2D channel [0,1]², hold the top/bottom walls still, and drive the flow with a pressure-gradient body force in the presence of a transverse field.

n = 32
mesh = grid_triangle_mesh(LevelSetGrid(n, n, 1.0 / (n - 1), 1.0 / (n - 1)))

μ, σ = 1.0, 1.0
B0 = [0.0, 0.0, 2.0]      # transverse field (out of plane)
G = 1.0
N = size(mesh.vertices, 2)

fixed_u = Int[]
for i in 1:N
    y = mesh.vertices[2, i]
    (y < 1e-12 || y > 1 - 1e-12) && append!(fixed_u, 2i - 1, 2i)   # no-slip walls
end

U = FESpace(mesh, ReferenceFE(2, 2); valuetype = VectorValue{2})
loads_u = assemble_vector(LinearForm(U, (x, v, ∇v) -> dot([G, 0.0], v)))

res = solve_mhd_stokes(mesh, μ, σ, B0; fixed_u = fixed_u, loads_u = loads_u, fixed_p = [1], order = 2)
(u = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.006536242929566008, 0.002375942780150479, 0.0014837217731477378, 0.003366858490333126, 0.007191278508969864, 0.004808246212999586, 0.010112910982309018, 0.005672767709562375, -0.0007785047449817965, 0.006780876704240938], p = [0.0, -0.6971470635907582, -0.2908276025263303, -0.18138063966160442, -0.12182313312808599, -0.07841951395850197, -0.050345280890192315, -0.031380364631157016, -0.018599022436997695, -0.010152041717906037  …  0.01107084550894451, 0.019350364930010154, 0.03191340090200722, 0.050665907938820746, 0.07844052422754841, 0.11957083141221392, 0.1893497010092649, 0.2787781947737096, 0.40617497958583976, 1.8526515471962386])

Verification — the Hartmann profile

The centreline velocity must match the analytic value, and it must sit below the field-free Poiseuille limit G/(8μ) — the field has damped the flow.

Ha = (B0[3] / 2) * sqrt(σ / μ)
u_analytic = (G / (σ * B0[3]^2)) * (1 - 1 / cosh(Ha))
center = argmin((mesh.vertices[1, i] - 0.5)^2 + (mesh.vertices[2, i] - 0.5)^2 for i in 1:N)
ux = res.u[2 * center - 1]

println("centreline u  = ", round(ux; digits = 5))
println("Hartmann u    = ", round(u_analytic; digits = 5))
println("Poiseuille    = ", round(G / (8μ); digits = 5))
centreline u  = 0.09797
Hartmann u    = 0.08799
Poiseuille    = 0.125
vtk(mesh; field = [res.u[2i-1] for i in 1:N], fieldname = "u_x (Hartmann)")

Next steps