diff --git a/tests/test_stream.py b/tests/test_stream.py index 1b18d4d6..744a3af3 100644 --- a/tests/test_stream.py +++ b/tests/test_stream.py @@ -304,6 +304,26 @@ def test_mixture(): Ts = s5.mixture.solve_T_at_SP(phase=s5.phase, mol=s5.mol, S=s5.S, T_guess=s5.T, P=s5.P) assert_allclose(T_expected, Ts, rtol=1e-3) pass + +def test_vle_critical_pure_component(): + N2 = tmo.Chemical('N2', cache=True) + tmo.settings.set_thermo([N2]) + s = tmo.Stream(None, N2=1) + + # Under critical point, liquid is fine + s.vle(T=N2.Tc - 1e-6, P=2 * N2.Pc) + assert s.phase == 'l' + + # Beyond critical point, always model as a gas + s.vle(T=N2.Tc, P=0.5 * N2.Pc) + assert s.phase == 'g' + s.vle(T=N2.Tc, P=2 * N2.Pc) + assert s.phase == 'g' + s.vle(P=N2.Pc, H=s.H + 10) + assert s.phase == 'g' + s.vle(P=N2.Pc, S=s.S + 10) + assert s.phase == 'g' + if __name__ == '__main__': test_stream() @@ -311,5 +331,6 @@ def test_mixture(): test_vlle() stream_methods() test_stream_property_cache() + test_vle_critical_pure_component() test_critical() test_mixture() \ No newline at end of file diff --git a/thermosteam/equilibrium/bubble_point.py b/thermosteam/equilibrium/bubble_point.py index ba3d48e4..d47ae139 100644 --- a/thermosteam/equilibrium/bubble_point.py +++ b/thermosteam/equilibrium/bubble_point.py @@ -197,7 +197,8 @@ def solve_Ty(self, z, P): if N == 0: raise ValueError('no components present') if N == 1: - T = self.chemicals[fn.first_true_index(positives)].Tsat(P, check_validity=False) + chemical = self.chemicals[fn.first_true_index(positives)] + T = chemical.Tsat(P, check_validity=False) if P < chemical.Pc else chemical.Tc y = z.copy() else: f = self._T_error @@ -251,7 +252,8 @@ def solve_Py(self, z, T): if N == 0: raise ValueError('no components present') if N == 1: - P = self.chemicals[fn.first_true_index(positives)].Psat(T) + chemical = self.chemicals[fn.first_true_index(positives)] + P = chemical.Psat(T) if T < chemical.Tc else chemical.Pc y = z.copy() else: if T > self.Tmax: T = self.Tmax diff --git a/thermosteam/equilibrium/vle.py b/thermosteam/equilibrium/vle.py index 605b255d..5b0e25e7 100644 --- a/thermosteam/equilibrium/vle.py +++ b/thermosteam/equilibrium/vle.py @@ -484,15 +484,20 @@ def thermal_condition(self): ### Single component equilibrium case ### def _set_thermal_condition_chemical(self, T, P): - # Either liquid or gas - Psat = self._chemical.Psat(T) - tol = 1e-3 - if P < Psat - tol: + chemical = self._chemical + if T >= chemical.Tc: self._liquid_mol[self._index] = 0 self._vapor_mol[self._index] = self._mol_vle - elif P > Psat + tol: - self._liquid_mol[self._index] = self._mol_vle - self._vapor_mol[self._index] = 0 + else: + # Either liquid or gas + Psat = chemical.Psat(T) + tol = 1e-3 + if P < Psat - tol: + self._liquid_mol[self._index] = 0 + self._vapor_mol[self._index] = self._mol_vle + elif P > Psat + tol: + self._liquid_mol[self._index] = self._mol_vle + self._vapor_mol[self._index] = 0 def _set_TV_chemical(self, T, V): # Set vapor fraction @@ -514,9 +519,13 @@ def _set_PH_chemical(self, P, H): thermo = self._thermo phase_data = self._phase_data mixture = thermo.mixture + chemical = self._chemical # Set temperature in equilibrium - self._T = self._thermal_condition.T = T = self._chemical.Tsat(P, check_validity=False) + self._T = self._thermal_condition.T = T = ( + # Maximally the critical temperature in equilibrium + chemical.Tsat(P, check_validity=False) if P < chemical.Pc else chemical.Tc + ) # Check if super heated vapor vapor_mol[index] = mol @@ -555,7 +564,8 @@ def _set_TH_chemical(self, T, H): liquid_mol[index] = 0 H_dew = mixture.xH(phase_data, T, P) if H >= H_dew: - self._thermal_condition.T = mixture.xsolve_T_at_HP(phase_data, H, T, P) + raise NotImplementedError('cannot solve for pressure yet') + self._thermal_condition.P = mixture.xsolve_P_at_HT(phase_data, H, T, P) return # Check if subcooled liquid @@ -563,7 +573,8 @@ def _set_TH_chemical(self, T, H): liquid_mol[index] = mol H_bubble = mixture.xH(phase_data, T, P) if H <= H_bubble: - self._thermal_condition.T = mixture.xsolve_T_at_HP(phase_data, H, T, P) + raise NotImplementedError('cannot solve for pressure yet') + self._thermal_condition.P = mixture.xsolve_P_at_HT(phase_data, H, T, P) return # Adjust vapor fraction accordingly @@ -579,9 +590,13 @@ def _set_PS_chemical(self, P, S): thermo = self._thermo phase_data = self._phase_data mixture = thermo.mixture + chemical = self._chemical # Set temperature in equilibrium - self._T = self._thermal_condition.T = T = self._chemical.Tsat(P, check_validity=False) + self._T = self._thermal_condition.T = T = ( + # Maximally the critical temperature in equilibrium + chemical.Tsat(P, check_validity=False) if P < chemical.Pc else chemical.Tc + ) # Check if super heated vapor vapor_mol[index] = mol @@ -620,7 +635,8 @@ def _set_TS_chemical(self, T, S): liquid_mol[index] = 0 S_dew = mixture.xS(phase_data, T, P) if S >= S_dew: - self._thermal_condition.T = mixture.xsolve_T_at_SP(phase_data, S, T, P) + raise NotImplementedError('cannot solve for pressure yet') + self._thermal_condition.P = mixture.xsolve_P_at_ST(phase_data, S, T, P) return # Check if subcooled liquid @@ -628,7 +644,8 @@ def _set_TS_chemical(self, T, S): liquid_mol[index] = mol S_bubble = mixture.xS(phase_data, T, P) if S <= S_bubble: - self._thermal_condition.T = mixture.xsolve_T_at_SP(phase_data, S, T, P) + raise NotImplementedError('cannot solve for pressure yet') + self._thermal_condition.P = mixture.xsolve_P_at_ST(phase_data, S, T, P) return # Adjust vapor fraction accordingly