begin adding reactive equilibrium

tests/test_reaction.py

@@ -214,6 +214,28 @@ def test_reaction_enthalpy_with_phases():
                               correct_atomic_balance=True)
     assert_allclose(combustion.dH, -2787650.3239119546, rtol=1e-3)
 
+def test_reactive_phase_equilibrium():
+    tmo.settings.set_thermo(['EthylLactate', 'LacticAcid', 'H2O', 'Ethanol'], cache=True)
+    rxn = tmo.Reaction('LacticAcid + Ethanol -> H2O + EthylLactate', reactant='LacticAcid', X=0.2)
+    stream = tmo.Stream(
+        H2O=1, Ethanol=5, LacticAcid=1
+    )
+    stream.vle(T=360, P=101325, liquid_reaction=rxn)
+    assert_allclose(
+        stream.imol['g'], 
+        [0.0,
+         0.1678023002166189, 
+         0.4995276555980901, 
+         3.010149544370034]
+    )
+    assert_allclose(
+        stream.imol[''], 
+        [0.1664395399566762,
+         0.6657581598267048,
+         0.6669118843585862,
+         1.8234109156732896]
+    )
+
 def test_repr():
     cal2joule = 4.184
     Glucan = tmo.Chemical('Glucan', search_db=False, formula='C6H10O5', Hf=-233200*cal2joule, phase='s', default=True)

thermosteam/_stream.py

@@ -41,6 +41,7 @@
 mass_units = indexer.ChemicalMassFlowIndexer.units
 vol_units = indexer.ChemicalVolumetricFlowIndexer.units
 
+
 class StreamData:
     __slots__ = ('_imol', '_T', '_P', '_phases')
 
@@ -49,8 +50,24 @@ def __init__(self, imol, thermal_condition, phases):
         self._T = thermal_condition._T
         self._P = thermal_condition._P
         self._phases = phases
+
 
+class TemporaryPhase:
+    __slots__ = ('stream', 'original', 'temporary')
 
+    def __init__(self, stream, original, temporary):
+        self.stream = stream
+        self.original = original
+        self.temporary = temporary
+
+    def __enter__(self):
+        self.stream._phase._phase = self.temporary
+
+    def __exit__(self, type, exception, traceback):
+        self.stream._phase._phase = self.original
+        if exception: raise exception
+
+
 # %%
 
 @utils.units_of_measure(UofM.stream_units_of_measure)
@@ -335,6 +352,9 @@ def __init__(self, ID: Optional[str]='',
             data = self._imol.data
             self.phases = [j for i, j in enumerate(self.phases) if data[i].any()]
 
+    def temporary_phase(self, phase):
+        return TemporaryPhase(self, self.phase, phase)
+
     @classmethod
     def from_data(cls, data, ID=None, price=0., characterization_factors=None, thermo=None):
         self = cls.__new__(cls)
@@ -347,6 +367,12 @@ def from_data(cls, data, ID=None, price=0., characterization_factors=None, therm
         self.set_data(data)
         return self
 
+    def __len__(self):
+        return 1
+
+    def __iter__(self):
+        yield self
+
     def __getitem__(self, key):
         phase = self.phase
         if key.lower() == phase.lower(): return self
@@ -355,6 +381,7 @@ def __getitem__(self, key):
     def __reduce__(self):
         return self.from_data, (self.get_data(), self._ID, self._price, self.characterization_factors, self._thermo)
 
+    # Phenomena-oriented simulation
     def equation(self, variable, f=None):
         if f is None: return lambda f: self.equation(variable, f)
         equations = self.equations
@@ -370,6 +397,14 @@ def _create_linear_equations(self, variable):
         else:
             return []
 
+    def _update_energy_departure_coefficient(self, coefficients):
+        source = self.source
+        if source is None: return
+        coeff = source._get_energy_departure_coefficient(self)
+        if coeff is None: return
+        key, value = coeff
+        coefficients[key] = value
+
     def _get_decoupled_variable(self, variable): pass
 
     def _update_decoupled_variable(self, variable, value, index=None):