Two miscible liquids across a step

The miscible-mixing family had nothing behind it. Its solver is interFoam with the interface compression taken out and a diffusivity put in: the two liquids mix instead of holding a surface. That makes its phase fraction obey an advection–diffusion equation, and with the fluid at rest and no turbulence model, pure diffusion — which has an exact solution for a step.

The point of the case is that the exact solution says three separate things, and they fail separately. Where the front is. How wide it is. And how much mixture there is altogether. A run can be right about the first and the third while being wrong by a third about the second, and one of the fixtures below is exactly that.

The phase fraction in this solver obeys

    d(alpha)/dt + div(phi alpha) = laplacian(Dab + alphatab nut, alpha)

With the fluid at rest, φ vanishes; with no turbulence model, nut vanishes; and what is left is the one-dimensional diffusion equation with constant D. A step at y₀ then has the similarity solution

    alpha(y, t) = 0.5 erfc( (y - y0) / (2 sqrt(D t)) )

Three consequences, and the case measures all three:

• The half-height does not move. Pure diffusion is symmetric about the step, so any displacement at all is something advecting the mixture — buoyancy from a density difference, or a spurious current at the interface.
• The width grows as √(Dt). Measured as the distance from α = 0.75 to α = 0.25, it is 1.908 √(Dt), which at t = 2 s is 0.1206564 m. This is the strong check: it is a rate, and it is the only thing on the page that says the diffusivity reaching the solve is the one the case was given. Numerical diffusion and physical diffusivity are indistinguishable in the answer itself.
• The integral of α is constant — to round-off, not to a discretisation error, because the fluid is at rest in a box it cannot leave. A scheme can sit on the erfc curve through the middle of the domain while losing material at a boundary.

The reference is the infinite-domain solution, so it is only usable while the transition stays away from the walls. Both boundaries are 6.3 √(Dt) from the step, where the exact solution differs from 0 and 1 by 3.872 × 10⁻⁶ — three orders below the agreement being measured. That is computed and checked by the harness rather than asserted.

The last two rows are there because the comparison stops meaning anything without them. If α varied across x and z, the profile being compared with a one-dimensional closed form would be an average over cells that disagree, and the agreement reported would be an average of disagreements. And the erfc solution is the pure diffusion one; with the fluid moving there would be an advection term in the equation being solved that is not in the equation being compared against.

One width within a tenth of a per cent of a closed form is a number that happens to be close. A sequence of them approaching it at the scheme’s own order is the claim.

The time step is refined separately on the finest mesh, and it barely moves the answer — which is the statement that what remains is spatial:

The turbulence model is why the width is the measurement. The product’s default for this family is k-omega SST, and the α equation diffuses with Dab + alphatab·nut. Left on, nothing about the run looks unwell: it converges, α stays inside [0, 1], every process exits 0, max|U| is 1.3 × 10⁻¹⁶ m/s so the fluid never moves at all, and the mixture is centred exactly where it started. The only thing that sees it is the width, which comes out 0.156471 m — 29.68 % wide. The L1 error also stops falling: 0.02688, 0.02718 and 0.02725 on the three meshes, an observed order of −0.016 and −0.004, so refining the mesh makes it very slightly worse rather than better.

The adaptive time step is why the order is here. With the product’s own adjustTimeStep yes, the width (−0.20 %), the position (0.000 %), the mass (−4.0 × 10⁻⁷) and the velocity all sit comfortably inside every band above, on every mesh. Only the refinement study says anything: the observed order is 1.57 then 0.94, against 1.93 and 1.82 for the same case at a fixed step. The run takes 207 steps at 32 cells and 207 steps at 128 — this solver’s α equation is explicit, so its step is bounded by a diffusion number and not only by a Courant number, and the product’s adaptive controls are convective. At rest they run to maxDeltaT and stay there however fine the mesh gets. Without a refinement study this benchmark would have passed that run and reported a verified solver.

The default densities are the one that passes. The staged phase pair is water over air, a density ratio of 815 held up by nothing, and the benchmark sets them equal so the fluid is genuinely at rest. Run with the defaults instead, every check still passes — width −0.047 %, order 1.97 on both mesh pairs — but max|U| is no longer round-off. It is 1.06 × 10⁻⁴, 1.61 × 10⁻⁴ and 2.70 × 10⁻⁴ m/s on the three meshes: small, but growing under refinement rather than shrinking, and the half-height moves 3.6 µm off the step. Both stay far inside the harness’s limits and the fixture is recorded as passing; it is reported here because a velocity that grows when the mesh is refined is not converging to zero.

• Only the mesh has a convergence study behind it. The time step was varied and the answer moved very little, which is consistent with the spatial error dominating, but no temporal order is claimed from three points where the residual spatial error is the larger term.
• The mass drift is 10⁻⁷, not round-off, and it is not fully explained. The analysis says the integral of α is constant to machine precision. Measured, it is −6.9 × 10⁻⁷, −4.3 × 10⁻⁷ and −3.4 × 10⁻⁷ on the three meshes — seven orders above round-off. It shrinks under refinement rather than staying put, which is how a discretisation loss behaves and not how a leak behaves, and the harness’s tolerance is set against what was measured rather than against what the mathematics promises. But no run here isolates where it comes from.
• One diffusivity, one end time, one geometry. D = 0.002 m²/s at t = 2 s in the box the application writes. The width check is a check on √(Dt) at a single point on that curve.
• The domain is finite and the reference is not. The comparison holds because the transition is 6.3 √(Dt) from each wall, where erfc is 3.9 × 10⁻⁶ from its limits. Run for longer and the reference stops applying — which is a property of this case, not of the solver.

One caveat on the orders themselves. This mesh study holds the time step fixed while the cell size halves, so the diffusion number D Δt/Δx² quadruples at every level — 0.0032, 0.0128, 0.0512 — and the sweep reads that climb as well as the mesh. Repeating it with the diffusion number held at 0.0032 gives L1 orders of 2.019 and 1.994, with an L1 lower by 5.9 % and 16.6 % on the two finer meshes.

So the 1.93 and 1.82 above are what a confounded study measures, and they understate the scheme rather than flatter it. They are left in place because they are what this benchmark’s fixture holds and what its gate judges. The same confound was found in two other benchmarks — in one of them it removes convergence altogether — and is recorded in the project’s gaps ledger.