Skip to content

Commit

Permalink
account for critical point in pure component vle
Browse files Browse the repository at this point in the history
  • Loading branch information
@yoelcortes
yoelcortes committed Aug 29, 2022
1 parent 7bbb637 commit 15c22cc
Show file tree
Hide file tree
Showing 3 changed files with 55 additions and 15 deletions.
21 changes: 21 additions & 0 deletions tests/test_stream.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,12 +304,33 @@ 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()
test_multistream()
test_vlle()
stream_methods()
test_stream_property_cache()
test_vle_critical_pure_component()
test_critical()
test_mixture()
6 changes: 4 additions & 2 deletions thermosteam/equilibrium/bubble_point.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
43 changes: 30 additions & 13 deletions thermosteam/equilibrium/vle.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -555,15 +564,17 @@ 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
vapor_mol[index] = 0
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
Expand All @@ -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
Expand Down Expand Up @@ -620,15 +635,17 @@ 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
vapor_mol[index] = 0
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
Expand Down

1 comment on commit 15c22cc

@yoelcortes
Copy link
Member Author

@yoelcortes yoelcortes commented on 15c22cc Aug 29, 2022 •
edited

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@BenPortner, after VLE, a pure component stream beyond the critical temperature is always a gas now.

Please sign in to comment.