-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathkempers01.py
283 lines (236 loc) · 11.8 KB
/
kempers01.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
'''
Author: Vegard Gjeldvik Jervell
Date: December 2021
Purpose: Implementation of the model proposed by Kempers (J. Chem. Phys. 115, 6330, 2001) for prediction of the Soret coefficient
doi : http://dx.doi.org/10.1063/1.1398315
Requires: numpy, scipy, ThermoPack (https://github.com/SINTEF/thermopack), KineticGas (https://github.com/vegardjervell/Kineticgas)
'''
import warnings, platform
import numpy as np
from scipy.constants import gas_constant, Avogadro, Boltzmann
from scipy.optimize import root
from pyctp import cubic
from pykingas import KineticGas
class Kempers(object):
def __init__(self, comps, eos, alpha_t0_N = 7, sigma=None, eps_div_k=None, mole_weights=None):
'''
:param comps (str): comma separated list of components
:param eos (ThermoPack): Initialized Equation of State object, initialized with components 'comp'
:param alpha_t0_N (int, optional): Order of approximation of Enskog solutions
:param sigma (array, optional): Hard sphere diameters [m]
:param eps_div_k (array, optional): Interaction potential well depth [k_B]
:param mole_weights (array, optional): Molar masses [g / mol]
'''
self.comps = comps
self.total_moles = 1 #Dummy value, because some ThermoPack methods require it. Does not effect output.
eoscomp_inds = [eos.getcompindex(comp) for comp in self.comps.split(',')]
if -1 in eoscomp_inds:
warnings.warn('Equation of state and model must be initialized with same components.\n'
"I'm initializing using SRK with "+comps+" now to avoid crashing")
self.eos = cubic.cubic()
self.eos.init(self.comps, 'SRK')
elif any(np.array(eoscomp_inds) != sorted(np.array(eoscomp_inds))):
eoscomps = ','.join(eos.get_comp_name(i) for i in sorted(eoscomp_inds))
warnings.warn('Equation of state and model must be initialized with same components in the same order\n'
'but are initialized with ' + eoscomps + ' and ' + self.comps+'.\n'
"I'm initializing using SRK with "+comps+" now to avoid crashing")
self.eos = cubic.cubic()
self.eos.init(self.comps, 'SRK')
else:
self.eos = eos
if mole_weights is None:
self.mole_weights = np.array([self.eos.compmoleweight(self.eos.getcompindex(comp)) for comp in comps.split(',')])
else:
self.mole_weights = np.copy(mole_weights)
self.kinetic_gas = KineticGas(comps, mole_weights=self.mole_weights, sigma=sigma, eps_div_k=eps_div_k)
self.alpha_t0_N = alpha_t0_N
def get_binary_soret_cov(self, T, p, x, phase, kin=False, BH=False):
'''
Get soret coefficients at specified phase point, for binary mixture, center of volume frame of reference
:param T (float): Temperature [K]
:param p (float): Pressure [Pa]
:param x (ndarray): Composition [mole fraction]
:param phase (int): ThermoPack phase key
:param kin (bool, optional): Return kinetic gas value?
:param BH (bool, optional): Use Barker-Henderson diameters?
:return: (ndarray) soret coefficients
'''
R = gas_constant
v, dvdn = self.eos.specific_volume(T, p, x, phase, dvdn=True)
h, dhdn = self.eos.enthalpy(T, p, x, phase, dhdn=True)
h0, dh0dn = self.eos.enthalpy(T, 1e-5, x, 2, dhdn=True)
dmudx = self.dmudx_TP(T, p, x, phase)
alpha_T0 = self.kinetic_gas.alpha_T0(T, v, x, BH=BH, N=self.alpha_t0_N)
v1, v2 = dvdn
h1, h2 = dhdn
h10, h20 = dh0dn
x1, x2 = x
dmu1dx1 = dmudx[0, 0]
alpha_1 = (v1 * v2 / (v1 * x1 + v2 * x2)) * (((h2 - h20) / v2) - ((h1 - h10) / v1)) / (x1 * dmu1dx1) + (R * T * alpha_T0[0]) / (x1 * dmu1dx1)
alpha = np.array([alpha_1, - alpha_1])
soret = alpha / T
if kin is True:
kin_contrib = alpha_T0 / T # * R /(self.x * dmudx.diagonal())
return soret, kin_contrib
else:
return soret
def get_soret_cov(self, T, p, x, phase, kin=False, BH=False):
'''
Get soret coefficients at specified phase point, center of volume frame of reference
:param T (float): Temperature [K]
:param p (float): Pressure [Pa]
:param x (ndarray): Composition [mole fraction]
:param phase (int): ThermoPack phase key
:param kin (bool, optional): Return kinetic gas value?
:param BH (bool, optional): Use Barker-Henderson diameters?
:return: (ndarray) soret coefficients
'''
if len(x) == 2:
return self.get_binary_soret_cov(T, p, x, phase, kin=kin)
R = gas_constant
v, dvdn = self.eos.specific_volume(T, p, x, phase, dvdn=True)
h, dhdn = self.eos.enthalpy(T, p, x, phase, dhdn=True)
h0, dh0dn = self.eos.enthalpy(T, 1e-5, x, 2, dhdn=True)
dmudx = self.dmudx_TP(T, p, x, phase)
alpha_T0 = self.kinetic_gas.alpha_T0(T, v, x, BH=BH, N=self.alpha_t0_N)
#using alpha_T0 as initial guess for root solver
initial_guess = alpha_T0
N = len(self.x)
#Defining the set of equations
def eq_set(alpha):
eqs = np.zeros(N)
for i in range(N-1):
eqs[i] = ((dhdn[-1] - dh0dn[-1])/dvdn[-1]) - ((dhdn[i] - dh0dn[i])/dvdn[i])\
+ R * T * ((alpha_T0[i] * (1 - x[i]) / dvdn[i])
- (alpha_T0[-1] * (1 - x[-1])/dvdn[-1]))\
- sum((dmudx[i, j]/dvdn[i] - dmudx[-1, j]/dvdn[-1]) * x[j] * (1 - x[j]) * alpha[j]
for j in range(N - 1))
eqs[N-1] = sum(x * (1 - x) * alpha)
return eqs
#Solve the set of equations, warn if non-convergent
solved = root(eq_set, initial_guess)
if solved.success is False:
warnings.warn('Solution did not converge for composition :' + str(x) + ', Temperature :'+ str(T) + ', Pressure : '+str(p/1e5))
alpha = np.full_like(solved.x, np.nan)
alpha = solved.x
else:
alpha = solved.x
soret = alpha / T
if kin is True:
kin_contrib = alpha_T0 / T #* R /(self.x * dmudx.diagonal())
return soret, kin_contrib
else:
return soret
def get_binary_soret_com(self, T, p, x, phase, kin=False, BH=False):
'''
Get soret coefficients at specified phase point, for binary mixture, center of mass frame of reference
:param T (float): Temperature [K]
:param p (float): Pressure [Pa]
:param x (ndarray): Composition [mole fraction]
:param phase (int): ThermoPack phase key
:param kin (bool, optional): Return kinetic gas value?
:param BH (bool, optional): Use Barker-Henderson diameters?
:return: (ndarray) soret coefficients
'''
R = gas_constant
v, = self.eos.specific_volume(T, p, x, phase)
h, dhdn = self.eos.enthalpy(T, p, x, phase, dhdn=True)
h0, dh0dn = self.eos.enthalpy(T, 1e-5, x, 2, dhdn=True)
dmudx = self.dmudx_TP(T, p, x, phase)
alpha_T0 = self.kinetic_gas.alpha_T0(T, v, x, BH=BH, N=self.alpha_t0_N)
m1, m2 = self.mole_weights
h1, h2 = dhdn
h10, h20 = dh0dn
x1, x2 = x
dmu1dx1 = dmudx[0, 0]
alpha_1 = (m1 * m2 / (m1 * x1 + m2 * x2)) * (((h2 - h20) / m2) - ((h1 - h10) / m1)) / (x1 * dmu1dx1) + (
R * T * alpha_T0[0]) / (x1 * dmu1dx1)
alpha = np.array([alpha_1, - alpha_1])
soret = alpha / T
if kin is True:
kin_contrib = alpha_T0 / T # * R /(self.x * dmudx.diagonal())
return soret, kin_contrib
else:
return soret
def get_soret_com(self, T, p, x, phase, kin=False, BH=False):
'''
Get soret coefficients at specified phase point, center of mass frame of reference
:param T (float): Temperature [K]
:param p (float): Pressure [Pa]
:param x (ndarray): Composition [mole fraction]
:param phase (int): ThermoPack phase key
:param kin (bool, optional): Return kinetic gas value?
:param BH (bool, optional): Use Barker-Henderson diameters?
:return: (ndarray) soret coefficients
'''
if len(x) == 2:
return self.get_binary_soret_com(T, p, x, phase, kin=kin, BH=BH)
R = gas_constant
M = self.mole_weights
v, = self.eos.specific_volume(T, p, x, phase)
h, dhdn = self.eos.enthalpy(T, p, x, phase, dhdn=True)
h0, dh0dn = self.eos.enthalpy(T, 1e-5, x, 2, dhdn=True)
dmudx = self.dmudx_TP(T, p, x, phase)
alpha_T0 = self.kinetic_gas.alpha_T0(T, v, x, BH=BH, N=self.alpha_t0_N)
# using alpha_T0 as initial guess for root solver
initial_guess = alpha_T0
N = len(self.x)
# Defining the set of equations
def eq_set(alpha):
eqs = np.zeros(N)
for i in range(N - 1):
eqs[i] = ((dhdn[-1] - dh0dn[-1]) / M[-1]) - ((dhdn[i] - dh0dn[i]) / M[i]) \
+ R * T * ((alpha_T0[i] * (1 - x[i]) / M[i])
- (alpha_T0[-1] * (1 - x[-1]) / M[-1])) \
- sum((dmudx[i, j] / M[i] - dmudx[-1, j] / M[-1]) * x[j] * (1 - x[j]) * alpha[j]
for j in range(N - 1))
eqs[N - 1] = sum(x * (1 - x) * alpha)
return eqs
# Solve the set of equations, warn if non-convergent
solved = root(eq_set, initial_guess)
if solved.success is False:
warnings.warn('Solution did not converge for composition :'+str(x)+ ', Temperature :' + str(T))
alpha = np.full_like(solved.x, np.nan)
else:
alpha = solved.x
soret = alpha / T
if kin is True:
kin_contrib = alpha_T0 / T #* R /(self.x * dmudx.diagonal())
return soret, kin_contrib
else:
return soret
def get_soret(self, T, p, x, phase, mode, kin=False, BH=False):
'''
Get soret coefficients at specified phase point
:param T (float): Temperature [K]
:param p (float): Pressure [Pa]
:param x (ndarray): Composition [mole fraction]
:param phase (int): ThermoPack phase key
:param mode (str): 'com' or 'cov' to determine mode
:param kin (bool, optional): Return kinetic gas value?
:param BH (bool, optional): Use Barker-Henderson diameters?
:return: (ndarray) Soret coefficients
'''
if mode == 'cov':
return self.get_soret_cov(T, p, x, phase, kin=kin, BH=BH)
elif mode == 'com':
return self.get_soret_com(T, p, x, phase, kin=kin, BH=BH)
def dmudn_TP(self, T, p, x, phase):
'''
Calculate chemical potential derivative with respect to number of moles at constant temperature and pressure
:return: ndarray, dmudn[i,j] = dmu_i / dn_j
'''
v, dvdn = self.eos.specific_volume(T, p, x, phase, dvdn=True)
mu, dmudn_TV = self.eos.chemical_potential_tv(T, v * self.total_moles,
x * self.total_moles, dmudn=True)
pres, dpdn = self.eos.pressure_tv(T, v * self.total_moles, x * self.total_moles, dpdn=True)
return dmudn_TV - np.tensordot(dpdn, dvdn, axes=0)
def dmudx_TP(self, T, p, x, phase):
'''
Calculate chemical potential derivative with respect to mole fraction of components
at constant temperature and pressure
:return: ndarray, dmudx[i,j] = dmu_i / dx_j
'''
x = np.array(x)
dmudn = self.dmudn_TP(T, p, x, phase)
return dmudn * self.total_moles