High: [Drying] Gas balance and Darcy velocity produce order-of-magnitude drying transient errors

[sergey-gusev94]

Maintainer triage (June 29, 2026)

• Status: Verified bug against the current master source in this fork.
• Severity: High.
• Area labels: area:drying.
• Tackle batch: P2 - shared thermo, phase, and basis correctness.

The source review confirmed that the cited code path and failure mode are still present. The original audit details are preserved below for reproduction notes and suggested fixes.

---

Original audit details

Severity: High · Part of #67 (round-2 audit batch)

All three defects live in PharmaPy/Drying_Model.py in the same drying subsystem and are exercised on the same code paths (material_balance RHS + solve_unit initialization). The gas-velocity (243-244) and gas-species-source (334) defects both corrupt the gas-phase balance and would be fixed and validated together in a single drying-correctness pass; the single-node init omission (472) is the same solve_unit setup that feeds that balance. Bundling avoids three overlapping PRs touching the same file and lets one regression test cover the coupled transient. The Wilke vapor-viscosity defect is intentionally kept separate (different file, general thermo property).

This issue bundles 3 defects that share one fix/PR.

1. (High) Darcy gas velocity discards relative permeability k_ra and divides by mean saturation (wrong gas flow at both saturation limits)

Location: PharmaPy/Drying_Model.py:243-244 (k_ra computed but unused at 239-241; correct form commented at 242)

Lines 239-241 compute the two-phase gas relative permeability k_ra = (1-sat_red)^2*(1-sat_red^1.4) but it is never used. The active line 243-244 computes vel_gas = self.dPg_dz/(alphavisc_gasrho_sol*(1-porosity))/np.mean(satur), i.e. single-phase Darcy velocity (= k_permdPg_dz/visc, since k_perm=1/alpha/rho_sol/(1-porosity) at line 520) divided by mean saturation. This (a) drops the k_ra throttle entirely and (b) substitutes a 1/mean(S) factor that runs backwards at the limits: at full saturation (S~1) gas cannot flow but the code returns full single-phase velocity, and as S decreases the velocity grows instead of approaching the finite k_ra->1 limit. The physically correct, intended form (k_permk_ra*dPg_dz/visc_gas) is the commented line 242.

• Trigger: Any normal Drying.solve_unit run: unit_model is the integrator RHS, so line 243-244 executes on every CVode step; vel_gas (u_gas) feeds gas convection (line 332), the mass-transfer denominator (line 334), and the energy-convection term (line 387).
• Wrong result: Gas superficial velocity is wrong in magnitude and trend throughout the transient (full single-phase speed even when the cake is liquid-saturated; growing as the cake dries instead of bounded by k_ra). This corrupts gas convection, mass-transfer, and energy-convection terms, distorting the entire moisture and temperature profile. Note: amplification is bounded (mean(satur) floors near the irreducible saturation ~0.155, giving ~6.5x, not the 1/eps blow-up originally claimed), but the physics is qualitatively wrong.
• Suggested fix: Restore the relative-permeability form by uncommenting line 242: vel_gas = self.k_perm * k_ra * self.dPg_dz / visc_gas, and remove the /np.mean(satur) factor on line 243-244. Ensure gas velocity -> 0 as S -> 1 and stays finite as S -> 0.
• Net-new (not a duplicate): No filed issue mentions k_ra, vel_gas, Darcy velocity, relative permeability, or np.mean(satur). Distinct from Critical: [Drying] Liquid saturation balance mixes molar drying rate with mass density #20/Critical: [Drying] Gas species balance adds molar source terms to mass-fraction states #21/High: [Drying] Energy terms multiply molar drying rates by mass-basis heat quantities #48 (molar-vs-mass dry_rate), High: [Drying] Condensed-phase energy balance doubles latent heat #24 (latent-heat factor), High: [Drying] Gas molecular weight is computed from mass fractions as mole fractions #28 (mean-MW), Medium: [Drying] Gas energy balance multiplies convective enthalpy by density twice #37 (extra rho_gas in convection), Medium: [Drying] Volatile liquid mass-fraction column is forced to zero each RHS evaluation #42 (x_liq hack).

2. (High) Gas-phase species evaporation source over-divided by an extra (1-satur) factor

Location: PharmaPy/Drying_Model.py:334

epsilon_gas (line 324) is already porosity*(1-satur), the gas volume fraction. The correct accumulation coefficient for a mass-fraction source is eps_grho_g, so the source should be N_i/(porosity(1-satur)rho_gas). But line 334 codes transfer_gas = dry_rate.T/epsilon_gas/dens_gas/(1-satur), dividing additionally by (1-satur), i.e. by porosity(1-satur)^2rho_gas. The companion total_mass_correction (line 336) already correctly carries the single 1/(1-satur) from d(eps_g)/dt, and the energy balance (line 362, denom_gas = cvg_mixepsilon_gas*rho_gas) uses exactly one (1-satur) for the same accumulation coefficient, directly contradicting line 334.

• Trigger: Any normal Drying.solve_unit run. Saturation starts near 1 for a fresh wet cake, so (1-satur) is small (~0.05-0.2), making the spurious factor inflate the evaporative source ~5x-20x in the early high-saturation regime that dominates the transient.
• Wrong result: The gas-phase y_i source (transfer_gas) is too large by 1/(1-satur), so the predicted gas-composition rise from evaporation is grossly overstated and the gas-phase species balance does not match the liquid loss; the whole coupled drying transient is biased.
• Suggested fix: Remove the trailing /(1 - satur) on line 334: transfer_gas = dry_rate.T / epsilon_gas / dens_gas. (Independently apply the MW conversion of Critical: [Drying] Gas species balance adds molar source terms to mass-fraction states #21.) The dynamic-saturation correction is already handled separately by total_mass_correction on line 336.
• Net-new (not a duplicate): Distinct from Critical: [Drying] Gas species balance adds molar source terms to mass-fraction states #21 (which is the mol-vs-mass MW error on the same line 334). Both must be corrected; fixing Critical: [Drying] Gas species balance adds molar source terms to mass-fraction states #21 does not remove this extra (1-satur).

3. (High) Uniform/single-node initialization omits temp_cond from the initial state vector

Location: PharmaPy/Drying_Model.py:470-475 (else branch; cf. correct if branch at 463-468)

When x_liq_init.ndim==1 and len(satur_init)==1, line 472 builds states_tuple = (satur_init, y_gas_init, x_liq_init, temp_gas_init), omitting temp_cond_init, giving a per-node width of 2+num_comp+num_volatiles. The parallel branch (lines 463-468) correctly stacks all five state groups (S, y_gas, x_liq, Tg, Tcond). unit_model later does states.reshape(-1, 3 + num_comp + self.num_volatiles) (line 225), which expects the temp_cond column.

• Trigger: Drying.solve_unit with a Cake whose saturation is a single-element array (single representative node) and a uniform (1-D mass_frac) liquid composition, which hits the else branch at line 470.
• Wrong result: model(0, states_prev.ravel()) at line 575 raises 'cannot reshape array of size N into shape (...)' (crash before the integrator starts) in ~87% of realistic dims; in the ~13% where divisibility coincidentally holds, state columns are silently misaligned (temp parsed from x_liq), corrupting every balance.
• Suggested fix: Include temp_cond_init in the tuple to match the other branch: states_tuple = (satur_init, y_gas_init, x_liq_init, temp_gas_init, temp_cond_init).
• Net-new (not a duplicate): No filed issue concerns the Drying y0/state-assembly in solve_unit; known Drying issues (Critical: [Drying] Liquid saturation balance mixes molar drying rate with mass density #20,Critical: [Drying] Gas species balance adds molar source terms to mass-fraction states #21,High: [Drying] Condensed-phase energy balance doubles latent heat #24,High: [Drying] Gas molecular weight is computed from mass fractions as mole fractions #28,Medium: [Drying] Gas energy balance multiplies convective enthalpy by density twice #37,Medium: [Drying] Volatile liquid mass-fraction column is forced to zero each RHS evaluation #42,High: [Drying] Energy terms multiply molar drying rates by mass-basis heat quantities #48) are all in the balance/property code, not solve_unit init.

Related existing issues (referenced, not modified): #20, #21, #24, #28, #37, #42, #48 (existing drying balance/property issues in the same module; reference only). #21 specifically shares line 334 (mol-vs-mass) but is a distinct mechanism from the extra (1-satur) factor here.

---

Found by a multi-agent source audit (round 2); confirmed by re-reading the source and cross-checked against all 80 existing open issues. Each defect was verified by an independent code-truth skeptic and a novelty skeptic.

[andres9403]

Triage: reproduced with failing regression tests (all 3 defects)

Root cause is as diagnosed in the audit above; all three defects are confirmed present on current base. The audit's cited line numbers have drifted (the file was refactored), so the current-master locations are: Darcy velocity Drying_Model.py:293-296 (k_ra computed but unused, correct form commented at :294); gas species source Drying_Model.py:408-409 (extra /(1 - satur), since epsilon_gas = porosity*(1-satur) at :400); single-node init Drying_Model.py:546 (temp_cond omitted vs. the multi-node branch that includes it), which unit_model then reshapes to 3 + num_comp + num_volatiles columns at :277.

Environment: base c7a38a412d1ab26c8d91ef975ff4a5b3a7393c3b (origin/master), Python 3.10, assimulo 3.4.3. Run: pytest tests/test_drying_gas_balance_issue81.py -m assimulo

Failing tests (self-authored, red on unchanged base) — tests/test_drying_gas_balance_issue81.py

Each test exercises real source: defect 2 via a direct material_balance call; defect 1 by capturing the vel_gas that unit_model hands to material_balance at full saturation; defect 3 by capturing the assembled state width solve_unit produces for a single-node uniform cake. Verified non-vacuous — all three turn green under the minimal fix (uncomment :294 and drop /np.mean(satur); remove the extra /(1 - satur); add temp_cond_init to the single-node tuple).

"""Regression tests for issue #81 (Drying gas balance / Darcy velocity).

Three coupled defects in ``PharmaPy/Drying_Model.py`` share one fix:

1. Darcy gas velocity drops the relative permeability ``k_ra`` and divides by
   ``np.mean(satur)`` instead of using ``k_perm * k_ra * dPg_dz / visc_gas``.
2. The gas-phase species evaporation source (``transfer_gas``) is over-divided
   by an extra ``(1 - satur)`` factor.
3. The single-node / uniform-init branch of ``solve_unit`` omits ``temp_cond``
   from the assembled state vector, so ``unit_model``'s reshape crashes.

Each test exercises real source and asserts the physically correct behavior.
"""
from types import SimpleNamespace

import numpy as np
import pytest

pytest.importorskip("assimulo")
import PharmaPy.Drying_Model as drying_mod
from PharmaPy.Drying_Model import Drying

pytestmark = pytest.mark.assimulo


class _BoundaryReached(Exception):
    """Raised by a stub to stop execution once the value of interest is seen."""


# --------------------------------------------------------------------------- #
# Defect 2: gas-phase species source over-divided by an extra (1 - satur)
# --------------------------------------------------------------------------- #
def test_gas_species_source_has_no_extra_one_minus_satur():
    """``transfer_gas`` must be ``dry_rate / epsilon_gas / rho_gas`` only.

    ``epsilon_gas = porosity * (1 - satur)`` already carries the single gas
    volume-fraction factor, matching the energy balance's ``denom_gas``. The
    active code divides ``transfer_gas`` by an extra ``(1 - satur)``.

    A drying rate placed only on the non-volatile (supercritical) component
    isolates the transfer term: it makes ``dsat_dt == 0`` (volatile rates are
    zero), so ``total_mass_correction`` vanishes and ``dygas_dt`` equals the
    transfer term alone.
    """
    dryer = Drying(number_nodes=2, supercrit_names=["nitrogen"])
    dryer.idx_volatiles = np.array([0, 2])   # component 1 is non-volatile
    dryer.porosity = 0.4
    dryer.rho_liq = np.array([1000.0, 1000.0])
    dryer.dz = np.ones(2)

    satur = np.array([0.5, 0.5])             # 1 - satur = 0.5
    dens_gas = np.array([1.0, 1.0])
    y_gas = np.array([[0.05, 0.90, 0.05],
                      [0.05, 0.90, 0.05]])
    x_liq = np.array([[0.3, 0.7],
                      [0.3, 0.7]])
    # Rate only on component 1 (non-volatile) -> dsat_dt == 0.
    dry_rate = np.array([[0.0, 1.0, 0.0],
                         [0.0, 1.0, 0.0]])
    u_gas = np.zeros(2)                       # kill convection
    inputs = {"mass_frac": np.array([0.05, 0.90, 0.05])}

    _, dygas_dt, _ = dryer.material_balance(
        time=0.0, satur=satur, temp_gas=np.array([320.0, 320.0]),
        temp_sol=np.array([318.0, 318.0]), y_gas=y_gas, x_liq=x_liq,
        u_gas=u_gas, dens_gas=dens_gas, dry_rate=dry_rate, inputs=inputs,
    )

    # dygas_dt is indexed [node, component]; component 1 across both nodes.
    # Correct source for component 1: rate / (porosity*(1-satur)) / rho_gas
    #   = 1.0 / (0.4 * 0.5) / 1.0 = 5.0  (NO extra 1/(1-satur))
    expected_comp1 = np.array([5.0, 5.0])
    np.testing.assert_allclose(dygas_dt[:, 1], expected_comp1)


# --------------------------------------------------------------------------- #
# Defect 1: Darcy gas velocity must vanish at full saturation
# --------------------------------------------------------------------------- #
def test_darcy_gas_velocity_vanishes_at_full_saturation():
    """At S -> 1 the gas relative permeability ``k_ra`` -> 0, so the gas
    superficial velocity must vanish. The active formula drops ``k_ra`` and
    divides by ``np.mean(satur)``, returning the full single-phase velocity.

    ``unit_model`` computes ``vel_gas`` inline and hands it to
    ``material_balance`` as ``u_gas``; we capture that boundary value.
    """
    num_nodes, num_comp, num_volatiles = 2, 3, 2
    dryer = Drying(number_nodes=num_nodes, supercrit_names=["nitrogen"])
    dryer.num_volatiles = num_volatiles
    dryer.idx_volatiles = np.array([0, 2])
    dryer.porosity = 0.4
    dryer.rho_sol = 1500.0
    dryer.s_inf = 0.155
    dryer.dPg_dz = 1.0e5
    dryer.pres_gas = np.full(num_nodes, 101325.0)
    dryer.mass_eta = False
    dryer.CakePhase = SimpleNamespace(alpha=1.0e10)
    # k_perm = 1 / (alpha * rho_sol * (1 - porosity)), as set in solve_unit.
    dryer.k_perm = 1.0 / (dryer.CakePhase.alpha * dryer.rho_sol * (1 - dryer.porosity))
    dryer.Liquid_1 = SimpleNamespace(
        num_species=num_comp, mw=np.array([18.0, 28.0, 46.0]),
        rho_liq=np.array([1000.0, 1000.0, 1000.0]),
    )
    dryer.Vapor_1 = SimpleNamespace(
        mw=np.array([18.0, 28.0, 46.0]),
        getViscosity=lambda temp, mass_frac: np.full(num_nodes, 1.8e-5),
    )

    # Short-circuit the heavy collaborators downstream of the velocity term.
    dryer.get_drying_rate = lambda x_liq, temp_sol, y_gas, p_gas: np.zeros(
        (num_nodes, num_comp))
    dryer.get_inputs = lambda time: {"Inlet": {"mass_frac": np.full(num_comp, 1 / 3)}}

    captured = {}

    def capture_material_balance(time, satur, temp_gas, temp_sol, y_gas,
                                 x_liq, u_gas, dens_gas, dry_rate, inputs):
        captured["u_gas"] = np.asarray(u_gas, dtype=float).copy()
        raise _BoundaryReached

    dryer.material_balance = capture_material_balance

    # Per node: [S, y0, y1, y2, xl0, xl1, Tg, Ts]  (width 3+num_comp+num_volatiles)
    node = np.array([0.999, 0.05, 0.90, 0.05, 0.3, 0.7, 320.0, 318.0])
    states = np.tile(node, (num_nodes, 1))

    with pytest.raises(_BoundaryReached):
        dryer.unit_model(0.0, states.ravel())

    u_gas = captured["u_gas"]
    # Full single-phase velocity here is ~6.2e-4 m/s; correct (k_ra-throttled)
    # velocity at S~1 is ~1e-12. Require it to be numerically negligible.
    assert np.max(np.abs(u_gas)) < 1.0e-6, (
        f"gas velocity should vanish at full saturation, got {u_gas}")


# --------------------------------------------------------------------------- #
# Defect 3: single-node init must include temp_cond in the state vector
# --------------------------------------------------------------------------- #
def test_single_node_init_includes_temp_cond(monkeypatch):
    """A single representative node with uniform liquid composition hits the
    ``len(satur_init) == 1`` branch, which builds
    ``(satur, y_gas, x_liq, temp_gas)`` and omits ``temp_cond``. The assembled
    per-node width must equal ``unit_model``'s reshape width
    ``3 + num_comp + num_volatiles`` (Drying_Model.py, unit_model).
    """
    num_nodes, num_comp = 3, 3
    dryer = Drying(number_nodes=num_nodes, supercrit_names=["nitrogen"])
    dryer.idx_supercrit = np.array([1])          # -> idx_volatiles = [0, 2]
    dryer.names_states_in = ["temp", "mass_frac"]
    dryer.z_centers = np.linspace(0.0, 0.02, num_nodes)
    dryer.cake_height = 0.02

    solid = SimpleNamespace(
        temp=318.0,
        getDensity=lambda: 1500.0,
        getCp=lambda: 900.0,
        getMoments=lambda mom_num: np.array([1.0, 1e-5, 1e-10, 1e-15, 1e-20]),
        x_distrib=np.linspace(1e-6, 1e-4, 10),
        distrib=np.ones(10),
        moments=np.array([1.0, 1e-5, 1e-10, 1e-15, 1e-20]),
    )
    dryer.Solid_1 = solid
    dryer.Liquid_1 = SimpleNamespace(
        num_species=num_comp, mw=np.array([18.0, 28.0, 46.0]),
        getDensity=lambda temp, mass_frac, basis: np.full(num_nodes, 1000.0),
        getSurfTension=lambda temp, mass_frac: np.full(num_nodes, 0.072),
    )
    dryer.Vapor_1 = SimpleNamespace(
        mass_frac=np.array([0.05, 0.90, 0.05]), temp=320.0)
    dryer.CakePhase = SimpleNamespace(
        saturation=np.array([0.8]),              # single node
        Solid_1=solid,
        Liquid_1=SimpleNamespace(mass_frac=np.array([0.3, 0.0, 0.7])),  # 1-D
        z_external=np.array([0.5]),
        alpha=1.0e10, porosity=0.4,
    )

    # Neutralize the irreducible-saturation solver (not under test here).
    monkeypatch.setattr(drying_mod, "get_sat_inf", lambda *a, **k: 0.15)

    captured = {}

    def probe_unit_model(t, y):
        captured["size"] = np.asarray(y).size
        raise _BoundaryReached

    dryer.unit_model = probe_unit_model

    with pytest.raises(_BoundaryReached):
        dryer.solve_unit(deltaP=1.0e5, runtime=10.0)

    num_volatiles = num_comp - 1
    per_node = 3 + num_comp + num_volatiles     # unit_model reshape width
    assert captured["size"] == num_nodes * per_node, (
        f"assembled state has {captured['size']} entries; expected "
        f"{num_nodes * per_node} (per-node width {per_node} incl. temp_cond)")

Red output on unchanged base:

FAILED test_gas_species_source_has_no_extra_one_minus_satur
E   AssertionError: Not equal to tolerance rtol=1e-07, atol=0
E    ACTUAL: array([10., 10.])
E    DESIRED: array([5., 5.])

FAILED test_darcy_gas_velocity_vanishes_at_full_saturation
E   AssertionError: gas velocity should vanish at full saturation, got [0.0006179 0.0006179]
E   assert np.float64(0.0006179018524697537) < 1e-06

FAILED test_single_node_init_includes_temp_cond
E   AssertionError: assembled state has 21 entries; expected 24 (per-node width 8 incl. temp_cond)
E   assert 21 == (3 * 8)

Observed vs expected

Defect	Observed (base)	Expected (fixed)
1 — gas velocity at S≈1 (should →0)	6.18e-4 m/s (full single-phase)	≈0 (k_ra-throttled)
2 — dygas_dt gas source, one component	10	5
3 — assembled single-node state size	21 (per-node width 7)	24 (width 8, incl. temp_cond)

Acceptance criteria

1. Gas velocity uses the relative-permeability form and vanishes as S → 1 while staying finite as S → s_inf (defect-1 test green).
2. Gas species source carries exactly one gas-holdup factor epsilon_gas = porosity*(1-satur), not (1-satur)^2 (defect-2 test green).
3. Single-node/uniform init assembles a per-node width of 3 + num_comp + num_volatiles (includes temp_cond), so unit_model's reshape succeeds (defect-3 test green).
4. Existing test suite still passes.

Next step

An open fix PR already implements all three fixes identically — #119 (fix/issue-81-drying-gas-balance, green CI, currently CHANGES_REQUESTED). So this is not a fresh gh-issue-to-pr target; the tests above are independent verification. Next step is to review/verify #119 against these acceptance criteria and address its outstanding review comments, rather than open a duplicate PR.