17 Laminar Cylinder Vortex Shedding — Explicit Classical RK4
Physical-time explicit integration on a blockage-controlled mesh, validated against Williamson and Henderson
17.1 Run Provenance and Artifacts
| Item | Value |
|---|---|
| Wake CFD baseline version | 0.22.0 |
| Solver | DBNS |
| Compute device | CPU |
| Parallel execution | 4 MPI ranks |
| Arithmetic precision | FP64 |
| Run date | 2026-08-24 |
This build contains a correction to the strong wall boundary condition on the classical-RK4 path. Classical RK4 keeps a weighted accumulator in addition to the base state, and the no-slip Dirichlet condition was previously applied to the base state only. On earlier builds the wall value drifted, which inflated the friction component of drag while leaving the pressure field — and therefore contour plots — visually unaffected.
Any explicit-RK4 result from a build earlier than 2026-08-23 is not comparable with the numbers below. Implicit dual-time integration is unaffected, because the weighted accumulator exists only on the RK4 path.
17.1.1 Downloadable Files
Not available: standalone input decks and solver result files have not yet been published for download.
17.2 Validation Objective
This case establishes the reference result for laminar cylinder vortex shedding in Wake CFD, using explicit physical-time integration: there is no pseudo-time loop, and the CFL number sets the physical time step directly.
Explicit integration was chosen for the reference because it removes a degree of freedom that implicit dual-time methods carry. A dual-time result depends on how far the inner pseudo-time loop is converged at each physical step, and that convergence is a setting rather than a property of the discretization. Explicit marching has no such setting: each step is the scheme, and the only numerical parameter is the step size itself.
Classical RK4 is fourth-order in time and introduces no amplitude damping at the step sizes used here, so the reported lift amplitude and shedding frequency are not shaped by time-integration error at the level being compared.
17.3 Reference Solution and Data
The Kármán vortex street at Re_D=100 is a standard unsteady validation case with well-established experimental values for Strouhal number, mean drag, and lift amplitude. At this Reynolds number the wake is laminar and two-dimensional — three- dimensional instability sets in near Re_D≈190 — so a 2D simulation is physically appropriate rather than a simplification.
17.4 Simulation Conditions
| Item | Setting |
|---|---|
| Geometry | Circular cylinder, D=1 m |
| Analysis | Unsteady, compressible laminar Navier–Stokes |
Re_D |
100 |
| Mach | 0.05 (mu=0.2081 is used to match Re_D) |
| Reference state | ISA sea level — P=101325 Pa, T=288.15 K → rho=1.2218, U=17.0387 m/s |
| Total simulated time | 12 s — about 34 shedding cycles |
| Readout window | t = 6→12 s (16 cycles), see Section 18.8.2 |
With explicit physical-time marching the time step is not prescribed; the CFL number determines it every step from the local mesh and flow state. The step therefore drifts slightly as the flow develops — here from 4.3663×10^-6 to 4.2857×10^-6 s, a 1.9% change over the run — which matters when frame times are assigned for animation (Section 18.8.2).
17.5 Mesh
| Item | Value |
|---|---|
| Cells | 52,428 (unstructured triangles with a structured near-wall layer) |
| Domain | 20D upstream, 40D downstream, ±30D lateral |
| Blockage | 1.67% (D / lateral extent) |
| Near-wall first cell | y_1 = 0.002D, growth ratio 1.2, about 13 cells inside the boundary layer |
| Wake resolution | D/35 (near wake), D/23 (2<x<5), D/14 (10<x<20) |
| Maximum cell-size jump in the wake | 5× |
Why blockage is reported. A finite lateral domain accelerates the flow past the cylinder and raises the shedding frequency relative to an unbounded experiment, so blockage is a direct contributor to any Strouhal offset. The commonly used guideline is to keep it below about 2%. Reference values from Williamson and Henderson correspond to an unbounded flow, so a validation report that omits its blockage cannot be checked against them meaningfully.
Abrupt cell-size jumps in the wake are reported for the same reason — they act as partially reflecting interfaces for the shed vortices and damp the wake artificially.
The mesh is generated from a parametric definition with a level parameter, so a sqrt(2)-ratio family for grid convergence can be produced from the same source. Grid convergence itself was not performed in this case and remains outstanding (Section 17.9).
17.6 Boundary and Initial Conditions
| Item | Setting |
|---|---|
| Wall | no-slip, adiabatic |
| Outer boundary | farfield |
Initial field: Not reported.
17.7 Numerical Methods
| Item | Setting |
|---|---|
| Convective scheme | SLAU2, second order, weighted least-squares gradient, no limiter |
| Time integration | explicit classical RK4 (four fixed stages), physical-time marching |
CFL |
4.0 (constant) → dt = 4.2857×10^-6 s |
17.8 Results
| Metric | This case | Literature | Deviation |
|---|---|---|---|
St |
0.1661 | 0.164 (Williamson) | +1.3% |
Mean C_D |
1.3473 ± 0.0067 | 1.35 (Henderson) | −0.2% |
C_L amplitude |
0.3327 | ≈0.33 | +0.8% |
| Drag/lift frequency ratio | 1.9993 | 2 (physical requirement) | −0.04% |
The two spectra at full size, since the peak identification is the basis of the reported Strouhal number:
17.8.1 Drag Decomposition
| Component | This case | Literature |
|---|---|---|
| Pressure | 1.0039 | ≈1.02 |
| Friction | 0.3391 | ≈0.33 |
| Friction share | 25.2% | ≈24% |
The decomposition is the stronger statement, not the total. A mean C_D within 0.2% of the reference can arise from compensating errors in the pressure and friction components. Here both components agree independently, which additionally confirms that the no-slip wall condition is enforced correctly — the wall shear stress is the quantity most sensitive to it.
17.8.2 Development History
Plotting against convective time t* = tU/D rather than seconds makes the record comparable across geometries and flow speeds, and the shedding period is read directly as T* = 1/St = 6.02. The run covers t* = 0→204, about 34 cycles.
| Stage | t* |
Behaviour |
|---|---|---|
| Impulsive start | <1 |
C_D spike of 75 (added mass and starting vortex) |
| Recirculation growth | 1→40 |
C_D falls to 1.10 — symmetric steady wake, no shedding yet |
| Instability growth | 15→70 |
C_L grows exponentially, sigma = 0.1091 per t* |
| Nonlinear saturation | 70→85 |
C_D recovers to 1.35 |
| Limit cycle | >85 |
Fully saturated |
The drag minimum at t*≈40 is lower than the saturated value of 1.347. This is the well-known result that vortex shedding increases mean drag relative to the symmetric steady wake, and it appears here without being sought.
The exponential growth rate sigma corresponds to the linear term of the Landau amplitude equation and can be compared with global linear stability analysis. That comparison was not made here; the value is recorded so that a future Reynolds-number sweep can construct sigma(Re) and extrapolate to the critical Reynolds number.
17.8.3 Flow Field
17.8.4 Readout Methodology
Two choices affect the reported numbers enough to be stated explicitly.
17.8.5 Readout window
The transient was initially assumed to end at t = 4 s. It does not. The window was selected by requiring the statistics to stop changing:
| Skip [s] | Window [s] | Cycles | Mean C_D |
Std. dev. | C_L amp. |
St |
|---|---|---|---|---|---|---|
| 4.0 | 7.99 | 22 | 1.3430 | 0.0132 | 0.3325 | 0.1657 |
| 5.0 | 6.99 | 19 | 1.3467 | 0.0069 | 0.3326 | 0.1660 |
| 6.0 | 5.99 | 16 | 1.3473 | 0.0067 | 0.3326 | 0.1661 |
| 7.0 | 4.99 | 14 | 1.3473 | 0.0068 | 0.3326 | 0.1662 |
| 8.0 | 3.99 | 11 | 1.3473 | 0.0068 | 0.3326 | 0.1662 |
The standard deviation halves at 6 s and the mean is then stable to four decimal places. The reported window is t = 6→12 s. Period-to-period variation of the limit cycle over that window is 0.00%.
17.8.6 Frequency estimation
A raw FFT peak was not accurate enough for this run. With a 6-second window the frequency resolution is df/f = 5.9%, so the discrete bin spacing alone bounds the Strouhal error at ±2.9% — the same size as the deviation being reported.
| Method | St |
|---|---|
| FFT maximum bin | 0.1665 |
| Parabolic interpolation of the log-magnitude peak | 0.1659 |
| Upward zero-crossing (adopted) | 0.1661 |
Zero-crossing is not limited by bin spacing. It agrees with the independent parabolic estimate to within 0.1%, and it is the estimate for which the drag/lift frequency ratio comes out as 1.9993 against the physical requirement of exactly 2 — a self-check on the method itself.
The bound scales as 1/(2·T_window), so the distinction matters less for longer runs — a 32-second window would bound the same error at 0.6%. It is stated here because this run traded window length for wall-clock time, and the resulting bound is not negligible against the deviation being reported.
17.9 Conclusions and Applicable Scope
- Explicit classical RK4 reproduces all three reference quantities within 1.3% on a blockage-controlled mesh, with the readout window and frequency estimator both shown to be converged. This is the reference result for this flow in Wake CFD.
- The pressure and friction components of drag agree with the reference separately, which is a stronger statement than agreement of the total and confirms correct enforcement of the no-slip wall.
- The residual
+1.3%Strouhal offset is not resolved into its causes. Three candidates remain open: compressibility atM=0.05(the reference values are incompressible), the remaining 1.67% blockage, and discretization error. Separating them requires a Mach sweep, a domain sweep, and a grid-convergence study respectively.
What this case does not establish.
- Grid convergence was not performed. A
sqrt(2)-ratio mesh family exists in the parametric source but was not run, so discretization error is not bounded. Agreement with experiment at a single grid level can include error cancellation. - The
+1.3%offset is not attributed. No Mach or domain sweep was run. - The linear growth rate is reported, not validated. It is not compared with global stability analysis.
- Results apply to one Reynolds number (
Re_D=100), one mesh, one convective scheme (SLAU2), and one time-integration scheme. They do not transfer to other conditions without repeating the exercise.
17.10 References
- Williamson, C. H. K., “Defining a Universal and Continuous Strouhal–Reynolds Number Relationship for the Laminar Vortex Shedding of a Circular Cylinder,” Physics of Fluids, 31(10), 1988, pp. 2742–2744.
- Henderson, R. D., “Details of the Drag Curve near the Onset of Vortex Shedding,” Physics of Fluids, 7(9), 1995, pp. 2102–2104.