Measuring Viscosity from a Poiseuille Profile

The task

Verify that a lattice Boltzmann solver produces the viscosity Chapman–Enskog predicts — \(\nu = c_s^2(\tau - \tfrac12)\) — rather than merely assuming it.

Step 1: the exact solution

Body-forced flow between parallel plates, at steady state:

\[\nu\frac{d^2u}{dy^2} = -g \qquad\Longrightarrow\qquad u(y) = \frac{g}{2\nu}\,y\,(H - y)\]

A parabola, zero at both walls, peaking at \(u_{\max} = gH^2/8\nu\).

Step 2: where is the wall?

The detail that decides whether you get 1% or 10%. With halfway bounce-back, populations reverse midway between the last fluid node and the first solid node. So if fluid rows are \(j = 0\ldots N_y-1\) at positions \(y_j = j + \tfrac12\), the walls sit at \(y = 0\) and \(y = N_y\):

\[H = N_y\]

Getting this wrong by one lattice unit — using \(H = N_y - 1\), say — introduces a systematic error of order \(2/N_y\), which at \(N_y = 21\) is about 10%. It will look like a failure of Chapman–Enskog and it is a failure of bookkeeping.

Step 3: fit, do not eyeball

Rather than reading off \(u_{\max}\) (one point, maximally sensitive to noise), least-squares fit the whole profile to the parabola shape. Writing \(u_j = C\,p_j\) with \(p_j = y_j(H - y_j)\):

\[C = \frac{\sum_j u_j p_j}{\sum_j p_j^2}, \qquad \nu_{\rm measured} = \frac{g}{2C}\]

Every row contributes, so the estimate is far more robust.

Step 4: converge properly

The parabola only holds at steady state. Run until the centreline velocity stops changing:

\[\left|\frac{u_{\rm mid}^{n+500} - u_{\rm mid}^{n}}{u_{\rm mid}^{n}}\right| < 10^{-10}\]

The time to reach steady state scales as \(H^2/\nu\), so low-\(\tau\) runs take longest — exactly the opposite of the intuition that low viscosity means fast. At \(\tau = 0.6\) the 21-row channel needs ~30 000 steps; at \(\tau = 1.5\) it needs ~4 000.

Step 5: the result

Measured on a \(6\times21\) channel with \(g = 10^{-6}\):

\(\tau\)	\(\nu\) theory	\(\nu\) measured	error	steps to converge
0.60	3.3333e−2	3.3423e−2	0.27%	30 000
0.70	6.6667e−2	6.6816e−2	0.22%	15 500
0.80	1.0000e−1	1.0015e−1	0.15%	11 000
1.00	1.6667e−1	1.6651e−1	0.09%	7 000
1.20	2.3333e−1	2.3227e−1	0.46%	5 000
1.50	3.3333e−1	3.2928e−1	1.21%	4 000

Step 6: read the residual

The error is not random — it is smallest near \(\tau = 1\) and grows at both ends. That structure is the physics of the test, and worth reading:

Rising \(\tau\): halfway bounce-back only places the wall exactly midway in the small-\(\tau\) limit. BGK develops a viscosity-dependent slip, so the effective \(H\) drifts and the fitted \(\nu\) with it. This is documented behaviour, and the reason MRT collision operators — which fix the slip — exist.
Falling \(\tau\) toward 0.5: \(\nu \to 0\), convergence takes forever, and the finite forcing raises the Mach number relative to the shrinking viscous scale.

A result that moves systematically when you change a numerical parameter is telling you something. Here it is telling you about your boundary condition, and the correct response is to report the dependence rather than to quote the best number.

Step 7: keep the Mach number down

Peak \(u/c_s\) stays around \(10^{-3}\) in the table above — safely incompressible. Raise the forcing until \(u/c_s \gtrsim 0.1\) and the measured \(\nu\) starts to drift for a different reason: the truncated equilibrium introduces spurious compressibility that no grid refinement removes.

What this exercise establishes

That the fluid is emergent. Nothing in the code computes a viscous stress; there is only streaming, local relaxation and a body force. Navier–Stokes — with the right viscosity, to three digits — falls out of that loop. It is the most convincing demonstration of Chapman–Enskog available, and it takes an afternoon.

Chapman–Enskog · D2Q9 lattice · LBM widget · Verification and validation · Hagen–Poiseuille (PC316)