Giving a closed solver a feature its API does not have
A ~1,200-line numerical library living inside a commercial CFD code, written because the things it needed were not there.
My dissertation models what happens when a gas turbine ingests sand: the particles hit the blades, weaken, and eventually shatter. Modelling the shattering means a comminution-based breakage model — the first time one has been put inside turbomachinery simulation software. Getting it there is where the engineering is.
The solver is closed. You get a user-subroutine hook, a fixed calling convention, and whatever the API hands you. Three things I needed did not exist — and one thing I needed was forbidden.
1. There is no maths library
The fragment-size statistics need an incomplete beta function. The solver's user-routine environment provides no such thing, so I implemented what the model needed from scratch: log-gamma, the continued-fraction incomplete beta function, and an inverse error function. Alongside them, cubic spline interpolation and spline inversion, to map between cumulative breakage distributions and fragment size classes in both directions.
2. There is no per-particle state
This is the interesting one.
Real particles weaken. An impact too weak to break a grain still damages it, and the next impact meets a weaker grain. To model that, each particle has to carry its accumulated damage with it — across a Lagrangian track, across MPI rank boundaries, for the whole simulation.
The solver carries one float per particle: the diameter. There is nowhere to put anything else.
So I put the damage inside the diameter. The accumulated fracture energy is log-quantised to 14 bits and packed into the low mantissa bits of the IEEE-754 diameter word, behind a flag bit, with guards for denormals, infinities and sign, and graceful fallback to an unmodified value if anything looks wrong. The diameter stays accurate to far more precision than the physics needs; the spare bits carry the history.
Without this, sub-threshold impacts do nothing and every particle meets each impact as good as new. The state channel is what makes history-dependent weakening possible at all inside a closed solver — not faster, not tidier: possible.
3. The return slot is a fixed width
When a particle breaks, the API accepts a fixed six-child return. A real fragment distribution does not have six bins. So the library does conservative rebinning: adaptive bin edges, mass- and number-conserving redistribution, cutoff filtering, renormalisation and adjacent-bin merging, mapped onto those six slots with inactive-slot handling — such that mass and number are conserved regardless of how the distribution actually falls.
A second constraint, and the choice it forced
Modelling breakage is only half the point. What the sponsors want to know is where the blade erodes. The solver ships two erosion models — Tabakoff and Finnie — and neither is particle-size dependent.
That is fine until particles start breaking. A fragment is smaller than its parent and erodes differently at the same impact conditions, so a size-independent model gets the answer wrong in precisely the situation the whole project exists to study. The fix is a size-dependent model (Oka), implemented in Fortran and coupled to the tracking engine so it evaluates in real time as particles move.
Except the solver does not allow a custom breakage model and a custom erosion model in the same routine. One of them has to go.
Breakage stayed inside. It has to: a fragment created in stage one must be tracked through stages two and three, and only the solver's own particle tracker can do that. Erosion moved outside — it can be computed afterwards from the impact conditions the tracker already records, and the in-solver breakage actually helps, because particle size at each impact comes along with the impact data.
Which is a clean answer, and it creates the next problem.
Rebuilding an algorithm the vendor does not document
Erosion computed at impact points gives you a scatter of values on a blade surface. That is close to useless to look at — the impacts are sparse and irregular, and no pattern is visible. What you want is erosion rate at the mesh nodes, which is what the solver's own post-processor shows.
Going from one to the other is not interpolation. You have to build a control volume around each node, work out which impact points fall inside it, and sum their contributions. The solver does this internally and does not document how.
So I reconstructed it: per-node control volumes, point-in-volume assignment via KD-tree spatial search, summation to the node. I could not reproduce the vendor's exact algorithm — I do not know their tie-breaking or their volume construction — but the replacement agrees with it in an average sense, and it writes nodal fields in a format the post-processor will import. Erosion patterns became something you could look at and act on.
I would rather say plainly that it matches on average than claim it is the same algorithm. It is not, and anyone who checked would find that out.
Finding an error nobody knew was there
Separately, while validating particle tracks, I found an undocumented particle-tracking error at transient rotor–stator interfaces in the commercial solver. It does not crash and it does not warn. It produces plausible, wrong answers.
That class of bug is the one worth catching, because nothing downstream tells you it happened. It had already reached conclusions in sponsor-facing work, which had to be corrected, and it changed how the group sets up this kind of simulation.
Why this is the work I want to be hired for
There is a difference between writing a user routine that swaps one drag correlation for another, and building a state channel that does not exist so a physical process can be represented at all. The first is a parameter study. The second is a decision about what the software is capable of.
The domain here is turbomachinery. The transferable part is not.