A surface-tension model that fails at high Reynolds number

Finding it took six experiments. Reporting it took two issues. Fixing it took the maintainers a paper.

In 2023 I was simulating supercooled large droplets hitting an aircraft surface — the physics behind in-flight icing certification. I used SPHinXsys, a well-regarded open-source smoothed-particle-hydrodynamics library, and the droplets kept coming apart.

The obvious conclusion was that I had set it up wrong. That turned out to be false, and establishing it was false is the actual work.

Six experiments, each killing one explanation

A bug report that says "it doesn't work" is worth nothing. What a maintainer needs is the smallest case that fails and the list of things it isn't. So I built a sequence where each step removed a candidate cause.

Setup	Result
3D droplet impact, Re 7154, We 259	Droplet breaks up instead of spreading
Square droplet relaxing to a circle, Re ≈ 1	Correct — a clean circle at rest
The same at Re ≈ 10³	Every fluid particle gone by t ≈ 0.4
The same with maximum numerical dissipation	Settles, but drifts off-centre and particles pass through the wall
Capillary oscillation against a published benchmark	Diverges even at low Re, once the interface is moving
2D impact, dissipative solver	Disintegrates; particles escape a domain 20× the droplet radius

The fifth line is the one that decides it. If the failure appears at low Reynolds number as soon as the interface has velocity, it cannot be a high-Re dissipation problem. It is how curvature is computed on a fast-moving interface.

The benign case: a square patch of water relaxing to a circle under surface tension at low Reynolds number, where the shipped model behaves correctly. Raising the density by three orders of magnitude — changing nothing else — makes the particles disappear.

Mass-centre position against time at three particle resolutions, showing the curves collapsing onto each other — The same oscillation at 900, 3600 and 14 400 particles. The three collapse onto one another, so the divergence from the benchmark is not a discretisation artefact — refining the mesh does not make it go away.

Why nobody had hit it

A failure this stark in a widely used library needs explaining. So I went through every example shipped with it and tabulated the Reynolds and Weber number each one actually exercises.

The high-Re examples are all single-phase, where surface tension never enters. The multiphase examples run at low Re, or drop viscosity and surface tension altogether. Nothing in the suite put high Re and surface tension in the same simulation. The gap sat in untested territory, which is exactly where this kind of thing lives.

This is the part I would want an interviewer to ask about. Not "I found a bug" — anyone can hit a bug. The question worth answering is how you tell the difference between a library being wrong and you being wrong, when the library is used by hundreds of people and you have been doing SPH for eight months.

What I built while chasing it

The leading suspect was the library's Riemann solver: it shipped a linearised, Roe-type solver whose numerical dissipation is implicit, so you cannot easily separate the physics from the damping. I implemented an HLLC approximate Riemann solver — wave-speed estimates, contact-wave state selection, equation-of-state recovery of the intermediate state — and wired it into both the fluid–fluid and fluid–wall integration paths.

It did not fix the problem, which is how I knew the surface-tension model was the culprit rather than the solver. Negative results are cheap to hide and expensive to discard.

What happened next

I reported it upstream in #378 (August 2023, the benchmark I could not reproduce) and #497 (December 2023, the high-Re failure with the parameter study).

The maintainers traced it to zero-surface-energy modes — the surface-tension contributions of symmetrically placed neighbours cancel, so the force is underestimated — and fixed it with a momentum-conserving penalty force, reaching Re = 10,000 and We = 25,000. That work was published in Computer Methods in Applied Mechanics and Engineering 444 (2025), and its acknowledgements read:

"Xiangyu Hu appreciates the discussions on high-speed drop impact with Rui Qiao and Nilotpal Chakraborty from Virginia Tech."

To be precise about what that means: I reported the failure and mapped where it appeared. The mechanism, the remedy and the high-Re results are the authors' work, and I am not an author on that paper. The fix shipped in SPHinXsys v1.2.

The code and data are in a repository I am tidying for release; it will be linked here.