Skip to content

Commit 1f0d036

Browse files
yoelcortesyoelcortes
committed
fix edge cases and improve plotting
1 parent 3a47254 commit 1f0d036

File tree

5 files changed

+88
-53
lines changed

5 files changed

+88
-53
lines changed

‎thermosteam/_settings.py‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,7 @@ def set_thermo(self, thermo: tmo.Thermo|Iterable[str|tmo.Chemical],
115115
skip_checks: Optional[bool]=False,
116116
ideal: Optional[bool]=False,
117117
db: Optional[str]='default',
118+
pkg=None,
118119
):
119120
"""
120121
Set the default :class:`~thermosteam.Thermo` object. If `thermo` is
@@ -140,11 +141,13 @@ def set_thermo(self, thermo: tmo.Thermo|Iterable[str|tmo.Chemical],
140141
algorithms.
141142
db : str, optional
142143
Database to load any chemicals.
144+
pkg : str, optional
145+
Name of property package. Defaults to 'Dortmund-UNIFAC'.
143146
144147
"""
145148
if not isinstance(thermo, (tmo.Thermo, tmo.IdealThermo)):
146149
thermo = tmo.Thermo(thermo, mixture=mixture, cache=cache, skip_checks=skip_checks,
147-
Gamma=Gamma, Phi=Phi, PCF=PCF, db=db)
150+
Gamma=Gamma, Phi=Phi, PCF=PCF, db=db, pkg=pkg)
148151
if ideal: thermo = thermo.ideal()
149152
self._thermo = thermo
150153

‎thermosteam/equilibrium/bubble_point.py‎

Lines changed: 13 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,11 @@ class BubblePoint:
115115
maxiter = 100
116116
T_tol = 1e-12
117117
P_tol = 1e-6
118+
Tmin_factor = 0.1
119+
Tmax_factor = 10
120+
Pmin_factor = 0.1
121+
Pmax_factor = 10
122+
118123
def __new__(cls, chemicals=None, thermo=None):
119124
thermo = settings.get_default_thermo(thermo)
120125
if chemicals is None:
@@ -276,7 +281,8 @@ def solve_Ty(self, z, P, liquid_conversion=None, guess=None):
276281
checkiter=False,
277282
maxiter=self.maxiter)
278283
except RuntimeError:
279-
Tmin = self.Tmin; Tmax = self.Tmax
284+
Tmin = max(self.Tmin_factor * T_guess, self.Tmin)
285+
Tmax = min(self.Tmax_factor * T_guess, self.Tmax)
280286
T = flx.IQ_interpolation(f, Tmin, Tmax,
281287
f(Tmin, *args), f(Tmax, *args),
282288
T_guess, self.T_tol, 5e-12, args,
@@ -296,7 +302,8 @@ def solve_Ty(self, z, P, liquid_conversion=None, guess=None):
296302
self.T_tol, 5e-12, args,
297303
checkiter=False, maxiter=self.maxiter)
298304
except RuntimeError:
299-
Tmin = self.Tmin; Tmax = self.Tmax
305+
Tmin = max(self.Tmin_factor * T_guess, self.Tmin)
306+
Tmax = min(self.Tmax_factor * T_guess, self.Tmax)
300307
T = flx.IQ_interpolation(f, Tmin, Tmax,
301308
f(Tmin, *args), f(Tmax, *args),
302309
T_guess, self.T_tol, 5e-12, args,
@@ -361,7 +368,8 @@ def solve_Py(self, z, T, liquid_conversion=None):
361368
P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9,
362369
args, checkiter=False, maxiter=self.maxiter)
363370
except RuntimeError:
364-
Pmin = self.Pmin; Pmax = self.Pmax
371+
Pmin = max(self.Pmin_factor * P_guess, self.Pmin)
372+
Pmax = min(self.Pmax_factor * P_guess, self.Pmax)
365373
P = flx.IQ_interpolation(f, Pmin, Pmax,
366374
f(Pmin, *args), f(Pmax, *args),
367375
P_guess, self.P_tol, 5e-12, args,
@@ -382,7 +390,8 @@ def solve_Py(self, z, T, liquid_conversion=None):
382390
args, checkiter=False,
383391
maxiter=self.maxiter)
384392
except RuntimeError:
385-
Pmin = self.Pmin; Pmax = self.Pmax
393+
Pmin = max(self.Pmin_factor * P_guess, self.Pmin)
394+
Pmax = min(self.Pmax_factor * P_guess, self.Pmax)
386395
P = flx.IQ_interpolation(f, Pmin, Pmax,
387396
f(Pmin, *args), f(Pmax, *args),
388397
P_guess, self.P_tol, 5e-12, args,

‎thermosteam/equilibrium/dew_point.py‎

Lines changed: 13 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -120,6 +120,11 @@ class DewPoint:
120120
maxiter = 50
121121
T_tol = 1e-9
122122
P_tol = 1e-3
123+
Tmin_factor = 0.1
124+
Tmax_factor = 10
125+
Pmin_factor = 0.1
126+
Pmax_factor = 10
127+
123128
def __new__(cls, chemicals=None, thermo=None):
124129
thermo = settings.get_default_thermo(thermo)
125130
if chemicals is None:
@@ -283,8 +288,8 @@ def solve_Tx(self, z, P, gas_conversion=None, guess=None):
283288
maxiter=self.maxiter,
284289
checkiter=False)
285290
except RuntimeError:
286-
Tmin = self.Tmin
287-
Tmax = self.Tmax
291+
Tmin = max(self.Tmin_factor * T_guess0, self.Tmin)
292+
Tmax = min(self.Tmax_factor * T_guess0, self.Tmax)
288293
T = flx.IQ_interpolation(f, Tmin, Tmax,
289294
f(Tmin, *args), f(Tmax, *args),
290295
T_guess0, self.T_tol, 5e-12, args,
@@ -304,7 +309,8 @@ def solve_Tx(self, z, P, gas_conversion=None, guess=None):
304309
self.T_tol, 5e-12, args,
305310
checkiter=False)
306311
except RuntimeError:
307-
Tmin = self.Tmin; Tmax = self.Tmax
312+
Tmin = max(self.Tmin_factor * T_guess, self.Tmin)
313+
Tmax = min(self.Tmax_factor * T_guess, self.Tmax)
308314
T = flx.IQ_interpolation(f, Tmin, Tmax,
309315
f(Tmin, *args), f(Tmax, *args),
310316
T_guess, self.T_tol, 5e-12, args,
@@ -364,8 +370,8 @@ def solve_Px(self, z, T, gas_conversion=None, guess=None):
364370
P = flx.aitken_secant(f, P_guess0, P_guess1, self.P_tol, 5e-12, args,
365371
checkiter=False, maxiter=self.maxiter)
366372
except RuntimeError:
367-
Pmin = self.Pmin
368-
Pmax = self.Pmax
373+
Pmin = max(self.Pmin_factor * P_guess0, self.Pmin)
374+
Pmax = min(self.Pmax_factor * P_guess0, self.Pmax)
369375
P = flx.IQ_interpolation(f, Pmin, Pmax,
370376
f(Pmin, *args), f(Pmax, *args),
371377
P_guess0, self.P_tol, 5e-12, args,
@@ -385,7 +391,8 @@ def solve_Px(self, z, T, gas_conversion=None, guess=None):
385391
P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9,
386392
args, checkiter=False)
387393
except RuntimeError:
388-
Pmin = self.Pmin; Pmax = self.Pmax
394+
Pmin = max(self.Pmin_factor * P_guess, self.Pmin)
395+
Pmax = min(self.Pmax_factor * P_guess, self.Pmax)
389396
P = flx.IQ_interpolation(f, Pmin, Pmax,
390397
f(Pmin, *args), f(Pmax, *args),
391398
P_guess, self.P_tol, 5e-12, args,

‎thermosteam/equilibrium/plot_equilibrium.py‎

Lines changed: 39 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
import flexsolve as flx
2222
from scipy import interpolate
2323
from scipy.ndimage.filters import gaussian_filter
24+
from typing import Iterable
2425

2526
__all__ = (
2627
'plot_vle_binary_phase_envelope',
@@ -59,17 +60,17 @@ def as_thermo(thermo, chemicals): # pragma: no cover
5960
# %% Plot functions
6061

6162
def plot_vle_binary_phase_envelope(chemicals, T=None, P=None, vc=None, lc=None, thermo=None,
62-
yticks=None): # pragma: no cover
63+
yticks=None, line_styles=None, N=None): # pragma: no cover
6364
"""
6465
Plot the binary phase envelope of two chemicals at a given temperature or pressure.
6566
6667
Parameters
6768
----------
6869
chemicals : Iterable[Chemical or str]
6970
Chemicals in equilibrium.
70-
T : float, optional
71+
T : float|list[float], optional
7172
Temperature [K].
72-
P : float, optional
73+
P : float|list[float], optional
7374
Pressure [Pa].
7475
vc : str, optional
7576
Color of vapor line.
@@ -86,53 +87,58 @@ def plot_vle_binary_phase_envelope(chemicals, T=None, P=None, vc=None, lc=None,
8687
.. figure:: ../images/water_ethanol_binary_phase_envelope.png
8788
8889
"""
90+
if line_styles is None: line_styles = ['-', '-.', '--']
8991
thermo = as_thermo(thermo, chemicals)
9092
chemical_a, chemical_b = chemicals = [thermo.as_chemical(i) for i in chemicals]
9193
BP = tmo.BubblePoint(chemicals, thermo)
92-
N = 50
94+
if N is None: N = 50
9395
zs_a = np.linspace(0, 1, N)
9496
zs_b = 1 - zs_a
9597
zs = np.vstack([zs_a, zs_b]).transpose()
9698
if P:
99+
if not isinstance(P, Iterable): P = [P]
100+
Xs = P
101+
X_units = 'Pa'
97102
assert not T, "must pass either T or P, but not both"
98-
bps = [BP(z, P=P) for z in zs]
99-
ms = [bp.T for bp in bps]
103+
bpss = [[BP(z, P=i) for z in zs] for i in P]
104+
mss = [[bp.T for bp in bps] for bps in bpss]
100105
ylabel = 'Temperature [K]'
101106
elif T:
107+
if not isinstance(T, Iterable): T = [T]
102108
assert not P, "must pass either T or P, but not both"
103-
bps = [BP(z, T=T) for z in zs]
104-
ms = [bp.P for bp in bps]
109+
Xs = T
110+
X_units = 'K'
111+
bpss = [[BP(z, T=i) for z in zs] for i in T]
112+
mss = [[bp.P for bp in bps] for bps in bpss]
105113
ylabel = 'Pressure [Pa]'
106114
else:
107115
raise AssertionError("must pass either T or P")
108-
ms = np.array(ms)
109-
ys_a = np.array([bp.y[0] for bp in bps])
110-
ms_liq = ms.copy()
111-
ms_gas = interpolate.interp1d(ys_a, ms, bounds_error=False, kind='slinear')(zs_a)
112-
if P:
113-
azeotrope = ms_liq > ms_gas
114-
elif T:
115-
azeotrope = ms_liq < ms_gas
116-
ms_liq[0] = 0.5 * (ms_liq[0] + ms_gas[0])
117-
ms_liq[-1] = 0.5 * (ms_liq[-1] + ms_gas[-1])
118-
if azeotrope.any():
119-
index = np.where(azeotrope)[0]
120-
left = index[0]
121-
right = index[-1]
122-
mid = int(np.median(index))
123-
azeotrope[left:right] = True
124-
azeotrope[mid] = False
125-
ms_gas[mid] = ms_liq[mid]
126-
ms_gas = interpolate.interp1d(zs_a[~azeotrope], ms_gas[~azeotrope], bounds_error=False, kind='slinear')(zs_a)
127-
128-
top, bottom = (chemical_a, chemical_b) if ys_a.mean() > 0.5 else (chemical_a, chemical_b)
129116
plt.figure()
117+
mss = np.array(mss)
118+
for X, ms, bps, ls in zip(Xs, mss, bpss, line_styles):
119+
ys_a = np.array([bp.y[0] for bp in bps])
120+
ms_liq = ms.copy()
121+
ms_gas = interpolate.interp1d(ys_a, ms, bounds_error=False, kind='slinear')(zs_a)
122+
# if P:
123+
# azeotrope = ms_liq > ms_gas
124+
# elif T:
125+
# azeotrope = ms_liq < ms_gas
126+
# ms_liq[0] = 0.5 * (ms_liq[0] + ms_gas[0])
127+
# ms_liq[-1] = 0.5 * (ms_liq[-1] + ms_gas[-1])
128+
# if azeotrope.any():
129+
# index = np.where(azeotrope)[0]
130+
# left = index[0]
131+
# right = index[-1]
132+
# mid = int(np.median(index))
133+
# azeotrope[left:right] = True
134+
# azeotrope[mid] = False
135+
# ms_gas[mid] = ms_liq[mid]
136+
# ms_gas = interpolate.interp1d(zs_a[~azeotrope], ms_gas[~azeotrope], bounds_error=False, kind='slinear')(zs_a)
137+
plt.plot(zs_a, ms_gas, c=vc if vc is not None else colors.red.RGBn, label=f'vapor [{int(X)} {X_units}]', ls=ls)
138+
plt.plot(zs_a, ms_liq, c=lc if lc is not None else colors.blue.RGBn, label=f'liquid [{int(X)} {X_units}]', ls=ls)
130139
plt.xlim([0, 1])
131-
plt.plot(zs_a, ms_gas, c=vc if vc is not None else colors.red.RGBn, label='vapor')
132-
plt.plot(zs_a, ms_liq, c=lc if lc is not None else colors.blue.RGBn, label='liquid')
133-
plt.ylim([ms.min(), ms.max()])
140+
plt.ylim([mss.min(), mss.max()])
134141
if yticks is None: yticks, ytext = plt.yticks()
135-
136142
plt.legend()
137143
plt.xlabel(f'{chemical_a} molar fraction')
138144
plt.ylabel(ylabel)

‎thermosteam/mixture/__init__.py‎

Lines changed: 19 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,8 @@
2727
from .._chemicals import Chemical, CompiledChemicals, chemical_data_array
2828

2929
__all__ = ('Mixture', 'IdealMixture')
30-
30+
31+
3132
# %% Convenience for EOS
3233

3334
def get_excess_property(eos, free_energy, phase):
@@ -570,9 +571,10 @@ def __init__(self, chemicals, eos_chemical_index,
570571
# Ensure consistent equation of state
571572
eos_chemicals = [i for _, i in eos_chemical_index.values()]
572573
eos_mix = self.eos_args('g', np.ones(len(eos_chemicals)), 298.15, 101325)[0]
573-
for chemical, eos in zip(eos_chemicals, eos_mix.pures()):
574-
chemical._eos = eos
575-
chemical.reset_free_energies()
574+
if eos_mix is not None:
575+
for chemical, eos in zip(eos_chemicals, eos_mix.pures()):
576+
chemical._eos = eos
577+
chemical.reset_free_energies()
576578

577579
# Populate excess energies
578580
for name in ('H_excess', 'S_excess'):
@@ -655,16 +657,21 @@ def eos_args(self, phase, mol, T, P):
655657
kijs = IPDB.get_ip_asymmetric_matrix(self.chemsep_db, data.CASs, 'kij')
656658
except:
657659
kijs = None
658-
self.cache[key] = eos = self.EOS(
659-
Tcs=data.Tcs, Pcs=data.Pcs, omegas=data.omegas, kijs=kijs,
660-
T=T, P=P, zs=zs, only_g=only_g, only_l=only_l,
661-
fugacities=False
662-
)
660+
if key == ():
661+
eos = None
662+
else:
663+
eos = self.EOS(
664+
Tcs=data.Tcs, Pcs=data.Pcs, omegas=data.omegas, kijs=kijs,
665+
T=T, P=P, zs=zs, only_g=only_g, only_l=only_l,
666+
fugacities=False
667+
)
668+
self.cache[key] = eos
663669
return eos, index, eos_index, eos_mol
664670

665671
def Hvap(self, mol, T, P):
666672
phase = 'l'
667673
eos, _, _, eos_mol = self.eos_args(phase, mol, T, P)
674+
if eos is None: return 0
668675
return eos.Hvap(T) * eos_mol
669676

670677
def dh_dep_dzs(self, phase, mol, T, P):
@@ -694,6 +701,7 @@ def Cn(self, phase, mol, T, P=101325):
694701
eos, index, eos_index, eos_mol = self.eos_args(
695702
phase, mol, T, P
696703
)
704+
if eos is None: return Cn
697705
excess_ref = sum([
698706
mol[eos_index[i]] * get_excess_property(
699707
eos.to_TPV_pure(T=T, P=101325, V=None, i=i),
@@ -715,6 +723,7 @@ def H(self, phase, mol, T, P):
715723
eos, _, eos_index, eos_mol = self.eos_args(
716724
phase, mol, T, P
717725
)
726+
if eos is None: return H
718727
H_dep = get_excess_property(eos, 'H', phase)
719728
H_excess = self.H_excess[phase]
720729
H_ref = 0
@@ -733,6 +742,7 @@ def S(self, phase, mol, T, P):
733742
eos, _, eos_index, eos_mol = self.eos_args(
734743
phase, mol, T, P
735744
)
745+
if eos is None: return S
736746
S_dep = get_excess_property(eos, 'S', phase)
737747
S_excess = self.S_excess[phase]
738748
S_ref = 0

0 commit comments

Comments
 (0)