Passive Scalar Transport¶

The passive scalar transport case validates the double-distribution-function (DDF) advection-diffusion module against exact analytical solutions for three complementary regimes:

12a Pure diffusion of a single Fourier mode on a triply periodic domain isolates the scalar collision-streaming kernel.
12b 1D advection-diffusion strip with a uniform inlet validates the ScalarUniformInlet BC plus the classic Peclet-controlled steady-state profile.
12c 1D conduction between two Dirichlet walls validates the ScalarRegularizedDirichlet BC and the linear steady-state profile.

The 12a case has no boundary condition, no inlet turbulence, no body force, and no fluid motion. The 12b and 12c cases exercise the scalar boundary conditions; the per-BC reconstruction (ScalarUniformInlet, ScalarRegularizedDirichlet) produces a finite, physically correct scalar. 12b (01c_scalar_advection_diffusion_strip.ipynb) overlays the exact advection-diffusion profile and shows clean second-order convergence; 12c (01d_scalar_wall_dirichlet_conduction.ipynb) recovers the linear conduction profile to machine precision.

Sinusoidal Scalar Diffusion (Case 12a)¶

A scalar field \(\phi(x, y, z, t)\) is initialised as a single cosine mode along \(x\) on a periodic box \([0, L)^3\), with the fluid at rest (\(u = 0\), \(\rho = \rho_0\)). Because the velocity field is identically zero, the advection term in \(\partial_t \phi + u \cdot \nabla \phi = D \nabla^2 \phi\) vanishes and the scalar evolves by pure diffusion. A single Fourier mode is an eigenfunction of the Laplacian, so the analytical solution is an exponential decay valid for all \(t > 0\).

Initial conditions¶

The scalar field at \(t = 0\) is:

(1)¶\[\phi(x, y, z, 0) = \phi_0 \cos(k x)\]

where \(k = 2\pi / L\) is the wave number, \(L\) is the domain length, and \(\phi_0\) is the amplitude. On a lattice of \(N\) nodes per side with \(\Delta x = 1\), the wave number in lattice units is \(k = 2\pi / N\). The fluid is initialised uniformly with \(\rho = 1\) and \(u = 0\).

The scalar IC is set through scalar_transports[*].initial_field as a SymPy expression in \(x, y, z\); the fluid IC defaults to the rest state and needs no equation initialisation.

Analytical solution¶

For the single-mode initial condition (1) and stationary fluid the diffusion equation admits the closed-form solution:

(2)¶\[\phi(x, y, z, t) = \phi_0 \cos(k x) \exp\!\left(-D k^2 t\right)\]

The amplitude decays exponentially at rate \(\lambda = D k^2\); the spatial profile is a rescaled copy of the IC at every time. The \(L^2\)-norm of the field decays at the same rate:

(3)¶\[\|\phi(\cdot, t)\|_{L^2} = \|\phi(\cdot, 0)\|_{L^2} \, \exp\!\left(-D k^2 t\right)\]

For a well-resolved simulation the solver must reproduce both the spatial profile and the decay rate to within the second-order truncation error of the DDF lattice Boltzmann scheme.

Simulation parameters¶

Four grid resolutions are run to assess spatial convergence. The fluid solves the standard isothermal RR-BGK equation on D3Q27 and stays at the rest equilibrium fixed point (\(u = 0\), \(\rho = 1\)); the scalar uses the minimal D3Q7 lattice (\(c_s^2 = 1/4\)) with the RR-BGK collision operator from the same Hermite hierarchy as the fluid.

The convergence study mirrors the TGV-2D protocol Simulation parameters: the scalar relaxation time is held close to \(1/2\) and follows \(\tau_\phi - 1/2 \propto N\), while \(n_{\text{steps}} \propto N\). This keeps the dimensionless Fourier number \(\mathrm{Fo} = D \, k^2 \, t = D_{\text{lbm}} \, k_{\text{lbm}}^2 \, n_{\text{steps}}\) constant across grids (~4.9% physical decay), so every resolution is compared at the same physical state.

Parameter	Value
Fluid velocity set	D3Q27
Fluid collision operator	RRBGK
Fluid relaxation time \(\tau\)	0.6 (constant; fluid is stationary, only damps numerical noise)
Fluid initial condition	\(u = 0\), \(\rho = 1\)
Scalar velocity set	D3Q7
Scalar collision operator	RRBGK
Scalar relaxation time \(\tau_\phi\)	0.501 - 0.508 (varies with grid, see table below)
Scalar diffusivity \(D\) (lattice units)	\((\tau_\phi - 1/2) c_{s,\phi}^2\), ramps with \(N\)
Scalar amplitude \(\phi_0\)	1.0
Domain	\([0, L)^3\), periodic in all three directions
Grid sizes \(N\)	16, 32, 64, 128 nodes per side
Boundary conditions	Fully periodic (`periodic_dims: [true, true, true]`)
Initial scalar field	`cos(2pix/N)`

The per-grid scaling is summarised below (using \(\tau_\phi - 1/2 = N / 16000\) and \(n_{\text{steps}} = 80 N\)):

\(N\)	\(\tau_\phi\)	\(D\) (lattice units)	\(n_{\text{steps}}\)
16	0.501	2.5e-4	1280
32	0.502	5.0e-4	2560
64	0.504	1.0e-3	5120
128	0.508	2.0e-3	10240

Note

The diffusivity \(D\) is specified directly in lattice units. Any SI conversion is left to a post-processing step on the saved snapshots.

Validation metrics¶

Two complementary metrics are extracted from the saved snapshots and compared against (2):

Spatial profile error. The normalised \(L^2\) error in the scalar field at the final step,

(4)¶\[E_{L^2} = \frac{1}{\phi_0} \sqrt{\frac{1}{N^3} \sum_{i,j,k} \left[\phi^{\text{lbm}}_{i,j,k} - \phi^{\text{exact}}_{i,j,k}\right]^2}\]

is plotted against $\Delta x$ on a log-log axis. A second-order scheme produces a slope of 2.

Decay rate. The \(L^2\)-norm \(\|\phi(\cdot, t)\|_{L^2}\) at every snapshot is fit on a log-axis, and the measured decay rate \(\lambda_{\text{lbm}}\) is compared with the analytical \(\lambda_{\text{exact}} = D k^2\). Both should agree to within the spatial truncation error of the scheme.

1D Advection-Diffusion Strip (Case 12b)¶

A scalar \(\phi\) is advected and diffused along the streamwise \(x\) axis on a 3-D box with periodic \(y, z\). A uniform fluid \(u = (U, 0, 0)\) carries the scalar from a Dirichlet inlet at the west face into a Dirichlet outlet at the east face, producing the textbook 1D advection-diffusion steady-state profile.

Setup¶

Fluid: \(u = (U, 0, 0)\), \(\rho = 1\). West face is a UniformFlow inlet at \((U, 0, 0)\), east face is a RegularizedNeumannOutlet (zero-gradient outflow).
Scalar: ScalarUniformInlet at the west face with \(\phi_{\text{inlet}} = 1\) (the advection velocity \(U\) comes from the fluid inlet); ScalarRegularizedDirichlet at the east face with \(\phi_w = 0\). \(y\) and \(z\) are periodic so the problem is effectively 1D in \(x\).

Analytical solution¶

The steady-state 1D advection-diffusion equation \(U\, d\phi/dx = D\, d^2\phi / dx^2\) with \(\phi(0) = 1\) and \(\phi(L) = 0\) has the closed-form solution

(5)¶\[\phi(x) = \frac{e^{\mathrm{Pe}\, x / L} - e^{\mathrm{Pe}}}{1 - e^{\mathrm{Pe}}}, \qquad \mathrm{Pe} = \frac{U L}{D}.\]

The Peclet number \(\mathrm{Pe}\) controls the curvature: at small \(\mathrm{Pe}\) (diffusion-dominated) \(\phi\) is approximately linear; at large \(\mathrm{Pe}\) (advection-dominated) \(\phi\) stays near 1 across most of the domain and collapses to 0 in a thin layer next to the outlet.

Simulation parameters¶

The convergence study holds \(\mathrm{Pe} = 8\) constant across four grid resolutions, putting the profile in the mixed regime where both the inlet and the outlet boundary layer contribute visibly to the steady-state shape.

Parameter	Value
Fluid velocity set	D3Q27
Fluid collision operator	RRBGK
Fluid relaxation time \(\tau\)	0.6
Streamwise velocity \(U\) (lattice units)	0.005 (constant across grids)
Scalar velocity set	D3Q7
Scalar collision operator	RRBGK
Scalar diffusivity \(D\) (lattice units)	\(N / 1600\), ramps with \(N\)
Grid sizes \(N\)	16, 32, 64, 128 (along \(x\)); 8x8 in \(y, z\)
Peclet number \(\mathrm{Pe} = U L / D\)	8 (constant across grids)
Boundary conditions (fluid)	West: `UniformFlow` \((U, 0, 0)\). East: `RegularizedNeumannOutlet`. \(y, z\) periodic.
Boundary conditions (scalar)	West: `ScalarUniformInlet` \((\phi_{\text{inlet}} = 1)\); the advection velocity comes from the fluid inlet. East: `ScalarRegularizedDirichlet` \((\phi_w = 0)\). \(y, z\) periodic.
Initial scalar field	0 everywhere; inlet drives the profile to steady state

Validation metrics¶

The normalised \(L^2\) error against the analytical profile above is taken at the final step on a 1D slice along \(x\) (the field is uniform in \(y, z\) to within machine precision):

\[E_{L^2} = \sqrt{\frac{1}{N} \sum_{i} \left[\phi^{\text{lbm}}_i - \phi^{\text{exact}}_i\right]^2}.\]

Plotting \(E_{L^2}\) against \(\Delta x = 1 / N\) on log-log axes is expected to give a slope of 2 for the regularized BC pair.

1D Wall Conduction (Case 12c)¶

A scalar \(\phi\) diffuses between two Dirichlet walls in a 3-D box with periodic \(y, z\). The fluid is held at rest (\(u = 0\)) so the dynamics reduce to pure 1D conduction, producing a linear steady-state profile.

Setup¶

Fluid: \(u = 0\), \(\rho = 1\). West and east faces use RegularizedHWBB no-slip walls; \(y, z\) are periodic.
Scalar: ScalarRegularizedDirichlet at the west face with \(\phi_w = 0\), ScalarRegularizedDirichlet at the east face with \(\phi_w = 1\). \(y, z\) are periodic.

Analytical solution¶

The steady-state 1D pure-diffusion equation \(D\, d^2\phi / dx^2 = 0\) with \(\phi(0) = 0\) and \(\phi(L) = 1\) is a linear ramp:

(6)¶\[\phi(x) = \frac{x}{L}.\]

The transient decays at rate \(\lambda = \pi^2 D / L^2\), so the profile reaches its asymptote on a time scale \(t_{\text{ss}} \sim L^2 / D\).

Simulation parameters¶

The case uses a constant scalar diffusivity \(D = 0.05\) (\(\tau_\phi = 0.7\), well-stable) across three grid sizes. The number of steps is scaled with \(N^2\) so the profile reaches steady state on every grid.

Parameter	Value
Fluid velocity set	D3Q27
Fluid collision operator	RRBGK
Fluid relaxation time \(\tau\)	0.6
Fluid initial condition	\(u = 0\), \(\rho = 1\)
Scalar velocity set	D3Q7
Scalar collision operator	RRBGK
Scalar diffusivity \(D\) (lattice units)	0.05 (constant; \(\tau_\phi = 0.7\))
Grid sizes \(N\)	16, 32, 64 (along \(x\)); 8x8 in \(y, z\)
Boundary conditions (fluid)	West, East: `RegularizedHWBB` no-slip. \(y, z\) periodic.
Boundary conditions (scalar)	West: `ScalarRegularizedDirichlet` \(\phi_w = 0\). East: `ScalarRegularizedDirichlet` \(\phi_w = 1\). \(y, z\) periodic.
Initial scalar field	0 everywhere; Dirichlet walls drive the profile to steady state

Validation metrics¶

The normalised \(L^2\) error against (6) is taken at the final step on a 1D slice along \(x\):

\[E_{L^2} = \sqrt{\frac{1}{N} \sum_{i} \left[\phi^{\text{lbm}}_i - i / (N - 1)\right]^2}.\]

A second-order accurate BC reconstruction gives slope 2 on the \(E_{L^2}\) vs. \(\Delta x\) log-log plot. The mid-domain value \(\phi(L/2)\) is also reported and should equal \(0.5\) to within the BC truncation error.

Note

Each sub-case runs end to end with a finite scalar field, and has its own YAML file and validation notebook in this case folder - 01_passive_scalar_transport.nassu.yaml (12a, pure diffusion; 01a/01b notebooks), 01.1_scalar_advection_diffusion_strip.nassu.yaml (12b; 01c_scalar_advection_diffusion_strip.ipynb), 01.2_scalar_wall_dirichlet_conduction.nassu.yaml (12c; 01d_scalar_wall_dirichlet_conduction.ipynb).