Skip to content

Commit 8d4dced

Browse files
HugoMValeHugoMVale
committed
Enhance equations of state with thermal expansion and compressibility methods
1 parent 63be42b commit 8d4dced

File tree

5 files changed

+268
-16
lines changed

5 files changed

+268
-16
lines changed

‎src/polykin/thermo/eos/base.py‎

Lines changed: 156 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -10,10 +10,11 @@
1010
from numpy import sqrt
1111
from scipy.constants import R
1212

13+
from polykin.math.derivatives import derivative_complex
1314
from polykin.utils.math import eps
14-
from polykin.utils.types import FloatVector
15+
from polykin.utils.types import FloatVector, Number
1516

16-
__all__ = ["EoS"]
17+
__all__ = ["EoS", "GasEoS"]
1718

1819

1920
class EoS(ABC):
@@ -104,7 +105,7 @@ def P(self, T: float, v: float, y: FloatVector) -> float:
104105
pass
105106

106107
@abstractmethod
107-
def Z(self, T: float, P: float, y: FloatVector) -> float:
108+
def Z(self, T: Number, P: Number, y: FloatVector) -> float:
108109
"""Calculate the compressibility factor of the fluid."""
109110
pass
110111

@@ -143,6 +144,78 @@ def v(
143144
"""
144145
return self.Z(T, P, y) * R * T / P
145146

147+
def beta(
148+
self,
149+
T: float,
150+
P: float,
151+
y: FloatVector,
152+
) -> float:
153+
r"""Calculate the thermal expansion coefficient.
154+
155+
$$ \beta
156+
= \frac{1}{v} \left( \frac{\partial v}{\partial T} \right)_P
157+
= \frac{1}{T}
158+
+ \frac{1}{Z} \left( \frac{\partial Z}{\partial T} \right)_P $$
159+
160+
where $P$ is the pressure, $T$ is the temperature, and $Z$ is the
161+
compressibility factor.
162+
163+
Parameters
164+
----------
165+
T : float
166+
Temperature [K].
167+
P : float
168+
Pressure [Pa].
169+
y : FloatVector (N)
170+
Mole fractions of all components [mol/mol].
171+
172+
Returns
173+
-------
174+
float
175+
Thermal expansion coefficient, $\beta$ [K⁻¹].
176+
"""
177+
dZdT, Z = derivative_complex(
178+
lambda x: self.Z(x, P, y),
179+
T,
180+
)
181+
return 1 / T + dZdT / Z
182+
183+
def kappa(
184+
self,
185+
T: float,
186+
P: float,
187+
y: FloatVector,
188+
) -> float:
189+
r"""Calculate the isothermal compressibility coefficient.
190+
191+
$$ \kappa
192+
= - \frac{1}{v} \left( \frac{\partial v}{\partial P} \right)_T
193+
= \frac{1}{P}
194+
- \frac{1}{Z} \left( \frac{\partial Z}{\partial P} \right)_T $$
195+
196+
where $P$ is the pressure, $T$ is the temperature, and $Z$ is the
197+
compressibility factor.
198+
199+
Parameters
200+
----------
201+
T : float
202+
Temperature [K].
203+
P : float
204+
Pressure [Pa].
205+
y : FloatVector (N)
206+
Mole fractions of all components [mol/mol].
207+
208+
Returns
209+
-------
210+
float
211+
Isothermal compressibility coefficient, $\kappa$ [Pa⁻¹].
212+
"""
213+
dZdP, Z = derivative_complex(
214+
lambda x: self.Z(T, x, y),
215+
P,
216+
)
217+
return 1 / P - dZdP / Z
218+
146219
def f(
147220
self,
148221
T: float,
@@ -190,6 +263,86 @@ def Z(self, T: float, P: float, z: FloatVector) -> FloatVector:
190263
"""
191264
pass
192265

266+
def beta(
267+
self,
268+
T: float,
269+
P: float,
270+
z: FloatVector,
271+
) -> FloatVector:
272+
r"""Calculate the thermal expansion coefficients of the possible phases
273+
of a fluid.
274+
275+
$$ \beta
276+
= \frac{1}{v} \left( \frac{\partial v}{\partial T} \right)_P
277+
= \frac{1}{T}
278+
+ \frac{1}{Z} \left( \frac{\partial Z}{\partial T} \right)_P $$
279+
280+
where $P$ is the pressure, $T$ is the temperature, and $Z$ is the
281+
compressibility factor.
282+
283+
Parameters
284+
----------
285+
T : float
286+
Temperature [K].
287+
P : float
288+
Pressure [Pa].
289+
z : FloatVector (N)
290+
Mole fractions of all components [mol/mol].
291+
292+
Returns
293+
-------
294+
FloatVector
295+
Thermal expansion coefficients of the possible phases [K⁻¹].
296+
If two phases are possible, the first result corresponds to the
297+
liquid.
298+
"""
299+
dT = 1.0
300+
Zp = self.Z(T + dT, P, z)
301+
Zm = self.Z(T - dT, P, z)
302+
dZdT = (Zp - Zm) / (2 * dT)
303+
Z = (Zp + Zm) / 2
304+
return 1 / T + dZdT / Z
305+
306+
def kappa(
307+
self,
308+
T: float,
309+
P: float,
310+
z: FloatVector,
311+
) -> FloatVector:
312+
r"""Calculate the isothermal compressibility coefficients of the
313+
possible phases of a fluid.
314+
315+
$$ \kappa
316+
= - \frac{1}{v} \left( \frac{\partial v}{\partial P} \right)_T
317+
= \frac{1}{P}
318+
- \frac{1}{Z} \left( \frac{\partial Z}{\partial P} \right)_T $$
319+
320+
where $P$ is the pressure, $T$ is the temperature, and $Z$ is the
321+
compressibility factor.
322+
323+
Parameters
324+
----------
325+
T : float
326+
Temperature [K].
327+
P : float
328+
Pressure [Pa].
329+
z : FloatVector (N)
330+
Mole fractions of all components [mol/mol].
331+
332+
Returns
333+
-------
334+
float
335+
Isothermal compressibility coefficients of the possible phases [Pa⁻¹].
336+
If two phases are possible, the first result corresponds to the
337+
liquid.
338+
"""
339+
dP = max(P * 1e-2, 1e1)
340+
Zp = self.Z(T, P + dP, z)
341+
Zm = self.Z(T, P - dP, z)
342+
dZdP = (Zp - Zm) / (2 * dP)
343+
Z = (Zp + Zm) / 2
344+
return 1 / P - dZdP / Z
345+
193346
@abstractmethod
194347
def phi(
195348
self,

‎src/polykin/thermo/eos/cubic.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -402,7 +402,7 @@ class RedlichKwong(CubicEoS):
402402
r"""[Redlich-Kwong](https://en.wikipedia.org/wiki/Redlich%E2%80%93Kwong_equation_of_state)
403403
equation of state.
404404
405-
This EoS is based on the following $P(v,T)$ relationship:
405+
This EoS is based on the following $P(T,v)$ relationship:
406406
407407
$$ P = \frac{RT}{v - b_m} -\frac{a_m}{v (v + b_m)} $$
408408
@@ -461,7 +461,7 @@ class SoaveRedlichKwong(CubicEoS):
461461
r"""[Soave-Redlich-Kwong](https://en.wikipedia.org/wiki/Cubic_equations_of_state#Soave_modification_of_Redlich%E2%80%93Kwong)
462462
equation of state.
463463
464-
This EoS is based on the following $P(v,T)$ relationship:
464+
This EoS is based on the following $P(T,v)$ relationship:
465465
466466
$$ P = \frac{RT}{v - b_m} -\frac{a_m}{v (v + b_m)} $$
467467
@@ -537,7 +537,7 @@ class PengRobinson(CubicEoS):
537537
r"""[Peng-Robinson](https://en.wikipedia.org/wiki/Cubic_equations_of_state#Peng%E2%80%93Robinson_equation_of_state)
538538
equation of state.
539539
540-
This EoS is based on the following $P(v,T)$ relationship:
540+
This EoS is based on the following $P(T,v)$ relationship:
541541
542542
$$ P = \frac{RT}{v - b_m} -\frac{a_m}{v^2 + 2 v b_m - b_m^2} $$
543543

‎src/polykin/thermo/eos/idealgas.py‎

Lines changed: 47 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616
class IdealGas(GasEoS):
1717
r"""Ideal gas equation of state.
1818
19-
This EoS is based on the following $P(v,T)$ relationship:
19+
This EoS is based on the following $P(T,v)$ relationship:
2020
2121
$$ P = \frac{R T}{v} $$
2222
@@ -64,6 +64,52 @@ def P(self, T: float, v: float, y=None) -> float:
6464
"""
6565
return R * T / v
6666

67+
def beta(self, T: float, P=None, y=None) -> float:
68+
r"""Calculate the thermal expansion coefficient.
69+
70+
$$ \beta
71+
= \frac{1}{v} \left( \frac{\partial v}{\partial T} \right)_P
72+
= \frac{1}{T} $$
73+
74+
Parameters
75+
----------
76+
T : float
77+
Temperature [K].
78+
P : float
79+
Pressure [Pa].
80+
y : FloatVector (N)
81+
Mole fractions of all components [mol/mol].
82+
83+
Returns
84+
-------
85+
float
86+
Thermal expansion coefficient, $\beta$ [K⁻¹].
87+
"""
88+
return 1 / T
89+
90+
def kappa(self, T: float, P: float, y=None) -> float:
91+
r"""Calculate the isothermal compressibility coefficient.
92+
93+
$$ \kappa
94+
= - \frac{1}{v} \left( \frac{\partial v}{\partial P} \right)_T
95+
= \frac{1}{P} $$
96+
97+
Parameters
98+
----------
99+
T : float
100+
Temperature [K].
101+
P : float
102+
Pressure [Pa].
103+
y : FloatVector (N)
104+
Mole fractions of all components [mol/mol].
105+
106+
Returns
107+
-------
108+
float
109+
Isothermal compressibility coefficient, $\kappa$ [Pa⁻¹].
110+
"""
111+
return 1 / P
112+
67113
def phi(self, T=None, P=None, y=None) -> FloatVector:
68114
r"""Calculate the fugacity coefficients of all components.
69115

‎src/polykin/thermo/eos/virial.py‎

Lines changed: 41 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@
2424
class Virial(GasEoS):
2525
r"""Virial equation of state truncated to the second coefficient.
2626
27-
This EoS is based on the following $P(v,T)$ relationship:
27+
This EoS is based on the following $P(T,v)$ relationship:
2828
2929
$$ P = \frac{R T}{v - B_m} $$
3030
@@ -129,6 +129,31 @@ def Bij(
129129
"""
130130
return B_mixture(T, self.Tc, self.Pc, self.Zc, self.w)
131131

132+
def P(
133+
self,
134+
T: float,
135+
v: float,
136+
y: FloatVector,
137+
) -> float:
138+
r"""Calculate the pressure of the fluid.
139+
140+
Parameters
141+
----------
142+
T : float
143+
Temperature [K].
144+
v : float
145+
Molar volume [m³/mol].
146+
y : FloatVector (N)
147+
Mole fractions of all components [mol/mol].
148+
149+
Returns
150+
-------
151+
float
152+
Pressure [Pa].
153+
"""
154+
Bm = self.Bm(T, y)
155+
return R * T / (v - Bm)
156+
132157
def Z(
133158
self,
134159
T: float,
@@ -164,30 +189,37 @@ def Z(
164189
Bm = self.Bm(T, y)
165190
return 1.0 + Bm * P / (R * T)
166191

167-
def P(
192+
def kappa(
168193
self,
169194
T: float,
170-
v: float,
195+
P: float,
171196
y: FloatVector,
172197
) -> float:
173-
r"""Calculate the pressure of the fluid.
198+
r"""Calculate the isothermal compressibility coefficient.
199+
200+
$$ \kappa
201+
= - \frac{1}{v} \left( \frac{\partial v}{\partial P} \right)_T
202+
= \frac{1}{P Z} $$
203+
204+
where $P$ is the pressure, $T$ is the temperature, and $Z$ is the
205+
compressibility factor.
174206
175207
Parameters
176208
----------
177209
T : float
178210
Temperature [K].
179-
v : float
180-
Molar volume [m³/mol].
211+
P : float
212+
Pressure [Pa].
181213
y : FloatVector (N)
182214
Mole fractions of all components [mol/mol].
183215
184216
Returns
185217
-------
186218
float
187-
Pressure [Pa].
219+
Isothermal compressibility coefficient, $\kappa$ [Pa⁻¹].
188220
"""
189-
Bm = self.Bm(T, y)
190-
return R * T / (v - Bm)
221+
Z = self.Z(T, P, y)
222+
return 1 / (P * Z)
191223

192224
def phi(
193225
self,

0 commit comments

Comments
 (0)