Case definition#

The Reynolds number is

Expected flow physics#

At $Re=100$, flow past a circular cylinder is expected to be laminar but unsteady. The symmetric steady wake is unstable, and the developed solution should exhibit periodic vortex shedding known as the von Kármán vortex street.

Governing equations#

The incompressible Navier–Stokes equations are

The LBM solver recovers these equations in the low-Mach hydrodynamic limit.

D2Q9 lattice Boltzmann formulation#

The D2Q9 discrete velocities are

The weights are

The equilibrium distribution is

Macroscopic density and velocity are reconstructed from

The physical-to-lattice conversion gives

TRT collision model#

The distribution and equilibrium functions are split into even and odd parts with respect to the opposite lattice direction $\bar{i}$:

The TRT collision step is

The even relaxation rate controls viscosity:

The odd rate is determined from the TRT magic parameter $\Lambda$:

Using

produces approximately

Why TRT is used here

TRT separates the viscous even modes from the odd kinetic modes. This significantly improves robustness for halfway bounce-back when $\tau$ is close to $0.5$.

1
let omega_plus = self.cfg.omega_plus();
2
let omega_minus = self.cfg.omega_minus();
3

4
for k in 0..Q {
5
    let ko = OPP[k];
6
    let f_plus = 0.5 * (self.f[base + k] + self.f[base + ko]);
7
    let f_minus = 0.5 * (self.f[base + k] - self.f[base + ko]);
8
    let eq_plus = 0.5 * (eq[k] + eq[ko]);
9
    let eq_minus = 0.5 * (eq[k] - eq[ko]);
10

11
    f_post[k] = self.f[base + k]
12
        - omega_plus * (f_plus - eq_plus)
13
        - omega_minus * (f_minus - eq_minus);
14
}

Boundary conditions#

1
let bi = self.base(nx - 2, j);
2
let rho_inner: f64 = (0..Q).map(|k| self.f[bi + k]).sum();
3
let scale = 1.0 / rho_inner;
4

5
for k in 0..Q {
6
    self.f[b + k] = self.f[bi + k] * scale;
7
}

Symmetry-breaking perturbation#

A perfectly symmetric grid, centered cylinder, symmetric boundaries, and $v=0$ initial field can preserve the unstable symmetric numerical branch. A small perturbation is therefore introduced behind the cylinder.

Version 3 uses a localized streamfunction perturbation:

Velocity perturbations are obtained from

with the code sign convention chosen consistently. Because the perturbation comes from a streamfunction,

This breaks wake symmetry without creating a strong compressibility pulse.

1
let g = (-((x - x0).powi(2) + (y - y0).powi(2))
2
    / (2.0 * sigma2))
3
    .exp();
4
let psi = amp * sigma * g;
5
let du = -(y - y0) * psi / sigma2;
6
let dv =  (x - x0) * psi / sigma2;
7
let ux = u_lu + du;
8
let uy = dv;

Force coefficients#

The momentum-exchange method calculates the force on the cylinder during bounce-back.

In lattice units, the denominator is

Version 3 averages force over each output_every interval:

This reduces high-frequency lattice noise without suppressing the physical shedding frequency.

Running the solver#

1
rm -rf results_test
2

3
cargo run --release -- \
4
  --out results_test \
5
  --t-end 4 \
6
  --output-every 100 \
7
  --print-every 500 \
8
  --frame-every 500 \
9
  --diagnostic-every 100 \
10
  --perturb 1e-3

1
rm -rf results
2

3
cargo run --release -- \
4
  --out results \
5
  --output-every 100 \
6
  --print-every 1000 \
7
  --frame-every 1000 \
8
  --diagnostic-every 100 \
9
  --perturb 1e-3

Output files#

1
case_summary.txt
2
field.csv
3
field_tecplot.dat
4
force_coefficients.csv
5
probe_wake.csv
6
wake_centerline_y5.csv
7
vertical_profile_x10.csv
8
frames/speed_00000.bmp
9
frames/speed_00001.bmp
10
frames/vorticity_00000.bmp
11
frames/vorticity_00001.bmp
12
images/speed.svg
13
images/u_velocity.svg
14
images/v_velocity.svg
15
images/pressure.svg
16
images/vorticity.svg
17
images/streamfunction.svg
18
images/force_coefficients.svg
19
images/wake_v_probe.svg

Animation workflow#

The signed vorticity field is the clearest visualization of alternating shedding.

1
ffmpeg -framerate 10 -i results/frames/vorticity_%05d.bmp \
2
  -c:v libx264 -pix_fmt yuv420p vorticity.mp4

Speed animation:

1
ffmpeg -framerate 10 -i results/frames/speed_%05d.bmp \
2
  -c:v libx264 -pix_fmt yuv420p speed.mp4

Frame interpretation

speed_00000.bmp and vorticity_00000.bmp represent the initial field. Subsequent files represent intervals of frame_every solver steps.

Plots#