Fundamental flaws in Gibbs and VCS algorithms for constant enthalpy / energy / entropy

Start with MultiPhase::equilibrate_MultiPhaseEquil and trace how the Gibbs and VCS solvers handle HP, TV, and related constraints. Reproduce the ice/liquid/gas example from the issue, then determine how convergence or phase-mole handling should account for discontinuous enthalpy; done means the case no longer returns an incorrect state silently and the reported failure is meaningful.

Both the "Gibbs" and "VCS" multiphase equilibrium solvers directly implement equilibrium at constant temperature and pressure, then wrap that in an outer iteration for all of the other property pairs (e.g. HP, TV, etc.) that they support to find the value of T and P that satisfy the constraint. This fails in any case where the equilibrium state has multiple phases in equilibrium.

The following example demonstrates this for a ice / liquid water / gas system in which burning hydrogen should melt some but not all of the ice. This will print:

• the initial state
• the equilibrium state according to the VCS solver
• the correct equilibrium state computed manually, assuming the gas-phase equilibrium is correct

gas = ct.Solution('h2o2.xml')
solid = ct.Solution('water.cti' ,'ice')
liquid = ct.Solution('water.cti', 'liquid_water')

# initial mole numbers
n0_gas = 1
n0_liquid = 0
n0_solid = 4

def setup():
    gas.TPX = 270, 101325, 'H2:0.04, O2:0.02, AR:0.94'
    solid.TP = 270, 101325
    liquid.TP = 270, 101325
    return ct.Mixture([(gas, n0_gas), (liquid, n0_liquid), (solid, n0_solid)])

def enthalpy(mix):
    ng, nl, ns = mix.phase_moles()
    return (ng * gas.enthalpy_mole
            + nl * liquid.enthalpy_mole
            + ns * solid.enthalpy_mole)

def print_state(mix):
    print('T = {}; H = {}; moles = {}'.format(mix.T, enthalpy(mix), mix.phase_moles()))

mix = setup()
H0 = enthalpy(mix)
print_state(mix)

mix.equilibrate('HP', solver='vcs')
print_state(mix)

N = sum(mix.phase_moles())
hl = liquid.enthalpy_mole
hg = gas.enthalpy_mole
hs = solid.enthalpy_mole

# solve for correct mole numbers:
# H0 = hg * ng + hl * nl + hs * ns
# N = ng + nl + ns
ng = mix.phase_moles()[0]
nl = (H0 - hg*ng - hs*(N-ng)) / (hl - hs)
ns = N - nl - ng

mix.set_phase_moles(gas, ng)
mix.set_phase_moles(solid, ns)
mix.set_phase_moles(liquid, nl)
print_state(mix)

The output from this code is:

initial state: T = 270.0; H = -1175977302.5; moles = [1.0, 0.0, 4.0]
VCS state: T = 273.144066429; H = -1162658273.81; moles = [0.9457838, 4.034216186 0.0]
correct state: T = 273.144066429; H = -1175977302.5; moles = [0.9457838, 1.81260397, 2.2216122]

The essential problem is this: An iterative method which only adjusts temperature cannot find this state because the enthalpy is not a continuous function of T, and the solver will always put all of the moles into either the liquid or solid phase as it jumps across the discontinuity. The only way to find this state is if the actual enthalpy constraint is solved simultaneously with the equations to determine the number of moles of each species in each phase.

At a minimum, this situation should raise an exception. Currently, the VCS solver returns an incorrect solution without any indication that anything is wrong. The Gibbs solver raises the following exception:

RuntimeError: 
***********************************************************************
CanteraError thrown by MultiPhaseEquil:
condensed-phase speciesH2O(S) is excluded since its thermo properties are 
not valid at this temperature, but it has non-zero moles in the initial state.
***********************************************************************

Which is in some ways an improvement, but also unhelpful -- "this temperature" is not the initial temperature of the Mixture. If we modify the input thermo data to claim that it is valid over an extended range, and add some debug output to `MultiPhase::equilibrate_MultiPhaseEquil to print out the temperature and enthalpy bounds of the search space at each iteration, we can see the failure mode:

Trying MultiPhaseEquil (Gibbs) equilibrium solver
n = 0, Tlow = 150.0000, Thigh = 7000.0000, Hlow = -9.9912e+02, Hhigh = -9.9912e+02
n = 1, Tlow = 270.0000, Thigh = 7000.0000, Hlow = -1.1875e+09, Hhigh = -9.9912e+02
n = 2, Tlow = 270.0000, Thigh = 1374.7727, Hlow = -1.1875e+09, Hhigh = 1.9124e+10
n = 3, Tlow = 270.0000, Thigh = 822.3864, Hlow = -1.1875e+09, Hhigh = -8.9048e+08
n = 4, Tlow = 270.0000, Thigh = 546.1932, Hlow = -1.1875e+09, Hhigh = -9.3753e+08
n = 5, Tlow = 270.0000, Thigh = 408.0966, Hlow = -1.1875e+09, Hhigh = -9.5975e+08
n = 6, Tlow = 270.0000, Thigh = 339.0483, Hlow = -1.1875e+09, Hhigh = -1.1276e+09
n = 7, Tlow = 270.0000, Thigh = 304.5241, Hlow = -1.1875e+09, Hhigh = -1.1507e+09
n = 8, Tlow = 270.0000, Thigh = 287.2621, Hlow = -1.1875e+09, Hhigh = -1.1576e+09
n = 9, Tlow = 270.0000, Thigh = 278.6310, Hlow = -1.1875e+09, Hhigh = -1.1607e+09
n = 10, Tlow = 270.0000, Thigh = 274.3155, Hlow = -1.1875e+09, Hhigh = -1.1622e+09
n = 11, Tlow = 272.1578, Thigh = 274.3155, Hlow = -1.1870e+09, Hhigh = -1.1622e+09
n = 12, Tlow = 273.1205, Thigh = 274.3155, Hlow = -1.1868e+09, Hhigh = -1.1622e+09
n = 13, Tlow = 273.1205, Thigh = 274.3155, Hlow = -1.1868e+09, Hhigh = -1.1622e+09
n = 14, Tlow = 273.1205, Thigh = 273.6486, Hlow = -1.1868e+09, Hhigh = -1.1625e+09
n = 15, Tlow = 273.1205, Thigh = 273.3846, Hlow = -1.1868e+09, Hhigh = -1.1626e+09
n = 16, Tlow = 273.1205, Thigh = 273.2525, Hlow = -1.1868e+09, Hhigh = -1.1626e+09
n = 17, Tlow = 273.1205, Thigh = 273.1865, Hlow = -1.1868e+09, Hhigh = -1.1626e+09
n = 18, Tlow = 273.1205, Thigh = 273.1535, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 19, Tlow = 273.1205, Thigh = 273.1535, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 20, Tlow = 273.1370, Thigh = 273.1535, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 21, Tlow = 273.1370, Thigh = 273.1535, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 22, Tlow = 273.1370, Thigh = 273.1444, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 23, Tlow = 273.1370, Thigh = 273.1444, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 24, Tlow = 273.1407, Thigh = 273.1444, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 25, Tlow = 273.1424, Thigh = 273.1444, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 26, Tlow = 273.1433, Thigh = 273.1444, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 27, Tlow = 273.1438, Thigh = 273.1444, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 28, Tlow = 273.1438, Thigh = 273.1444, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 29, Tlow = 273.1438, Thigh = 273.1441, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 30, Tlow = 273.1438, Thigh = 273.1441, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 31, Tlow = 273.1440, Thigh = 273.1441, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 32, Tlow = 273.1440, Thigh = 273.1441, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 33, Tlow = 273.1440, Thigh = 273.1441, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 34, Tlow = 273.1440, Thigh = 273.1441, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 35, Tlow = 273.1440, Thigh = 273.1441, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 36, Tlow = 273.1440, Thigh = 273.1441, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
n = 37, Tlow = 273.1441, Thigh = 273.1441, Hlow = -1.1868e+09, Hhigh = -1.1627e+09
(repeat through n = 100)

***********************************************************************
CanteraError thrown by MultiPhase::equilibrate_MultiPhaseEquil:
No convergence for T
***********************************************************************

So at least the Gibbs solver knows that it hasn't found the solution.