OpenLB 1.7
Loading...
Searching...
No Matches
interactionPotential.h
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2024 Tim Bingert, Luiz Czelusniak,
4 * Maximilian Schecher, Adrian Kummerlaender
5 * E-mail contact: info@openlb.net
6 * The most recent release of OpenLB can be downloaded at
7 * <http://www.openlb.net/>
8 *
9 * This program is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU General Public License
11 * as published by the Free Software Foundation; either version 2
12 * of the License, or (at your option) any later version.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public
20 * License along with this program; if not, write to the Free
21 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
22 * Boston, MA 02110-1301, USA.
23*/
24
25#ifndef INTERACTION_POTENTIAL_H
26#define INTERACTION_POTENTIAL_H
27
28
35namespace olb {
36
37//Calculate VLEs with PR EoS with non-classical Huron-Vidal mixing rule for n-component mixtures
39public:
40 MultiComponentPengRobinson(double p_, double T_, std::vector<double> z_, std::vector<double> a_L, std::vector<double> b_L,
41 std::vector<double> M_L, std::vector<double> Tc_L, std::vector<double> pc_L, std::vector<double> omega_, std::vector<double> devi,
42 std::vector<double> alpha_, std::vector<double> gI_L, std::vector<double> gII_L);
43 double G_Excess(std::vector<double> x_i, std::vector<double> g_jk, std::vector<double> beta_jk);
44 double b_mix(std::vector<double> x_i);
45 double a_mix(std::vector<double> x_i, std::vector<double> a_i, double b, double G_E);
46 double p_PR(double v, double a_m, double b_m);
47 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);
48 double lnfugacity_i(std::vector<double> x_i, double v, double a_m, double b_m, double gamma_i, double _T, int i);
49 std::vector<double> iterate_VLE(double residuum, double beta0);
50 std::vector<double> getChis(int n);
51private:
52 double p, T;
53 std::vector<double> z;
54 const double G{-1.};
55 const double R{1.};
56 const double lambda{0.6232252401};
57 int num_comp;
58 std::vector<double> a_c;
59 std::vector<double> b;
60 std::vector<double> M;
61 std::vector<double> T_c;
62 std::vector<double> p_c;
63 std::vector<double> omega;
64 std::vector<double> m;
65 std::vector<double> alpha;
66 std::vector<double> g_I;
67 std::vector<double> g_II;
68 std::vector<double> init_v = {0.105,120.0};
69 std::vector<double> volatilities = {1.,0.,0.,0.};
70 std::vector<double> chis = {0.5,0.5,0.5,0.5};
71};
72
73namespace interaction {
74
75// established -- only multicomponent flow
78
79 template <typename V, typename PARAMETERS>
80 V compute(V rho, PARAMETERS& params) any_platform {
81 return rho;
82 };
83};
84
85struct ShanChen94 {
87
88 template <typename V, typename PARAMETERS>
89 V compute(V rho, PARAMETERS& params) any_platform {
90 return 4 * util::exp(-200 / rho);
91 };
92};
93
95 struct G : public descriptors::FIELD_BASE<1> { };
96 struct A : public descriptors::FIELD_BASE<1> { };
97 struct B : public descriptors::FIELD_BASE<1> { };
98 struct T : public descriptors::FIELD_BASE<1> { };
99
101
102 template <typename V, typename PARAMETERS>
103 V compute(V rho, PARAMETERS& params) any_platform {
104 V c = params.template get<B>() * rho / V{4};
105 V p = rho
106 * (V(0.18727 / 0.4963) * params.template get<A>() / params.template get<B>())
107 * params.template get<T>()
108 * ((V{1} +c+c*c-c*c*c) / (V{1}-c) / (V{1}-c) / (V{1}-c))
109 - params.template get<A>() * rho*rho;
110 return util::sqrt(V{6}*(p-rho/V{3})/ params.template get<G>());
111 };
112};
113
114template <unsigned N_COMPONENTS>
116 template <typename X_I, typename G_JK, typename BETA_JK, typename V=typename X_I::value_t>
117 V G_Excess(const X_I& x_i, const G_JK& g_jk, const BETA_JK& beta_jk) any_platform {
118 V G_E = 0;
119 for(unsigned i = 0; i < N_COMPONENTS; i++){
120 V sum1 = 0;
121 V sum2 = 0;
122 for(unsigned j = 0; j < N_COMPONENTS; j++){
123 sum1 = sum1 + x_i[j]*g_jk[i+N_COMPONENTS*j]*beta_jk[i+N_COMPONENTS*j];
124 sum2 = sum2 + x_i[j]*beta_jk[i+N_COMPONENTS*j];
125 }
126 G_E = G_E + x_i[i] * sum1/sum2;
127 }
128 return G_E;
129 };
130
131 template <typename X_I, typename B, typename V=typename X_I::value_t>
132 V b_mix(const X_I& x_i, const B& b) any_platform {
133 V b_m = 0;
134 for(unsigned i = 0; i < N_COMPONENTS; i++){
135 b_m += b[i]*x_i[i];
136 }
137 return b_m;
138 };
139
140 template <typename X_I, typename A_I, typename B, typename V=typename X_I::value_t>
141 V a_mix(const X_I& x_i, const A_I& a_i, const B& b, V b_m, V G_E) any_platform {
142 V a_sum = 0;
143 for(unsigned i = 0; i < N_COMPONENTS; i++){
144 a_sum = a_sum + x_i[i]*a_i[i]/b[i];
145 }
146 /*double a_m = 0;
147 for(unsigned i = 0; i < N_COMPONENTS; i++){
148 for(unsigned j = 0; j < N_COMPONENTS; j++){
149 a_m += x_i[i]*x_i[j]*util::pow(a_i[i]*a_i[j],0.5);
150 }
151 }
152 return a_m;*/
153 return b_m * (a_sum - G_E / 0.6232252401);
154 };
155
156 template <typename RHO_FIELD, typename __T, typename K, typename A_C, typename B, typename T_C, typename M, typename ALPHA, typename G_I, typename G_II, typename BIG_M, typename V=typename RHO_FIELD::value_t>
157 V compute(const RHO_FIELD& rhoField, const __T& _T, const K& k, const A_C& a_c, const B& b, const T_C& T_c, const M& m, const ALPHA& alpha, const G_I& g_I, const G_II& g_II, const BIG_M& _M) any_platform {
158 Vector<V,N_COMPONENTS> molarRhoField{}, a_i{}, x_i{};
159 Vector<V,N_COMPONENTS*N_COMPONENTS> g_jk{}, beta_jk{};
160 V R = 1;
161 V molarRho = 0;
162 V Rho = 0;
163 for(unsigned i = 0; i < N_COMPONENTS; i++){
164 a_i[i] = a_c[i] * util::pow((1 + m[i] * (1 - util::pow((_T/T_c[i]),0.5))),2);
165 for(unsigned j = 0; j < N_COMPONENTS; j++){
166 g_jk[N_COMPONENTS*i+j] = g_I[N_COMPONENTS*i+j] + g_II[N_COMPONENTS*i+j] * _T;
167 beta_jk[N_COMPONENTS*i+j] = b[j] * util::exp(-alpha[N_COMPONENTS*i+j]*g_jk[N_COMPONENTS*i+j]/(R*_T));
168 }
169 molarRhoField[i] = rhoField[i]/_M[i];
170 molarRho += molarRhoField[i];
171 Rho += rhoField[i];
172 }
173 for(unsigned i = 0; i < N_COMPONENTS; i++){
174 x_i[i] = molarRhoField[i]/molarRho;
175 }
176 V g_E = G_Excess(x_i, g_jk, beta_jk);
177 V b_m = b_mix(x_i, b);
178 V a_m = a_mix(x_i, a_i, b, b_m, g_E);
179 V p = molarRho*R*_T/(1-b_m*molarRho) - a_m*molarRho*molarRho/(1+2*b_m*molarRho-b_m*molarRho*b_m*molarRho);
180 V psi = util::sqrt(-6.*(k*p - Rho/3.));
181 return psi;
182 };
183
184 template <typename RHO_FIELD, typename A_C, typename B, typename T_C, typename M, typename ALPHA, typename G_I, typename G_II, typename BIG_M, typename V=typename RHO_FIELD::value_t>
185 V computeP(const RHO_FIELD& rhoField, V _T, const A_C& a_c, const B& b, const T_C& T_c, const M& m, const ALPHA& alpha, const G_I& g_I, const G_II& g_II, const BIG_M& _M) any_platform {
186 olb::Vector<V,N_COMPONENTS> molarRhoField{}, a_i{}, x_i{};
188 V R = 1;
189 V molarRho = 0;
190 for(unsigned i = 0; i < N_COMPONENTS; i++){
191 a_i[i] = a_c[i] * util::pow((1 + m[i] * (1 - util::pow((_T/T_c[i]),0.5))),2);
192 for(unsigned j = 0; j < N_COMPONENTS; j++){
193 g_jk[N_COMPONENTS*i+j] = g_I[N_COMPONENTS*i+j] + g_II[N_COMPONENTS*i+j] * _T;
194 beta_jk[N_COMPONENTS*i+j] = b[j] * util::exp(-alpha[N_COMPONENTS*i+j]*g_jk[N_COMPONENTS*i+j]/(R*_T));
195 }
196 molarRhoField[i] = rhoField[i]/_M[i];
197 molarRho += molarRhoField[i];
198 }
199 for(unsigned i = 0; i < N_COMPONENTS; i++){
200 x_i[i] = molarRhoField[i]/molarRho;
201 }
202 V g_E = G_Excess(x_i, g_jk, beta_jk);
203 V b_m = b_mix(x_i, b);
204 V a_m = a_mix(x_i, a_i, b, b_m, g_E);
205 V p = molarRho*R*_T/(1-b_m*molarRho) - a_m*molarRho*molarRho/(1+2*b_m*molarRho-b_m*molarRho*b_m*molarRho);
206 return p;
207 }
208
209};
210
211}
212
213
214// established -- original for both single- and multicomponent flow
215
216template <typename T, typename S>
217class ShanChen93 : public AnalyticalF<1,T,S> {
218private:
219 T _rhoZero;
220public:
221 ShanChen93(T rhoZero=1.);
222 bool operator() (T psi[], const S rho[]);
223};
224
225// established -- only multicomponent flow
226
227template <typename T, typename S>
228class PsiEqualsRho : public AnalyticalF<1,T,S> {
229private:
230public:
231 PsiEqualsRho();
232 bool operator() (T psi[], const S rho[]) override;
233};
234
235// established -- only singlecomponent flow
236
237template <typename T, typename S>
238class ShanChen94 : public AnalyticalF<1,T,S> {
239private:
240 T _rhoZero;
241 T _psiZero;
242public:
243 ShanChen94(T rhoZero=200., T psiZero=4.);
244 bool operator() (T psi[], const S rho[]) override;
245};
246
247template <typename T, typename S>
248class PengRobinson : public AnalyticalF<1,T,S> {
249private:
250 T _G;
251 T _acentricFactor;
252 T _a;
253 T _b;
254 T _R;
255 T _alpha;
256 T _t;
257 T _tc;
258public:
259 PengRobinson(T G, T acentricFactor=0.334, T a=2./49., T b=2./21., T tr=.8);
260 bool operator() (T psi[], const S rho[]);
261 // second operator allows to incorporate temperature changes
262 bool operator() (T psi[], const S rho[], const S t[]);
263};
264
265template <typename T, typename S>
266class CarnahanStarling : public AnalyticalF<1,T,S> {
267private:
268 T _G;
269 T _a;
270 T _b;
271 T _R;
272 T _t;
273public:
274 CarnahanStarling(T G, T a=1., T b=4., T tr=.7);
275 bool operator() (T psi[], const S rho[]) override;
276 // second operator allows to incorporate temperature changes
277 bool operator() (T psi[], const S rho[], const S t[]);
278};
279
280// under development -- for singlecomponent flow
281
282// 0.5 -> psiZero=0.65
283// 1 -> psiZero=1.9
284// 1.5 -> psiZero=3.5
285// 2. -> psiZero=5,45
286template <typename T, typename S>
287class Krause : public AnalyticalF<1,T,S> {
288private:
289 T _rhoZero;
290 T _psiZero;
291public:
292 Krause(T rhoZero=1., T psiZero=1.9);
293 bool operator() (T psi[], const S rho[]);
294};
295
296template <typename T, typename S> // density of liquid phase always equals rhoZero for G=-1
297class WeisbrodKrause : public AnalyticalF<1,T,S> {
298private:
299 T _rhoZero;
300 T _sigmu;
301public:
302 WeisbrodKrause(T rhoZero=1., T sigmu=1.);
303 bool operator() (T psi[], const S rho[]);
304};
305
306template <typename T, typename S> // not very good
307class Normal : public AnalyticalF<1,T,S> {
308private:
309 T _sigma;
310 T _mu;
311public:
312 Normal(T sigma=1., T mu=1.);
313 bool operator() (T psi[], const S rho[]);
314};
315
316} // end namespace olb
317
318#endif
319
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)
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
Plain old scalar vector.
Definition vector.h:47
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
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.
#define any_platform
Define preprocessor macros for device-side functions, constant storage.
Definition platform.h:78
Base of a field whose size is defined by [C,U_1,...,U_N]^T * [1,V_1,...V_N].
V compute(V rho, PARAMETERS &params) any_platform
V b_mix(const X_I &x_i, const B &b) any_platform
V compute(const RHO_FIELD &rhoField, const __T &_T, const K &k, const A_C &a_c, const B &b, const T_C &T_c, const M &m, const ALPHA &alpha, const G_I &g_I, const G_II &g_II, const BIG_M &_M) any_platform
V a_mix(const X_I &x_i, const A_I &a_i, const B &b, V b_m, V G_E) any_platform
V computeP(const RHO_FIELD &rhoField, V _T, const A_C &a_c, const B &b, const T_C &T_c, const M &m, const ALPHA &alpha, const G_I &g_I, const G_II &g_II, const BIG_M &_M) any_platform
V G_Excess(const X_I &x_i, const G_JK &g_jk, const BETA_JK &beta_jk) any_platform
V compute(V rho, PARAMETERS &params) any_platform
V compute(V rho, PARAMETERS &params) any_platform
Plain wrapper for list of types.
Definition meta.h:276