Multiphase and interfacial flows

A fifth of the space, but how much of the mass?

Imagine a 1 L representative portion of a gas–liquid flow. Gas occupies 0.2 L and liquid occupies 0.8 L. The supplied gas density is 2 kg/m³; the supplied liquid density is 1000 kg/m³. Does 20% gas by volume mean 20% gas by mass? And if the gas travels faster than the liquid, does the stored mass fraction tell us the fraction of mass flowing past a window each second?

The gas in this portion has mass 0.4 g; the liquid has mass 800 g. Gas therefore accounts for about 0.04998% of the stored mass. A picture dominated by visible bubbles can still represent a very small gas mass. Density converts occupied volume into mass; motion is a separate question.

In this supplied model the gas moves at 1.5 m/s and the liquid at 0.5 m/s in the same direction. Across the whole 0.001 m² flow window, gas carries 0.6 g/s while liquid carries 400 g/s. The gas fraction of the mass flow is about 0.1498%, roughly three times its stored mass fraction, because gas moves three times faster. Neither fraction is 20%.

On a first reading, follow the window example, compare the two sphere sizes, and study the drop and evaporation examples. Return to the sections marked “Further” for averaged equations and distributions. The published lessons “Conservation of mass, momentum, and energy” and “Pressure, buoyancy, and fluid statics” prepare the balance and buoyancy reasoning. Essential quantities, including the drag model's dimensionless ratios, are defined here before they are used.

Name the phases and their arrangement

A phase is a region with a distinct material state or composition separated from another region by an interface. Liquid drops in gas, gas bubbles in liquid and solid grains carried by water are multiphase examples. Two immiscible liquids can form two phases even though both are liquid. Dissolved salt in a uniform liquid solution does not by itself create a second moving liquid phase.

The continuous phase forms the surrounding connected fluid; a dispersed phase occurs as separated bubbles, drops or particles. These names describe the arrangement at the chosen scale. In a stratified pipe, two connected layers may be more useful descriptions than continuous and dispersed. In a slug flow, large gas pockets alternate with liquid regions; in annular flow a liquid film and a central gas region can coexist with entrained drops.

Arrangement changes the available interface, wall contact, paths for momentum transfer and opportunities for phases to separate. The same overall gas fraction can describe many small bubbles or one long pocket. Fraction alone therefore selects neither a flow pattern nor a pressure drop. Geometry, orientation, inlet history, flow rates, material properties and wetting conditions can change the outcome.

This lesson follows transport between and within phases. “Surface tension and capillarity” explains static interface work, curvature pressure and wetting. “Fluid instabilities and convection” studies whether a specified small disturbance grows. Neither a static capillary relation nor a linear instability calculation supplies a complete developed two-phase flow.

Separate what is present from what crosses

The volume fraction αₖ is the volume of phase k divided by the total sampling volume. It lies between 0 and 1, and the fractions sum to 1 when the phases fill the volume. For gas and liquid, αG + αL = 1. The sampling scale matters: a tiny point can lie entirely in gas while a larger representative volume contains both phases.

For the window calculation we explicitly assign an axially uniform representative cell, so the gas and liquid area fractions equal their volume fractions. Each phase has a uniform velocity uₖ over its own occupied window area. Its volume rate is Qₖ = αₖAuₖ, where A is the whole window area. The superficial flux jₖ = Qₖ/A = αₖuₖ uses the whole area, not only the area occupied by that phase. It has units m/s but need not be a particle's speed.

Here jG = 0.2 × 1.5 = 0.3 m/s and jL = 0.8 × 0.5 = 0.4 m/s. Multiplication by the whole area gives 0.3 L/s of gas and 0.4 L/s of liquid. The gas volume-flow fraction is 0.3/0.7 ≈ 42.86%, larger than its 20% occupied fraction. A moving stream samples both the amount available and its crossing speed.

For constant phase densities, mixture density is ρm = αGρG + αLρL = 800.4 kg/m³. Stored gas mass fraction is YG = αGρG/ρm. Flowing gas mass quality is xG = ṁG/(ṁG + ṁL), with ṁₖ = ρₖQₖ. We use quality for these positive co-current mass flows; signed countercurrent fluxes require explicit conventions and may not behave as fractions between 0 and 1. An arbitrary photograph or an unresolved cross-section does not justify the uniform-cell assumptions.

Two means do not erase slip

Slip is relative phase motion: uG − uL = 1.0 m/s in the supplied cell. A single mean can hide that motion. The total volume flux j = αGuG + αLuL is 0.700 m/s. It is the volume-weighted mean for this uniform case. The mass-average, or barycentric, velocity um = (αGρGuG + αLρLuL)/ρm is about 0.500500 m/s, close to the liquid speed because liquid carries nearly all the mass.

The gas drift flux relative to j is jG − αGj = αG(uG − j). It is 0.160 m/s here, with an equal and opposite liquid drift flux. It describes excess gas volume transport relative to a frame moving at j. The gas's relative velocity uG − j is 0.800 m/s; multiplying by αG gives a flux, so the two quantities should not share a label.

Further: the axial convective momentum flux per total area is Π = αGρGuG² + αLρLuL² = 200.900 Pa. Writing only ρmum² gives about 200.500200 Pa. Expansion of each uₖ = um + (uₖ − um) gives the exact identity Π = ρmum² + D, where D = Σₖαₖρₖ(uₖ − um)² ≈ 0.399800 Pa. The cross term vanishes because um is mass weighted. D is nonnegative in this one-direction calculation.

This is a momentum-transport identity, not a new equation of state or a computed pressure. Pressure and viscous stress were excluded from Π. Even with no fluctuations inside either phase, interphase slip adds transport beyond the mean term. More general averaging also introduces correlations within each phase; replacing all of them by a guessed mixture viscosity is an additional model, not an algebraic consequence.

The same amount can expose twice the interface

Imagine a separate 1 L sample containing a dispersed volume fraction of 0.02 as equal, nonoverlapping spherical drops. A sphere of diameter D has volume πD³/6 and area πD². Dividing the total dispersed volume by one sphere's volume gives the representative number of drops. Multiplying that number by one sphere's area gives interfacial area per mixture volume ai = 6αd/D, in m²/m³, or m⁻¹.

At D = 1 mm, ai = 120 m⁻¹, so this 1 L sample has 0.120 m² of interface. At D = 0.5 mm, ai = 240 m⁻¹ and the total interface is 0.240 m². Halving diameter requires eight times as many spheres to preserve dispersed volume but doubles total interface. These are geometric inventories, not predictions that a mixer will produce those sizes.

Heat or species transfer often depends on both available area and transfer per unit area. More interface can therefore matter even when phase amounts are unchanged. It does not guarantee twice the device performance: gradients, residence time, circulation, crowding and transfer coefficients can also change. Breakup creates new interfaces and coalescence joins them; their rates require dynamics and material conditions.

Further: for spherical drops with number distribution n(D), αd = (π/6)∫D³n(D)dD and ai = π∫D²n(D)dD. Define the surface–volume, or Sauter, diameter D32 = ∫D³n(D)dD / ∫D²n(D)dD. Then ai = 6αd/D32. An arithmetic mean diameter generally cannot replace D32. For strongly distorted interfaces, the spherical identity itself must be replaced by measured or resolved area.

A material boundary and a phase-changing front differ

An interface is a place to impose conservation and mechanical conditions. Let n point from liquid L into vapor V, and let VΓ be the interface's normal speed. Define m″ as mass crossing unit interfacial area per second, positive from liquid to vapor. With no mass stored on the ideal zero-thickness interface, m″ = ρL(uL·n − VΓ) = ρV(uV·n − VΓ). This is the same mass rate seen from the moving boundary on both sides.

If there is no mass transfer, m″ = 0 and both adjacent normal fluid velocities equal VΓ. Such an interface moves with the material in the normal direction. For ordinary viscous fluids without an interfacial slip model, the tangential velocities also match, although their velocity gradients and stresses need not match separately. A clean liquid interface is consequently different from the no-slip surface of a rigid particle.

Traction is the force per area transmitted by stress across a surface. Its normal balance includes curvature and surface tension; its tangential balance includes gradients of surface tension along the interface. Uniform surface tension removes that gradient term, not the interface itself. Temperature or composition variations can create Marangoni stresses: unequal tangential pulls that couple surface transport to bulk motion.

A moving or phase-changing interface therefore needs more than a picture of its position. Depending on the question, it needs stress conditions, thermal or species conditions, material properties and boundary data. The sign of a stress jump depends on the chosen normal and jump convention. The phase-change mass convention above will also fix the signs of the next heat balance.

A clean drop can settle faster than a rigid sphere

Consider a separate isolated liquid drop of radius a = 50 μm settling steadily downward in a quiescent liquid. Assign drop density 1050 kg/m³, surrounding density 1000 kg/m³, both dynamic viscosities 0.001 Pa·s, surface tension 0.030 N/m and g = 9.81 m/s². The domain is unbounded; the interface is clean with uniform surface tension and no mass transfer. This is a small, nearly spherical creeping-flow model.

The downward gravitational force minus buoyancy is (4π/3)a³(ρd − ρc)g. For a rigid sphere, balancing that force with 6πμcaU gives US = 2(ρd − ρc)ga²/(9μc) = 0.2725 mm/s. The subscript c denotes the continuous outer liquid. A clean fluid drop has internal motion and different interfacial stress conditions. Its Hadamard–Rybczynski result is UHR = US × 3(1 + λ)/(2 + 3λ), where λ = μd/μc.

Here λ = 1, so UHR = 1.2US = 0.3270 mm/s. The clean-drop drag coefficient in this force law is 2πμca(2 + 3λ)/(1 + λ), giving 5πμca at λ = 1, compared with 6πμca for the rigid reference. Both terminal states balance the same net body force, about 2.568 × 10⁻¹⁰ N; the smaller drag at a given speed permits the larger terminal speed.

Check the approximation with diameter D = 2a = 100 μm. The clean result has ReD = ρcUHRD/μc = 0.0327, WeD = ρcUHR²D/σ ≈ 3.56 × 10⁻⁷ and Ca = μcUHR/σ ≈ 1.09 × 10⁻⁵. The density-contrast Bond number BoD = |ρd − ρc|gD²/σ = 1.635 × 10⁻⁴ compares buoyancy with capillarity. These small ratios support the chosen low-inertia, weak-deformation approximation; they are not a universal error bound.

Surfactants can modify surface mobility through their distribution and tangential stresses, so an experimentally prepared drop need not follow the clean result. Walls, neighboring drops, deformation and inertial wakes can change drag too. The rigid and clean values here are two specified boundary-condition models, not universal bounds for every driven, contaminated or deforming inclusion. No transient settling time has been calculated.

Heat supplies latent enthalpy; mass need not follow the front

Now use a separate local planar evaporation example. The normal n points from liquid to vapor. Assign heat flux qL·n = 12 000 W/m² arriving from the liquid side and qV·n = 2000 W/m² leaving into the vapor side. The difference is energy supplied to phase change. Take latent enthalpy hV − hL = 2.0 × 10⁶ J/kg, liquid density 1000 kg/m³ and vapor density 2 kg/m³.

With negligible kinetic, viscous and interfacial-storage energy terms, the local heat balance is m″(hV − hL) = qL·n − qV·n. Thus m″ = 10 000/(2.0 × 10⁶) = 0.005 kg/(m²·s), positive evaporation. Across 0.010 m² of interface the net heat is 100 W and the converted mass rate is 0.000050 kg/s, or 0.050 g/s. These are supplied heat fluxes; no temperature field or saturation pressure has been solved.

Choose a stationary interface, VΓ = 0. The mass jump then gives liquid normal speed 5 × 10⁻⁶ m/s and vapor normal speed 0.0025 m/s, both along n. Different densities require different speeds for the same crossing mass. The normal velocity jump is 0.002495 m/s. A stationary boundary can therefore have mass passing through it, unlike a material boundary.

A different laboratory motion can have liquid at rest. For the same relative mass flux, its front recedes into the liquid at VΓ = −5 × 10⁻⁶ m/s and its vapor speed is 0.002495 m/s. This is an alternative frame of the local kinematics, not an additional velocity in the stationary-front case. In either description the vapor moves relative to the front at 2.5 mm/s.

The front-relative kinetic-enthalpy difference is about 3.125 × 10⁻⁶ J/kg, tiny compared with the assigned latent enthalpy. That checks one omitted energy term in this example. A full solution still needs compatible temperature, pressure, stress and interface conditions; the heat and mass balances alone do not determine them. Condensation has m″ < 0 with the same convention, so reversing the net heat transfer reverses the sign rather than changing the definition.

Bubbles can have different origins

An injected air bubble is not automatically vapor made from the surrounding liquid. A bubble may contain vapor, noncondensable gas or both. Its response can involve gas compression, heat transfer, species transfer and phase change in different proportions. Naming it “gas” does not select an equation of state or an evaporation rate.

Boiling commonly follows heating relative to the local saturation condition; cavitation commonly follows a pressure reduction relative to the liquid's vapor-pressure condition. Both can form vapor, but their thermal and pressure histories differ. Nucleation sites, impurities, dissolved gas, boundaries and the time available affect inception. A saturation curve is an equilibrium guide, not a universal rule that every initially smooth liquid immediately produces a specified bubble at its crossing.

Condensation transfers vapor into liquid and releases latent enthalpy. A liquid film and a collection of drops offer different heat-transfer paths and drainage. In boiling equipment a change in wall wetting or the formation of an insulating vapor region can alter the heat removed for a given wall condition. Those device outcomes need geometry and heat-transfer models beyond our assigned planar flux balance.

Further: phase balances reveal the missing closures

For a representative phase-average description, write ∂t(αₖρₖ) + ∇·(αₖρₖuₖ) = Γₖ. Here ρₖ is the intrinsic mean density of phase k, and uₖ is its density-weighted intrinsic mean velocity: the phase mass flux divided by αₖρₖ where that phase is present. If density is uniform within the phase, this equals its ordinary intrinsic volume-average velocity. With variable density, using the volume-average velocity in this mass balance would require a density–velocity correlation term. Γₖ is mass gained by phase k per mixture volume per time. Without interfacial mass storage, ΓG + ΓL = 0: the total mass balance contains no source from merely changing phase. The individual phase balances still have equal and opposite sources.

Separate phase momentum equations include transport, phase stresses, body forces and interphase momentum exchange. Drag is one contribution; acceleration of surrounding fluid, shear-induced lateral forces, fluctuating motion and mass crossing the interface can matter in other regimes. Consistent internal exchanges cancel in the combined whole-system balance when interfacial forces, storage and surface stresses are treated consistently. Their cancellation in a total balance does not determine either phase's motion.

A two-fluid model retains separate phase velocities and usually separate thermal variables. It needs constitutive or closure relations for quantities not determined by the averaged conservation equations: interfacial drag and area, turbulent correlations, heat transfer and phase-change rates, for example. A mixture model reduces some of these unknowns by prescribing relative motion or relaxation. A homogeneous common-velocity approximation suppresses slip; a homogeneous equilibrium model adds thermodynamic-equilibrium assumptions. Common velocity alone does not imply common temperature or instantaneous phase equilibrium.

The averaging scale must be stated. A cell that resolves one interface and a large cell containing many unresolved bubbles can both report a fraction between 0 and 1, but that number has a different role in each model. A successful average must retain the information or closures needed by its intended question, and should be checked against measurements or a more resolved calculation in that regime.

Further: a drift identity is not yet a closure

In the assigned uniform cell, jG = αGj + αGαL(uG − uL) is an exact rearrangement. It becomes predictive only when the slip is supplied or determined by additional dynamics. Calling the last term drift flux does not tell us what it should be in a different pipe.

Across a nonuniform section, let ⟨·⟩ denote an area average of locally defined two-phase fields. Since jG = αGuG pointwise, ⟨jG⟩ = ⟨αGj⟩ + ⟨αG(uG − j)⟩. For nonzero ⟨αG⟩ and ⟨j⟩, define C0 = ⟨αGj⟩/(⟨αG⟩⟨j⟩) and Vgj = ⟨αG(uG − j)⟩/⟨αG⟩. Then ⟨jG⟩ = ⟨αG⟩[C0⟨j⟩ + Vgj]. These are definitions until profiles or a calibrated model supply C0 and Vgj.

Our uniform cell has C0 = 1 and Vgj = 0.800 m/s, recovering jG = 0.300 m/s. Real section profiles need not give C0 = 1. A fitted relation for one flow arrangement does not automatically apply to a stratified flow, a different orientation or a phase-changing condition. Countercurrent or zero-total-volume-flux cases also demand care because a normalized expression can become singular even when physical phase fluxes remain finite.

Further: a size population has its own transport

When size controls drag, residence time or transfer area, one average diameter may discard essential information. Let n(D, x, t)dD be the number of spherical inclusions per mixture volume with diameters between D and D + dD; n has units m⁻⁴ when D is in metres. A schematic population balance is ∂tn + ∇·(un) + ∂D(Gn) = births − deaths, where u is the size-conditioned transport velocity and G = dD/dt is continuous growth or shrinkage in m/s. Nucleation, breakup and coalescence contribute to the birth and death terms.

The diameter derivative is transport through size space: an evaporating drop can leave one diameter interval and enter a smaller one without being newly born. Breakup and coalescence redistribute the population differently. For constant-density spheres without phase transfer, admissible breakup/coalescence kernels conserve total dispersed volume even though they change number and interfacial area. Phase change may change the dispersed volume through the growth term and the coupled mass balance.

This framework states what must be counted; it does not supply universal collision, breakup or coalescence rates. Surface contamination, turbulence, film drainage and local concentration can influence those rates. The spherical area identity and the general balance remain useful organizing tools while the actual closures and their applicable regime remain explicit.

Choose the resolution that serves the question

An interface-resolved calculation follows the geometry of the boundary and solves adjacent fluid motion. Volume-of-fluid methods transport a cell's occupied fraction; level-set methods represent a boundary through a level surface of a field; front-tracking methods move an explicit interface representation. Each requires a numerical treatment of geometry, transport and interfacial forces. Naming a method does not establish that its grid resolves a thin film, a small bubble or a breakup event.

An Eulerian–Eulerian average treats phases as interpenetrating averaged fields over the same region and uses closures for unresolved interactions. An Eulerian–Lagrangian model follows individual inclusions or representative parcels while the carrier is a continuum; a point-inclusion treatment still needs an appropriate force and transfer law. A parcel may represent many physical drops. Resolving their tracks does not necessarily resolve each surrounding flow or interface.

Useful checks include phase fractions remaining bounded, total and phase mass consistency, coupled energy conservation with latent enthalpy, appropriate force balance, and refinement of predicted quantities. Artificial currents near a supposedly static interface or grid-dependent loss of a film are warning signs, not physical mixing evidence. Agreement in one terminal-drop test cannot establish correctness for a boiling channel.

Further: constant density can give ∇·u = 0 inside each bulk phase while phase change creates a jump at the interface. In a sharp-interface one-field description, ∇·u = m″(1/ρV − 1/ρL)δΓ in the distributional sense, where δΓ concentrates the source on the interface. Integrating across the interface recovers the normal velocity jump. Imposing one globally zero-divergence velocity without accommodating this expansion would contradict the supplied unequal-density evaporation kinematics.

Ask what would change the conclusion

In a separator, relative motion helps phases leave through different outlets; droplet size and interfacial mobility can change residence requirements. In a gas–liquid contactor, area and residence time affect transfer while coalescence changes the area inventory. In a spray, size-conditioned trajectories and evaporation compete. In a heat-removal device, phase conversion carries latent enthalpy but also changes volume flux and interface arrangement.

The practical question determines which effects to retain. A density average may suffice to calculate stored mass but cannot select the drag on a contaminated drop. A clean isolated-drop law may explain a mobility effect but cannot select a crowded-flow regime. The local evaporation balance converts an assigned heat supply to mass flux, but cannot predict boiling inception or safe device capacity without the missing thermal and mechanical problem.

For a new case, identify the phases and topology, choose a sampling scale, distinguish amounts from fluxes, list allowed interfacial transfers, and check the model's dimensionless ratios and boundary conditions. Then ask whether the required closures have evidence for that arrangement. The next planned subject, “Microfluidic, reacting, biological, and geophysical flows”, connects these tools to additional small-scale, chemical, biological and rotating or stratified settings; those extra mechanisms are not silently added to the models here.

Common shortcuts worth testing

“Gas occupies 20%, so it carries 20% of the mass.” First choose inventory or flow, then include density and, for flow, velocity. “The mixture has one velocity.” State which average is intended and whether slip has been discarded or separately modeled. “Density averaging gives every mixture property.” Stored mass is additive; effective viscosity and transfer laws require more information.

“Smaller bubbles always improve transfer by the same factor.” Smaller equal spheres at fixed amount expose more area, but transfer per area and residence can change. “A clean drop obeys rigid Stokes drag.” The interface conditions differ, as the two terminal calculations show. “An incompressible liquid–vapor model has one continuous divergence-free velocity.” Phase-change mass conservation with unequal densities gives a normal jump even when both bulk densities stay constant.

“A plausible bubble picture validates the simulation.” Conservation, force and energy checks, refinement, resolved scales and suitable experimental comparisons provide the needed evidence. Keep the physical question, the assigned inputs and the model's unanswered questions visible alongside the result.

Worked examples

Try it yourself

1. In the supplied moving cell, is the gas volume-flow fraction 20%, 42.86% or 0.1498%? Explain which measurement each number describes.

   Show the answer and reasoning

   The gas volume-flow fraction is about 42.86%.

   Gas volume flow is 0.3 L/s of 0.7 L/s total. The 20% is occupied volume fraction; about 0.1498% is the positive mass-flow quality after density weighting.

2. For the supplied cell, Π = 200.900 Pa and ρmum² ≈ 200.500200 Pa. Find the remaining convective transport and decide whether it is a calculated thermodynamic pressure.

   Show the answer and reasoning

   D ≈ 0.399800 Pa; it is relative-motion momentum transport, not a calculated pressure.

   Subtract the mean term from the sum of phase convective fluxes. Its pressure unit follows momentum transported per area per time. Pressure and viscous stress were excluded, so dimensional agreement does not identify the physical quantity.

3. With j = 0.700 m/s and uG = 1.500 m/s, find the gas drift velocity and drift flux at αG = 0.2.

   Show the answer and reasoning

   Drift velocity is 0.800 m/s; drift flux is 0.160 m/s.

   The velocity is uG − j. The flux is αG times that relative velocity. Both have m/s units, but only the flux includes the amount of gas available to cross.

4. At fixed total dispersed volume, equal spheres double in diameter. By what factors do representative count and total interface change?

   Show the answer and reasoning

   Count becomes one eighth and total interface becomes one half.

   One sphere's volume scales as D³ and its area as D². Holding the total volume fixed requires count proportional to D⁻³, so total area scales as D⁻¹. This is a geometric inventory, not a coalescence-time prediction.

5. A small drop has an unknown surfactant distribution. Is the clean λ = 1 speed 0.3270 mm/s automatically its terminal speed?

   Show the answer and reasoning

   No; the interfacial condition is not established.

   The clean result assumes uniform surface tension and the corresponding mobile fluid interface. Surfactant gradients can change tangential stresses and circulation. Measurements or a justified surface-transport model are needed; an unspecified interface is not proven clean by its size.

6. Keep the same densities and latent enthalpy but let the signed net heat qL·n − qV·n be −10 000 W/m². What are m″ and the stationary-front normal velocities?

   Show the answer and reasoning

   m″ = −0.005 kg/(m²·s), uL·n = −0.005 mm/s and uV·n = −2.500 mm/s.

   The unchanged positive direction is liquid to vapor. A negative mass flux denotes condensation, with both normal velocities toward the liquid in the stationary-front description. Changing the process does not require silently reversing the normal.

7. For the original positive evaporation flux, choose the laboratory liquid velocity to be zero. Find the interface and vapor speeds along n.

   Show the answer and reasoning

   VΓ = −0.005 mm/s and uV·n = 2.495 mm/s.

   From m″ = ρL(0 − VΓ), the interface recedes into liquid. The vapor must move at m″/ρV = 2.500 mm/s relative to that moving interface, giving 2.495 mm/s in this laboratory frame. These are alternative-frame kinematics, not extra velocities in the stationary-front case.

8. A simulation reports the correct total mass and a bounded mean gas fraction but uses no size or interfacial-area information. Has it thereby established an accurate gas–liquid heat-transfer rate?

   Show the answer and reasoning

   No; the relevant interface and thermal closure remain unverified.

   Conservation and bounded fractions are necessary checks. Different size distributions can have the same amount and different area, while temperature gradients and transfer coefficients can differ too. The model needs appropriate resolved information or validated closures for the heat-transfer question.

What to carry forward

• Occupied fraction, stored mass fraction and crossing-flow fractions answer different questions; density and relative phase speeds can make them very different.
• Slip survives averaging as additional transport. A mean, an exact identity and a predictive closure are distinct steps.
• Phase arrangement and interface area matter: equal spherical inclusions at fixed amount have area proportional to inverse diameter.
• Clean drops, rigid particles and contaminated interfaces impose different mechanical conditions; use drag laws with their stated regime.
• Phase-change mass crosses relative to the moving front. Unequal phase densities can require a velocity jump even with constant bulk densities.
• Choose resolved or averaged models for the question, retain the required closures, and verify conservation, conditions, resolution and physical comparisons.