Halving the mesh made the blow-up arrive twice as fast — complex eigenvalues in the 6-equation two-fluid model

Halving the mesh made the blow-up arrive twice as fast#

When a solver blows up, the usual first move is to refine the mesh. If the error drops, the problem is in the discretization. If nothing changes, the problem is in the physical model. That is roughly the order most of us work through.

Sometimes the needle moves the other way. Halve the cell size and the blow-up arrives exactly twice as early. Quadruple the cell count and it arrives four times as early. Shrinking the time step does nothing to the growth rate.

Push sigma from 0 to 1.1 and the dots in the left half-plane drop onto the real axis. Drop slip u_r to 0 and the disturbance stops growing.

The dots turn green at the same instant the right-hand panel stops amplifying. The finer-mesh lane hits the ceiling first. The numbers come later.

That symptom is not a discretization bug. It says the governing equations are ill-posed as an initial-value problem — the state where a perturbation grows at a rate inversely proportional to its wavelength, without bound. When Pandare and Luo built a density-based finite-volume two-fluid solver in their 2018 AIAA paper, this was the first thing they had to deal with.

The eigenvalues of the single-pressure six-equation model were pulled directly, to see where the complex pair comes from and what coefficient the interfacial-pressure term needs before the pair sits on the real axis again. The answer is exactly 1.

A pair of slow waves leaves the real axis#

The two-fluid model treats both phases as interpenetrating continua. Mass, momentum and energy are solved separately for each phase. Collapsing the two pressures into one ($p_g = p_l \equiv p$) leaves six PDEs. That is the Wallis model, also called the single-pressure six-equation model.

Switch off compressibility for a moment, quasi-linearize in one dimension with the primitive variables $(\alpha_g, p, u_g, u_l)$, and the pair of slow waves has a closed-form eigenvalue.

Here $\alpha_k$ is the volume fraction, $\rho_k$ the phasic density, $u_k$ the phasic velocity, and $\sigma$ the coefficient of the interfacial-pressure term introduced below. The first term is a density-weighted mean velocity; the second is how far the two waves separate.

Everything hangs on what sits under the radical. For $\sigma < 1$ it goes negative and the two eigenvalues become a complex-conjugate pair. As long as the slip $|u_l - u_g|$ is nonzero, this happens without fail. The moment the phases move at different speeds, the model is ill-posed. Zero slip wipes the radicand out.

One table — seven equations, six, and six with a correction#

Three options compete for this spot. Standing them side by side shows what each one buys and what it sells.

The seven-equation model gives the volume fraction its own transport equation. That buys hyperbolicity, at the cost of pressure-relaxation terms that bring stiffness with them. I laid out that structure in an earlier post on flux splitting for Baer–Nunziato. The catch is that the model is physically justified mainly for densely packed particles and suspensions. It fits a stratified air-water pipe far less well.

The 4×4 eigenvalues, from Python#

Before trusting the closed form, solve the original system as it stands. Build $A W_t + B W_x = 0$ with compressibility retained and take the eigenvalues of $A^{-1}B$. Air and water, $\alpha_g = 0.5$, with the gas running 10 m/s ahead.

import numpy as np
 
def interfacial_dp(a, rg, rl, ur, sigma):
    """Stuhmiller correction: p_int = p - dp"""
    return sigma * a * (1 - a) * rg * rl * ur**2 / (a * rl + (1 - a) * rg)
 
def two_fluid_matrices(a, rg, rl, cg, cl, ug, ul, sigma):
    """A W_t + B W_x = 0,  W = (alpha_g, p, u_g, u_l)"""
    dp = interfacial_dp(a, rg, rl, ul - ug, sigma)
    kg, kl = a / (rg * cg**2), (1 - a) / (rl * cl**2)
    A = np.array([[ 1.0, kg,  0.0,    0.0],
                  [-1.0, kl,  0.0,    0.0],
                  [ 0.0, 0.0, a * rg, 0.0],
                  [ 0.0, 0.0, 0.0,    (1 - a) * rl]])
    B = np.array([[ ug,  ug * kg, a,          0.0],
                  [-ul,  ul * kl, 0.0,        1 - a],
                  [ dp,  a,       a * rg * ug, 0.0],
                  [-dp,  1 - a,   0.0,        (1 - a) * rl * ul]])
    return A, B
 
def char_speeds(sigma, a=0.5, rg=1.2, rl=1000.0, cg=340.0, cl=1500.0, ug=10.0, ul=0.0):
    A, B = two_fluid_matrices(a, rg, rl, cg, cl, ug, ul, sigma)
    return np.linalg.eigvals(np.linalg.solve(A, B))
 
print("air/water, alpha_g=0.5, u_g=10, u_l=0 m/s")
print("sigma   max|Im lambda|   slow pair Re")
for s in [0.0, 0.5, 0.9, 1.0, 1.1, 1.5]:
    lam = char_speeds(s)
    slow = np.sort(lam.real)[1:3]
    print("%5.2f   %12.5f   %8.4f %8.4f" % (s, np.abs(lam.imag).max(), slow[0], slow[1]))

air/water, alpha_g=0.5, u_g=10, u_l=0 m/s
sigma   max|Im lambda|   slow pair Re
 0.00        0.34614     0.0120   0.0120
 0.50        0.24481     0.0120   0.0120
 0.90        0.10967     0.0120   0.0120
 1.00        0.00718     0.0120   0.0120
 1.10        0.00000    -0.0972   0.1212
 1.50        0.00000    -0.2326   0.2566

At $\sigma = 0$ the imaginary part is 0.346 m/s. The real parts of the two slow waves sit together at 0.0120. The gas moves at 10 m/s, yet the wave crawls at 0.012 m/s, because of density weighting: water is 830 times heavier than air, so the mean is dragged onto the liquid.

Raising $\sigma$ shrinks the imaginary part. At 1.1 it hits zero and the pair splits to $-0.097$ and $0.121$. The residual 0.00718 at $\sigma = 1.0$ comes from compressibility. The closed form was derived in the incompressible limit, so a finite speed of sound nudges the threshold a hair above 1.

Write the density-weighted mean out again and 0.0120 stops looking mysterious. With $\alpha_g = 0.5$, $\rho_g = 1.2$, $\rho_l = 1000$ and $u_g = 10$, the denominator is $0.5\times 1.2 + 0.5\times 1000 = 500.6$. The numerator is $0.5\times 1.2\times 10 = 6$, so the mean speed is $6/500.6 = 0.0120$. Water is $1000/1.2 \approx 833$ times heavier than air, and that ratio is the whole division. The gas runs at 10 m/s. The real part of the slow pair still sits at 0.012 m/s.

The imaginary parts follow the closed form just as tightly. The ratio of 0.34614 at $\sigma = 0$ to 0.24481 at $\sigma = 0.5$ is $\sqrt{2} \approx 1.414$. At $\sigma = 0.9$ the factor is $\sqrt{0.1}$ and the imaginary part is 0.10967. $|\mathrm{Im}\,\lambda|$ scales as $\sqrt{1-\sigma}$. At $\sigma = 1.1$ the real parts split to $-0.0972$ and $0.1212$, and the midpoint is still 0.012. The half-width 0.109 matches the imaginary part at $\sigma = 0.9$. Only the sign inside $(\sigma-1)$ flipped.

The closed form puts the threshold at exactly 1#

Now measure the same numbers from the closed form, and bisect for the critical $\sigma$.

from math import sqrt, pi
 
def material_pair(a, sigma, rg=1.2, rl=1000.0, ug=10.0, ul=0.0):
    """Incompressible-limit closed form for the two slow (material) waves"""
    al = 1.0 - a
    den = al * rg + a * rl
    mean = (al * rg * ug + a * rl * ul) / den
    disc = (sigma - 1.0) * a * al * rg * rl * (ul - ug) ** 2 / den**2
    if disc >= 0.0:
        return (mean - sqrt(disc), mean + sqrt(disc)), 0.0
    return (mean, mean), sqrt(-disc)
 
print("closed form vs the 4x4 eigenvalues above")
for s in [0.0, 0.5, 0.9, 1.1, 1.5]:
    (r1, r2), im = material_pair(0.5, s)
    print("sigma=%4.2f  Re = %8.4f %8.4f   |Im| = %8.5f" % (s, r1, r2, im))
 
print()
print("growth rate of the shortest resolved mode, L = 1 m, sigma = 0")
_, im0 = material_pair(0.5, 0.0)
for n in [50, 100, 200, 400, 800]:
    k = pi * n          # k = pi / dx, dx = 1/n
    print("N=%4d  dx=%7.5f  k=%8.1f 1/m  growth=%8.2f 1/s" % (n, 1.0 / n, k, k * im0))
 
print()
print("critical sigma (incompressible limit) for a few states")
for a in [0.1, 0.5, 0.9]:
    for ur in [1.0, 30.0]:
        lo, hi = 0.0, 5.0
        for _ in range(60):
            mid = 0.5 * (lo + hi)
            _, im = material_pair(a, mid, ug=ur)
            if im > 0.0: lo = mid
            else: hi = mid
        print("alpha_g=%.1f  u_r=%4.1f  ->  sigma_c = %.6f" % (a, ur, hi))

closed form vs the 4x4 eigenvalues above
sigma=0.00  Re =   0.0120   0.0120   |Im| =  0.34599
sigma=0.50  Re =   0.0120   0.0120   |Im| =  0.24466
sigma=0.90  Re =   0.0120   0.0120   |Im| =  0.10941
sigma=1.10  Re =  -0.0974   0.1214   |Im| =  0.00000
sigma=1.50  Re =  -0.2327   0.2566   |Im| =  0.00000
 
growth rate of the shortest resolved mode, L = 1 m, sigma = 0
N=  50  dx=0.02000  k=   157.1 1/m  growth=   54.35 1/s
N= 100  dx=0.01000  k=   314.2 1/m  growth=  108.70 1/s
N= 200  dx=0.00500  k=   628.3 1/m  growth=  217.40 1/s
N= 400  dx=0.00250  k=  1256.6 1/m  growth=  434.79 1/s
N= 800  dx=0.00125  k=  2513.3 1/m  growth=  869.58 1/s
 
critical sigma (incompressible limit) for a few states
alpha_g=0.1  u_r= 1.0  ->  sigma_c = 1.000000
alpha_g=0.1  u_r=30.0  ->  sigma_c = 1.000000
alpha_g=0.5  u_r= 1.0  ->  sigma_c = 1.000000
alpha_g=0.5  u_r=30.0  ->  sigma_c = 1.000000
alpha_g=0.9  u_r= 1.0  ->  sigma_c = 1.000000
alpha_g=0.9  u_r=30.0  ->  sigma_c = 1.000000

Leave sigma at 0 and the fine-mesh lane hits the blow-up line first. Push it past 1 and all four lanes lie down together. The 54.35 1/s at N = 50 and the 869.58 1/s at N = 800 split from the same start.

The closed form matches the 4×4 eigenvalues to three decimals. Sweep the volume fraction from 0.1 to 0.9 and the slip from 1 to 30 m/s, and the threshold stays at 1.000000. That is why the coefficient in Stuhmiller's correction,

is not an arbitrary tuning knob. $\sigma = 1$ is the smallest coefficient that drives the radicand to zero. Production codes leave some margin and use a value slightly above 1.

Instability and ill-posedness are different animals#

A numerically unstable scheme improves when the time step shrinks. An ill-posed problem does not, because the growth rate scales with the wavenumber.

Halve $\Delta x$ and the shortest representable wavelength halves with it, so the growth rate doubles. In the output above, 54.35 1/s at $N = 50$ becomes 869.58 1/s at $N = 800$ — exactly a factor of 16. The finer the mesh, the faster the answer dies. The repair belongs to the closure, not to the mesh width.

The table multiplies the same way. At $N = 50$, $k = 157.1$ 1/m and $|\mathrm{Im}\,\lambda| = 0.34599$, so $157.1 \times 0.34599 \approx 54.35$ 1/s. Double $N$ and $k$ doubles with it: 54.35, 108.70, 217.40, 434.79, 869.58, each exactly twice the last. From $N = 50$ to $N = 800$ the mesh is 16 times finer and the growth is 16 times larger. A disturbance that amplifies 54 times a second at $\Delta x = 0.02$ m is near 870 times a second at $\Delta x = 0.00125$ m. Cutting the time step leaves the wavenumber 157.1 alone. What shrinks is $\Delta t$, not $k$.

Real codes often hide the symptom. First-order upwinding contributes numerical diffusion of order $O(k^2 \Delta x)$, which cancels the growth well enough that the run limps along. That is why a code that behaves at low order explodes the moment you raise the order. Numerical diffusion had been quietly paying a debt the model owed.

The rest of the table — how a density-based solver survives low Mach#

Recovering hyperbolicity is not the end of it. Most real multiphase applications sit at very low Mach number, where a density-based solver is chained to the acoustic CFL and the time step collapses.

That territory has traditionally belonged to pressure-based methods. Assuming a solenoidal velocity field erases the speed of sound from the equations, so the CFL depends only on the flow velocity. The price is that compressibility can never be treated rigorously. Bring in a high-temperature phenomenon such as boiling and the errors grow.

Pandare and Luo take a different route: keep the density-based formulation but transform to the primitive variables $[p, v, T]$ and solve fully implicitly. Making pressure an unknown conditions the system far better at low Mach. Interfacial force terms such as drag and virtual mass are also treated implicitly, relaxing the time-step limits further.

The flux side carries the same kind of compromise. When a strong shock meets a material interface, AUSM$^+$-up produces negative pressures. The established remedy was to call an exact Riemann solver just at those faces, but the Newton iterations are expensive. The paper instead adds a volume-fraction coupling term to the mass flux, buying the same robustness — essentially Lax–Friedrichs-type dissipation proportional to the volume-fraction jump. The requirement that a stationary interface must not be disturbed is the same constraint as pulling a smeared interface back with anti-diffusion every step.

Do not raise the order while σ sits below 1#

If a single 4×4 with leftover slip already shows an imaginary part, halving the mesh is a diagnosis, not a fix. A run whose growth goes from 54.35 1/s at $N = 50$ to 869.58 1/s at $N = 800$ is sick in the model, not in the scheme. Read $\sigma$ first.

Below 1, the run is living on the $O(k^2 \Delta x)$ numerical diffusion of first-order upwind. Raise the order and that cushion drops out. Only after $\sigma$ sits near 1.1 does mesh refinement start cutting the error instead of feeding it. Once a density-based solver has survived the low-Mach regime, the next choice is whether to build the pressure explicitly, without a Poisson solve.