OpenLB 1.7
Loading...
Searching...
No Matches
interactionPotential.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2015 Peter Weisbrod
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22*/
23
24#ifndef INTERACTION_POTENTIAL_HH
25#define INTERACTION_POTENTIAL_HH
26
27
29
30
31namespace olb {
32
33
34
35template <typename T, typename S>
36ShanChen93<T,S>::ShanChen93(T rhoZero) : AnalyticalF<1,T,S>(1), _rhoZero(rhoZero)
37{
38 this->getName() = "ShanChen93";
39}
40
41template <typename T, typename S>
42bool ShanChen93<T,S>::operator()(T psi[], const S rho[])
43{
44 psi[0]=util::sqrt(_rhoZero)*(1-util::exp(-(rho[0]/_rhoZero)));
45 return true;
46}
47
48
49template <typename T, typename S>
50ShanChen94<T,S>::ShanChen94(T rhoZero, T psiZero) : AnalyticalF<1,T,S>(1), _rhoZero(rhoZero), _psiZero(psiZero)
51{
52 this->getName() = "ShanChen94";
53}
54
55template <typename T, typename S>
56bool ShanChen94<T,S>::operator()(T psi[], const S rho[])
57{
58 psi[0]=_psiZero*util::exp(-_rhoZero/rho[0]);
59 return true;
60}
61
62
63template <typename T, typename S>
64PengRobinson<T,S>::PengRobinson(T G, T acentricFactor, T a, T b, T tr) : AnalyticalF<1,T,S>(1), _G(G), _acentricFactor(acentricFactor), _a(a), _b(b)
65{
66 _R = 1.;
67 //a=0.45724*R*R*tc*tc/pc;
68 //b=0.0778*R*tc/pc;
69 _tc = 0.0778/0.45724*_a/_b/_R;
70 //T pc = 0.0778*_R*tc/_b;
71 //T rhoc = pc/0.307/_R/tc;
72 _t = _tc*tr;
73 //Zc=0.307 Tc=0.072922004 pc=0.059569985 rhoc=2.6609121
74 _alpha = 1. + (0.37464+1.54226*_acentricFactor-0.26992*_acentricFactor*_acentricFactor)*(1.-util::sqrt(_t/_tc));
75 _alpha = _alpha*_alpha;
76 this->getName() = "PengRobinson";
77}
78
79template <typename T, typename S>
80bool PengRobinson<T,S>::operator()(T psi[], const S rho[])
81{
82 T p = (rho[0]*_R*_t/(1.-_b*rho[0]))-(_a*_alpha*rho[0]*rho[0]/(1.+2.*_b*rho[0]-_b*_b*rho[0]*rho[0]));
83 psi[0] = util::sqrt(6.*(p-rho[0]/3.)/_G);
84 return true;
85}
86
87// second operator allows to incorporate temperature changes
88template <typename T, typename S>
89bool PengRobinson<T,S>::operator()(T psi[], const S rho[], const S t[])
90{
91 _t = t[0];
92 _alpha = 1. + (0.37464+1.54226*_acentricFactor-0.26992*_acentricFactor*_acentricFactor)*(1.-util::sqrt(_t/_tc));
93 _alpha = _alpha*_alpha;
94 T p = (rho[0]*_R*_t/(1.-_b*rho[0]))-(_a*_alpha*rho[0]*rho[0]/(1.+2.*_b*rho[0]-_b*_b*rho[0]*rho[0]));
95 psi[0] = util::sqrt(6.*(p-rho[0]/3.)/_G);
96 return true;
97}
98
99
100template <typename T, typename S>
101CarnahanStarling<T,S>::CarnahanStarling(T G, T a, T b, T tr) : AnalyticalF<1,T,S>(1), _G(G), _a(a), _b(b)
102{
103 _R = 1.;
104 //a=0.4963*tc*tc*R*R/pc;
105 //b=0.18727*R*tc/pc;
106 T tc = 0.18727/0.4963*_a/_b/_R;
107 //T pc = 0.18727*_R*tc/_b;
108 //T rhoc = pc/0.35930763/_R/tc;
109 _t = tc*tr;
110 //Zc=0.35930763 Tc=0.094333065 pc=0.0044164383 rhoc=0.13029921
111 this->getName() = "CarnahanStarling";
112}
113
114template <typename T, typename S>
115bool CarnahanStarling<T,S>::operator()(T psi[], const S rho[])
116{
117 T c = _b*rho[0]/4.;
118 T p = rho[0]*_R*_t*((1.+c+c*c-c*c*c)/(1.-c)/(1.-c)/(1.-c))-_a*rho[0]*rho[0];
119 psi[0] = util::sqrt(6.*(p-rho[0]/3.)/_G);
120 return true;
121}
122
123// second operator allows to incorporate temperature changes
124template <typename T, typename S>
125bool CarnahanStarling<T,S>::operator()(T psi[], const S rho[], const S t[])
126{
127 _t = t[0];
128 T c = _b*rho[0]/4.;
129 T p = rho[0]*_R*_t*((1.+c+c*c-c*c*c)/(1.-c)/(1.-c)/(1.-c))-_a*rho[0]*rho[0];
130 psi[0] = util::sqrt(6.*(p-rho[0]/3.)/_G);
131 return true;
132}
133
134
135template <typename T, typename S>
137{
138 this->getName() = "PsiEqualsRho";
139}
140
141template <typename T, typename S>
142bool PsiEqualsRho<T,S>::operator()(T psi[], const S rho[])
143{
144 psi[0]=rho[0];
145 return true;
146}
147
148
149template <typename T, typename S>
150Krause<T,S>::Krause(T rhoZero, T psiZero) : AnalyticalF<1,T,S>(1), _rhoZero(rhoZero), _psiZero(psiZero)
151{
152 this->getName() = "Krause";
153}
154
155template <typename T, typename S>
156bool Krause<T,S>::operator()(T psi[], const S rho[])
157{
158 psi[0]=_psiZero/1.77*1.414/_rhoZero*util::exp(-(_rhoZero-rho[0])*(_rhoZero-rho[0])/_rhoZero/_rhoZero);
159 return true;
160}
161
162
163template <typename T, typename S>
164WeisbrodKrause<T,S>::WeisbrodKrause(T rhoZero, T sigmu) : AnalyticalF<1,T,S>(1), _rhoZero(rhoZero), _sigmu(sigmu)
165{
166 _rhoZero=_rhoZero/1.005088;
167 this->getName() = "WeisbrodKrause";
168}
169
170template <typename T, typename S>
171bool WeisbrodKrause<T,S>::operator()(T psi[], const S rho[])
172{
173 psi[0]=util::sqrt(_rhoZero)*1.5179/_sigmu*util::exp(-(_sigmu-rho[0]/_rhoZero)*(_sigmu-rho[0]/_rhoZero)/_sigmu/_sigmu);
174 return true;
175}
176
177
178template <typename T, typename S>
179Normal<T,S>::Normal(T sigma, T mu) : AnalyticalF<1,T,S>(1), _sigma(sigma), _mu(mu)
180{
181 this->getName() = "Normal";
182}
183
184template <typename T, typename S>
185bool Normal<T,S>::operator()(T psi[], const S rho[])
186{
187 psi[0]=1./2.507/_sigma*util::exp(-(rho[0]-_mu)*(rho[0]-_mu)/_sigma/_sigma/2.);
188 return true;
189}
190
191MultiComponentPengRobinson::MultiComponentPengRobinson(double p_, double T_, std::vector<double> z_, std::vector<double> a_L, std::vector<double> b_L,
192 std::vector<double> M_L, std::vector<double> Tc_L, std::vector<double> pc_L, std::vector<double> omega_, std::vector<double> devi,
193 std::vector<double> alpha_, std::vector<double> gI_L, std::vector<double> gII_L) :
194p(p_), T(T_), z(z_), num_comp(z.size()), a_c(a_L), b(b_L), M(M_L), T_c(Tc_L), p_c(pc_L), omega(omega_), m(devi), alpha(alpha_), g_I(gI_L), g_II(gII_L)
195{ }
196double MultiComponentPengRobinson::G_Excess(std::vector<double> x_i, std::vector<double> g_jk, std::vector<double> beta_jk){
197 double G_E = 0;
198 for(int i = 0; i < num_comp; i++){
199 double sum1 = 0;
200 double sum2 = 0;
201 for(int j = 0; j < num_comp; j++){
202 sum1 = sum1 + x_i[j]*g_jk[i+num_comp*j]*beta_jk[i+num_comp*j];
203 sum2 = sum2 + x_i[j]*beta_jk[i+num_comp*j];
204 }
205 G_E = G_E + x_i[i] * sum1/sum2;
206 }
207 return G_E;
208}
209double MultiComponentPengRobinson::b_mix(std::vector<double> x_i){
210 double b_m = 0;
211 for(int i = 0; i < num_comp; i++){
212 b_m = b_m + b[i]*x_i[i];
213 }
214 return b_m;
215}
216double MultiComponentPengRobinson::a_mix(std::vector<double> x_i, std::vector<double> a_i, double b_m, double G_E){
217 double a_sum = 0;
218 for(int i = 0; i < num_comp; i++){
219 a_sum = a_sum + x_i[i]*a_i[i]/b[i];
220 }
221 /*double a_m = 0;
222 for(int i = 0; i < num_comp; i++){
223 a_i[i] = a_c[i] * util::pow((1 + m[i] * (1 - util::pow((_T/T_c[i]),0.5))),2);
224 for(int j = 0; j < num_comp; j++){
225 a_m += x_i[i]*x_i[j]*util::pow(a_i[i]*a_i[j],0.5);
226 }
227 }
228 return a_m;*/
229 return b_m * (a_sum - G_E / lambda);
230}
231double MultiComponentPengRobinson::p_PR(double v, double a_m, double b_m){
232 return R*T/(v - b_m) - a_m/(util::pow(v,2) + 2*v*b_m - util::pow(b_m,2));
233}
234double MultiComponentPengRobinson::gamma_i(std::vector<double> x_i, double a_i, std::vector<double> g_jk, std::vector<double> beta_jk, double _T, int i){
235 double sum1 = 0;
236 double sum2 = 0;
237 double sum3 = 0;
238 double sum4 = 0;
239 for(int j = 0; j < num_comp; j++){
240 sum1 = sum1 + x_i[j]*g_jk[num_comp*j+i]*beta_jk[num_comp*j+i];
241 sum2 = sum2 + x_i[j]*beta_jk[num_comp*j+i];
242 double sub31 = 0;
243 double sub41 = 0;
244 double sub42 = 0;
245 for(int k = 0; k < num_comp; k++){
246 sub31 = sub31 + x_i[k]*beta_jk[j+num_comp*k];
247 sub41 = sub41 + x_i[k]*g_jk[j+num_comp*k]*beta_jk[j+num_comp*k];
248 sub42 = sub42 + x_i[k]*beta_jk[j+num_comp*k];
249 }
250 sum3 = sum3 + x_i[j]*g_jk[num_comp*i+j]*beta_jk[num_comp*i+j] / sub31;
251 sum4 = sum4 + x_i[j]*beta_jk[num_comp*i+j]*sub41/util::pow(sub42,2);
252 }
253 return 1/(R*_T) * (a_i / b[i] - 1/lambda * (sum1/sum2 + sum3 - sum4));
254}
255double MultiComponentPengRobinson::lnfugacity_i(std::vector<double> x_i, double v, double a_m, double b_m, double gamma_i, double _T, int i){
256 //TODO change inputs for non-isothermal cases!!!
257 double term1 = b[i]/(v-b_m) + util::log(x_i[i]*R*_T/(v-b_m));
258 double term2 = -gamma_i / util::pow(2,1.5) * util::log((v + (util::pow(2,0.5)+1)*b_m)/(v-(util::pow(2,0.5)-1)*b_m));
259 double term3 = -b[i]*a_m/(util::pow(2,1.5)*b_m*R*_T)*
260 ((util::pow(2,0.5)+1)/(v+(util::pow(2,0.5)+1)*b_m)+(util::pow(2,0.5)-1)/(v-(util::pow(2,0.5)-1)*b_m));
261 /*std::vector<double> a_i(num_comp);
262 a_m = 0;
263 for(int i = 0; i < num_comp; i++){
264 a_i[i] = a_c[i] * util::pow((1 + m[i] * (1 - util::pow((_T/T_c[i]),0.5))),2);
265 for(int j = 0; j < num_comp; j++){
266 a_m += x_i[i]*x_i[j]*util::pow(a_i[i]*a_i[j],0.5);
267 }
268 }
269 gamma_i = -a_m/(b_m*R*_T)*(b[i]/b_m-2/a_m*(x_i[0]*util::pow(a_i[0]*a_i[i],0.5)+x_i[1]*util::pow(a_i[1]*a_i[i],0.5)));
270 term1 = b[i]/(v-b_m) + util::log(x_i[i]*R*_T/(v-b_m));
271 term2 = -gamma_i / util::pow(2,1.5) * util::log((v + (util::pow(2,0.5)+1)*b_m)/(v-(util::pow(2,0.5)-1)*b_m));
272 term3 = -b[i]*a_m/(util::pow(2,1.5)*b_m*R*_T)*((util::pow(2,0.5)+1)/(v+(util::pow(2,0.5)+1)*b_m)+(util::pow(2,0.5)-1)/(v-(util::pow(2,0.5)-1)*b_m));*/
273 return term1 + term2 + term3;
274}
275std::vector<double> MultiComponentPengRobinson::iterate_VLE(double residuum, double beta0){
276 std::vector<double> a_i(num_comp), g_jk(num_comp*num_comp), beta_jk(num_comp*num_comp);
277 for(int i = 0; i < num_comp; i++){
278 a_i[i] = a_c[i] * util::pow((1 + m[i] * (1 - util::pow((T/T_c[i]),0.5))),2);
279 for(int j = 0; j < num_comp; j++){
280 g_jk[num_comp*i+j] = g_I[num_comp*i+j] + g_II[num_comp*i+j] * T;
281 beta_jk[num_comp*i+j] = b[j] * util::exp(-alpha[num_comp*i+j]*g_jk[num_comp*i+j]/(R*T));
282 }
283 }
284 std::vector<double> vx(2+2*num_comp), K(num_comp), x(num_comp), y(num_comp), G_E(2), b_m(2), a_m(2), gammas(2*num_comp);
285 std::vector<double> deltaLogFugacity(num_comp), errorFugacity(num_comp);
286 double totalError = 0;
287 for(int i = 0; i < num_comp; i++){
288 //adjust prefactor or initial values for partition coefficients, 0.9 for hydrocarbons
289 K[i] = 1.0*(p_c[i]/(p))*util::exp(5.373*(1+omega[i])*(1-T_c[i]/T));
290 errorFugacity[i] = 1;
291 }
292 for(int i = 0; i < 2; i++){
293 vx[i] = init_v[i];
294 }
295 double beta = beta0;
296 double f = 0;
297 do{
298 totalError = 0;
299 f = 0;
300 for(int i = 0; i < num_comp; i++){
301 f += (K[i]-1)*z[i]/(1+(K[i]-1)*beta);
302 }
303 while(abs(f) > residuum){
304 double dfdbeta = 0;
305 for(int i = 0; i < num_comp; i++){
306 dfdbeta += -1*(K[i]-1)*(K[i]-1)*z[i]/((1+(K[i]-1)*beta)*(1+(K[i]-1)*beta));
307 }
308 beta = beta - f / dfdbeta;
309 f = 0;
310 for(int i = 0; i < num_comp; i++){
311 f += (K[i]-1)*z[i]/(1+(K[i]-1)*beta);
312 }
313 }
314 double x_sum = 0;
315 double y_sum = 0;
316 for(int i = 0; i < num_comp; i++){
317 x[i] = z[i]/(1+(K[i]-1)*beta);
318 y[i] = K[i]*x[i];
319 x_sum += x[i];
320 y_sum += y[i];
321 }
322 if(x_sum != 1 || y_sum != 1){
323 for(int i = 0; i < num_comp; i++){
324 vx[i+2] = x[i]/x_sum;
325 x[i] = vx[i+2];
326 vx[i+2+num_comp] = y[i]/y_sum;
327 y[i] = vx[i+2+num_comp];
328 }
329 }
330
331 G_E[0] = G_Excess(x, g_jk, beta_jk);
332 G_E[1] = G_Excess(y, g_jk, beta_jk);
333 b_m[0] = b_mix(x);
334 b_m[1] = b_mix(y);
335 a_m[0] = a_mix(x, a_i, b_m[0], G_E[0]);
336 a_m[1] = a_mix(y, a_i, b_m[1], G_E[1]);
337 for(int i = 0; i < num_comp; i++){
338 gammas[i] = gamma_i(x, a_i[i], g_jk, beta_jk, T, i);
339 gammas[i+num_comp] = gamma_i(y, a_i[i], g_jk, beta_jk, T, i);
340 }
341
342 for(int i = 0; i < 2; i++){
343 while(abs(p-p_PR(vx[i], a_m[i], b_m[i])) > residuum){
344 double dpdv = -R*T/util::pow(vx[i]-b_m[i],2) +
345 2*a_m[i]*(vx[i]+b_m[i])/util::pow(util::pow(vx[i],2) + 2*vx[i]*b_m[i] - util::pow(b_m[i],2),2);
346 vx[i] = vx[i] - (p - p_PR(vx[i], a_m[i], b_m[i])) / (-dpdv);
347 }
348 }
349 for(int i = 0; i < num_comp; i++){
350 deltaLogFugacity[i] = lnfugacity_i(x, vx[0], a_m[0], b_m[0], gammas[i], T, i) -
351 lnfugacity_i(y, vx[1], a_m[1], b_m[1], gammas[i+num_comp], T, i);
352 errorFugacity[i] = util::exp(deltaLogFugacity[i]) - 1;
353 K[i] = util::exp(deltaLogFugacity[i]) / (x[i]/y[i]);
354 totalError += util::pow(errorFugacity[i],2);
355 }
356 }while(util::pow(totalError,0.5) > residuum);
357 for(int i = 1; i<num_comp; i++){
358 volatilities[i] = (M[i]*(vx[i+2]/vx[0] - vx[i+2+num_comp]/vx[1])) / (M[0]*(vx[2]/vx[0] - vx[2+num_comp]/vx[1]));
359 }
360 double vol_sum = 0;
361 for(int i = 0; i<num_comp; i++){
362 vol_sum += volatilities[i];
363 }
364 for(int i = 0; i<num_comp; i++){
365 chis[i] = volatilities[i] / vol_sum;
366 }
367 return vx;
368}
369std::vector<double> MultiComponentPengRobinson::getChis(int n){
370 std::vector<double> n_chis(n);
371 for(int i = 0; i<n; i++){
372 n_chis[i] = chis[i];
373 }
374 return n_chis;
375}
376
377} // end namespace olb
378
379#endif
AnalyticalF are applications from DD to XD, where X is set by the constructor.
bool operator()(T psi[], const S rho[]) override
has to be implemented for 'every' derived class
CarnahanStarling(T G, T a=1., T b=4., T tr=.7)
std::string & getName()
read and write access to name
Definition genericF.hh:51
bool operator()(T psi[], const S rho[])
has to be implemented for 'every' derived class
Krause(T rhoZero=1., T psiZero=1.9)
std::vector< double > iterate_VLE(double residuum, double beta0)
double a_mix(std::vector< double > x_i, std::vector< double > a_i, double b, double G_E)
double b_mix(std::vector< double > x_i)
double p_PR(double v, double a_m, double b_m)
double G_Excess(std::vector< double > x_i, std::vector< double > g_jk, std::vector< double > beta_jk)
MultiComponentPengRobinson(double p_, double T_, std::vector< double > z_, std::vector< double > a_L, std::vector< double > b_L, std::vector< double > M_L, std::vector< double > Tc_L, std::vector< double > pc_L, std::vector< double > omega_, std::vector< double > devi, std::vector< double > alpha_, std::vector< double > gI_L, std::vector< double > gII_L)
std::vector< double > getChis(int n)
double lnfugacity_i(std::vector< double > x_i, double v, double a_m, double b_m, double gamma_i, double _T, int i)
double gamma_i(std::vector< double > x_i, double a_i, std::vector< double > g_jk, std::vector< double > beta_jk, double _T, int i)
bool operator()(T psi[], const S rho[])
has to be implemented for 'every' derived class
Normal(T sigma=1., T mu=1.)
bool operator()(T psi[], const S rho[])
has to be implemented for 'every' derived class
PengRobinson(T G, T acentricFactor=0.334, T a=2./49., T b=2./21., T tr=.8)
bool operator()(T psi[], const S rho[]) override
has to be implemented for 'every' derived class
bool operator()(T psi[], const S rho[])
has to be implemented for 'every' derived class
ShanChen94(T rhoZero=200., T psiZero=4.)
bool operator()(T psi[], const S rho[]) override
has to be implemented for 'every' derived class
bool operator()(T psi[], const S rho[])
has to be implemented for 'every' derived class
WeisbrodKrause(T rhoZero=1., T sigmu=1.)
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
ADf< T, DIM > log(const ADf< T, DIM > &a)
Definition aDiff.h:475
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
ADf< T, DIM > exp(const ADf< T, DIM > &a)
Definition aDiff.h:455
Top level namespace for all of OpenLB.
std::enable_if_t< std::is_arithmetic< T >::type::value, T > abs(T x)
Definition util.h:396