Model comminution kinetics and specific surface area

verifiedAgainst e41a048 · verifiedOn 2026-09-11.

What you will build

In traditional sandbox crafting, processing machinery (like a crusher, pulverizer, or ball mill) exists as an artificial gate: a recipe takes 10 seconds in a tier-1 machine and 5 seconds in a tier-2 machine.

In Voxamine, crushing exists because chemical physics requires reactive surface area (MASTERPLAN.md §28.5, §32). For heterogeneous solid–gas (e.g. roasting chalcopyrite in oxygen) or solid–liquid reactions (e.g. acid leaching of malachite), reaction occurs strictly at the exposed solid interface:

$$r = k(T) \cdot \prod [reactant]^\nu \cdot S_{\text{area}} \cdot \left(1 - \frac{Q}{K}\right)$$

A $1\text{ cm}$ chunk of dense ore has negligible surface area per mole of mineral. Milling that ore down to $100\ \mu\text{m}$ powder increases the exposed reactive interfacial area by a factor of 100, accelerating the reaction rate by two full orders of magnitude. The crusher is not a balance knob; it is what the chemical rate law directly demands.

You will compute solid molar volume $V_m = \frac{M / 1000}{\rho}$ and reactive molar interfacial area $A_m = \frac{6 V_m}{d}$ in SpecificSurfaceAreaMath.cs, plug the heterogeneous surface term into ReactionKineticsMath.cs, evaluate temperature-dependent Arrhenius rate constants $k(T) = A \exp(-E_a / RT)$ and thermodynamic driving force $(1 - Q/K)$, and verify zero-allocation performance on solver hot paths.

// 1. Evaluate reactive surface area for crushed copper ore (d = 250 µm):
SpecificSurfaceAreaMath.TryCalculateHeterogeneousSurfaceTerm(
    grainSizeMeters: 2.5e-4d,
    molarMassGramsPerMol: 63.546d,
    densityKilogramsPerCubicMeter: 8960d,
    out double surfaceTermMolarArea);

// 2. Evaluate net reaction rate in the vessel solver:
ReactionKineticsMath.TryCalculateNetRate(
    rateConstant: k,
    activityOrderProduct: activityProd,
    specificSurfaceArea: surfaceTermMolarArea,
    reactionQuotient: Q,
    equilibriumConstant: K,
    out double netRate);

Where this sits

Comminution kinetics connects raw physical mass to the chemical differential equation solver:

• Particle & Charge Profiles (ItemDefinition / CarriedBulkMatter) defines solid grain size $d$, density $\rho$, and molar mass $M$.
• Surface Area Math (SpecificSurfaceAreaMath) derives molar volume $V_m$ and reactive interface $A_m = \frac{6 V_m}{d}$ for uniform spherical geometries.
• Kinetic Identities (ReactionKineticsMath) combines the Arrhenius rate constant, activity products, and distance from equilibrium $(1 - Q/K)$ into a signed net rate $r$.
• Vessel Solver (VesselSolver & DeterministicSubstepping) advances sub-tick extents $\xi$ and updates reactant masses inside metallurgical furnaces and leaching tanks.

flowchart TD
  subgraph SOLID["Solid Charge Properties"]
    GRAIN["Particle Grain Size d (m)\n(Lump 10mm → Sand 1mm → Pulp 75µm)"]
    PROPS["Molar Mass M (g/mol)\nSolid Density ρ (kg/m³)"]
  end

  subgraph AREA["SpecificSurfaceAreaMath"]
    VOL["Molar Volume\nVm = (M / 1000) / ρ (m³/mol)"]
    MOLAR["Molar Surface Area\nAm = 6 · Vm / d (m²/mol)"]
    HOMO["Homogeneous Reaction\nS_area = 1.0 (dimensionless)"]
  end

  subgraph KINETICS["ReactionKineticsMath"]
    ARRHENIUS["k(T) = A · exp(-Ea / RT)\n(Arrhenius Rate Constant)"]
    EQUILIB["Equilibrium Factor\n(1 - Q / K)"]
    NET["Net Rate r = k · Π[a]^ν · S_area · (1 - Q/K)"]
    AFFINITY["Affinity A = RT · (ln K - ln Q)"]
  end

  subgraph SOLVER["Chemistry Vessel Solver"]
    TICK["DeterministicSubstepping\n(Advances reaction extent dξ = r · dt)"]
    MASS["Phase Transition & Vessel Mass Balance"]
  end

  PROPS --> VOL
  VOL --> MOLAR
  GRAIN --> MOLAR
  MOLAR --> NET
  HOMO --> NET
  ARRHENIUS --> NET
  EQUILIB --> NET
  NET --> TICK
  AFFINITY --> TICK
  TICK --> MASS

Before you start

Read Docs/CHEMISTRY_SOURCES.md and MASTERPLAN.md §28.5 (“Heterogeneous reactions & rate laws”), then inspect:

• SpecificSurfaceAreaMath.cs, geometry and surface area derivations.
• ReactionKineticsMath.cs, deterministic rate and affinity equations.
• SpecificSurfaceAreaMathTests.cs, unit tests asserting inverse diameter scaling and zero-allocation performance.

Keep the distinction between reaction mechanisms clear:

Regime	Example	Surface Term $S_{\text{area}}$	Units	Mechanism
Homogeneous	$2,\text{CO} + \text{O}_2 \to 2,\text{CO}_2$ (gas phase)	$1.0$ (HomogeneousSurfaceTerm)	Dimensionless	Molecules collide freely throughout the phase volume.
Heterogeneous	$\text{Cu}2\text{O}{(s)} + \text{CO}{(g)} \to 2,\text{Cu}{(s)} + \text{CO}_{2(g)}$	$\frac{6 V_m}{d}$	$\text{m}^2/\text{mol}$	Gas molecules must strike the physical solid crystal surface.

The build, step by step

1. Derive Solid Molar Volume

In SpecificSurfaceAreaMath.TryCalculateMolarVolumeCubicMetersPerMol, calculate the molar volume $V_m$ ($\text{m}^3/\text{mol}$) from molar mass $M$ ($\text{g/mol}$) and solid crystal density $\rho$ ($\text{kg/m}^3$):

public static bool TryCalculateMolarVolumeCubicMetersPerMol(
    double molarMassGramsPerMol,
    double densityKilogramsPerCubicMeter,
    out double molarVolumeCubicMetersPerMol)
{
    if (!IsFinitePositive(molarMassGramsPerMol) || !IsFinitePositive(densityKilogramsPerCubicMeter))
    {
        molarVolumeCubicMetersPerMol = 0d;
        return false;
    }

    molarVolumeCubicMetersPerMol = molarMassGramsPerMol / GramsPerKilogram
        / densityKilogramsPerCubicMeter;
    return IsFinitePositive(molarVolumeCubicMetersPerMol);
}

For elemental copper ($M = 63.546\ \text{g/mol}$, $\rho = 8960\ \text{kg/m}^3$): $$V_m = \frac{0.063546\ \text{kg/mol}}{8960\ \text{kg/m}^3} \approx 7.0922 \times 10^{-6}\ \text{m}^3/\text{mol}$$

2. Derive Reactive Molar Interfacial Area

For uniform spherical particles of diameter $d$:

• Particle volume: $v_p = \frac{\pi}{6} d^3$
• Particle surface area: $a_p = \pi d^2$
• Area-to-volume ratio: $\frac{a_p}{v_p} = \frac{\pi d^2}{\frac{\pi}{6} d^3} = \frac{6}{d}$

In SpecificSurfaceAreaMath.TryCalculateMolarSurfaceAreaSquareMetersPerMol:

public static bool TryCalculateMolarSurfaceAreaSquareMetersPerMol(
    double grainSizeMeters,
    double molarVolumeCubicMetersPerMol,
    out double molarSurfaceAreaSquareMetersPerMol)
{
    if (!IsFinitePositive(grainSizeMeters) || !IsFinitePositive(molarVolumeCubicMetersPerMol))
    {
        molarSurfaceAreaSquareMetersPerMol = 0d;
        return false;
    }

    molarSurfaceAreaSquareMetersPerMol = 6d * molarVolumeCubicMetersPerMol / grainSizeMeters;
    return IsFinitePositive(molarSurfaceAreaSquareMetersPerMol);
}

Similarly, mass-specific surface area $a_s$ ($\text{m}^2/\text{kg}$) equals: $$a_s = \frac{6}{\rho \cdot d}$$ Notice that $A_m = a_s \cdot \frac{M}{1000}$, providing a rigorous unit consistency bridge between bulk mass and stoichiometric moles.

3. Single-Step Heterogeneous Surface Term

Combine both calculations into a single zero-allocation call in SpecificSurfaceAreaMath.TryCalculateHeterogeneousSurfaceTerm:

public static bool TryCalculateHeterogeneousSurfaceTerm(
    double grainSizeMeters,
    double molarMassGramsPerMol,
    double densityKilogramsPerCubicMeter,
    out double surfaceTermSquareMetersPerMol)
{
    if (!TryCalculateMolarVolumeCubicMetersPerMol(
            molarMassGramsPerMol, densityKilogramsPerCubicMeter, out double molarVolume))
    {
        surfaceTermSquareMetersPerMol = 0d;
        return false;
    }

    return TryCalculateMolarSurfaceAreaSquareMetersPerMol(
        grainSizeMeters, molarVolume, out surfaceTermSquareMetersPerMol);
}

4. Evaluate Arrhenius Rate Constant $k(T)$

In ReactionKineticsMath.TryCalculateArrheniusRateConstant, calculate the temperature-dependent forward rate constant:

$$k(T) = A \exp\left(-\frac{E_a \cdot 1000}{R \cdot T}\right)$$

double exponent = -activationEnergyKilojoulesPerMol * ChemicalConstants.JoulesPerKilojoule
    / (ChemicalConstants.MolarGasConstantJoulesPerMolKelvin * temperatureKelvin);
rateConstant = preExponentialFactor * DeterministicMath.Exp(exponent);

Using deterministic fixed-point or bounded exponential math guarantees platform-identical simulation results on Windows, Linux, and macOS.

5. Calculate Net Rate with Equilibrium Reversal

In ReactionKineticsMath.TryCalculateNetRate, compute the signed reaction rate:

public static bool TryCalculateNetRate(
    double rateConstant,
    double activityOrderProduct,
    double specificSurfaceArea,
    double reactionQuotient,
    double equilibriumConstant,
    out double netRate)
{
    if (!IsFiniteNonNegative(rateConstant)
        || !IsFiniteNonNegative(activityOrderProduct)
        || !IsFiniteNonNegative(specificSurfaceArea)
        || !IsFinitePositive(reactionQuotient)
        || !IsFinitePositive(equilibriumConstant))
    {
        netRate = 0d;
        return false;
    }

    double equilibriumFactor = 1d - reactionQuotient / equilibriumConstant;
    netRate = rateConstant * activityOrderProduct * specificSurfaceArea * equilibriumFactor;
    return IsFinite(netRate);
}

Key invariant: The net rate changes sign naturally when $Q > K$:

• When $Q < K$: $(1 - Q/K) > 0 \implies r > 0$ (forward reaction proceeds).
• When $Q = K$: $(1 - Q/K) = 0 \implies r = 0$ (chemical equilibrium).
• When $Q > K$: $(1 - Q/K) < 0 \implies r < 0$ (reverse reaction occurs spontaneously).

No separate reverse reaction definition is ever required.

6. Calculate Reaction Affinity without Ratio Overflow

In ReactionKineticsMath.TryCalculateReactionAffinityKilojoulesPerMol, calculate chemical affinity $A = -\Delta_r G = R T \ln(K / Q)$:

// Subtract logs instead of forming K/Q, avoiding an otherwise needless
// overflow when a strongly favourable reaction is far from equilibrium:
affinityKilojoulesPerMol = ChemicalConstants.MolarGasConstantJoulesPerMolKelvin
    * temperatureKelvin
    * (DeterministicMath.Ln(equilibriumConstant) - DeterministicMath.Ln(reactionQuotient))
    / ChemicalConstants.JoulesPerKilojoule;

Subtracting logarithms ($\ln K - \ln Q$) prevents double overflow even when $K = 10^{300}$ and $Q = 10^{-300}$.

Screenshot pending

Capture crushed-ore-surface-area-curve.png: Diagnostic plot showing inverse grain-size scaling $d$ vs $S_{\text{area}}$ and comparative reaction rate curves for lump ore (10 mm) vs ball-mill pulp (75 µm). Editor dark theme.

Verify

1. Execute EditMode Tests Headlessly

Run the dedicated kinetics test suites:

/opt/unity/Editor/Unity -batchmode -nographics \
  -projectPath /home/soulwax/workspace/engines/unity/minecraft/Minecraft-HD \
  -runTests -testPlatform EditMode \
  -testFilter VoxelSandbox.Tests.SpecificSurfaceAreaMathTests

/opt/unity/Editor/Unity -batchmode -nographics \
  -projectPath /home/soulwax/workspace/engines/unity/minecraft/Minecraft-HD \
  -runTests -testPlatform EditMode \
  -testFilter VoxelSandbox.Tests.ReactionKineticsMathTests

Verify that all key assertions pass:

• MolarSurfaceArea_ScalesInverselyWithGrainSize: asserts that a $10\times$ decrease in grain diameter yields exactly a $10\times$ increase in reactive surface area.
• MassSpecificAndMolarSurfaceArea_AgreeThroughMolarMass: verifies $A_m = a_s \cdot (M / 1000)$ across double precision bounds.
• NetRate_ChangesSignAtEquilibriumWithoutAnAuthoredReverseReaction: proves forward ($+18$), equilibrium ($0$), and reverse ($-24$) kinetics.
• Affinity_HasTheExpectedDirectionAndAvoidsRatioOverflow: verifies stability at extreme quotients like $Q = 10^{-300}$ and $K = 10^{300}$.
• HotPath_DoesNotAllocateManagedMemory: verifies zero GC allocation over $10{,}000$ iterations using GC.GetAllocatedBytesForCurrentThread().

Now do your own

To incorporate comminution into a new ore reduction or hydrometallurgy process (e.g. Fine Bauxite Digestion or Pyrite Pressure Oxidation):

1. Define Grain Size Stages: Set realistic grain sizes for processing stages:
   • Run-of-mine lump: $d = 20\text{ mm}$ ($2 \times 10^{-2}\text{ m}$)
   • Jaw crushed gravel: $d = 2\text{ mm}$ ($2 \times 10^{-3}\text{ m}$)
   • Ball-mill pulp: $d = 75\ \mu\text{m}$ ($7.5 \times 10^{-5}\text{ m}$)
2. Calculate $S_{\text{area}}$: Evaluate TryCalculateHeterogeneousSurfaceTerm using mineral crystal density and formula molar mass.
3. Plug into Vessel Solver: Pass $S_{\text{area}}$ to ReactionKineticsMath.TryCalculateNetRate.
4. Observe Roasting Time: Watch run-of-mine lumps take hours of game time while fine pulp roasts in seconds!

Pitfalls

Symptom	Cause	Fix
Crushing an ore increases the final equilibrium yield of metal.	Mistook reaction kinetics ($r$) for chemical equilibrium ($K$).	Comminution only increases reaction speed $r$, never the equilibrium constant $K$. Final equilibrium yield depends solely on thermodynamics ($\Delta_r G^\circ$) and stoichiometry.
Reaction rate is identical regardless of whether ore is crushed or whole.	Reaction was authored with HomogeneousSurfaceTerm ($1.0$) instead of calling TryCalculateHeterogeneousSurfaceTerm.	Use the heterogeneous surface calculation for solid–gas and solid–liquid reactions.
Solver throws NaN or Infinity during vigorous reactions.	Direct division equilibriumConstant / reactionQuotient overflowed double range.	Use DeterministicMath.Ln(K) - DeterministicMath.Ln(Q) when computing chemical affinity.
GC allocation spikes during vessel reaction ticks.	Created temporary arrays or class instances inside kinetic evaluation loops.	Keep all math functions in static classes with out double parameters, asserting zero allocations in tests.
Molar surface area does not match mass-specific area.	Forgot to convert molar mass from grams to kilograms ($M / 1000$).	Ensure density is in $\text{kg/m}^3$ and molar mass is divided by $1000$ to obtain consistent SI units ($\text{m}^3/\text{mol}$).