25#ifndef INTERACTION_POTENTIAL_H
26#define INTERACTION_POTENTIAL_H
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);
53 std::vector<double> z;
56 const double lambda{0.6232252401};
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};
73namespace interaction {
79 template <
typename V,
typename PARAMETERS>
88 template <
typename V,
typename PARAMETERS>
102 template <
typename V,
typename PARAMETERS>
104 V c = params.template get<B>() * rho / V{4};
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>());
114template <
unsigned N_COMPONENTS>
116 template <
typename X_I,
typename G_JK,
typename BETA_JK,
typename V=
typename X_I::value_t>
119 for(
unsigned i = 0; i < N_COMPONENTS; i++){
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];
126 G_E = G_E + x_i[i] * sum1/sum2;
131 template <
typename X_I,
typename B,
typename V=
typename X_I::value_t>
134 for(
unsigned i = 0; i < N_COMPONENTS; i++){
140 template <
typename X_I,
typename A_I,
typename B,
typename V=
typename X_I::value_t>
143 for(
unsigned i = 0; i < N_COMPONENTS; i++){
144 a_sum = a_sum + x_i[i]*a_i[i]/b[i];
153 return b_m * (a_sum - G_E / 0.6232252401);
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 {
163 for(
unsigned i = 0; i < N_COMPONENTS; i++){
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));
169 molarRhoField[i] = rhoField[i]/_M[i];
170 molarRho += molarRhoField[i];
173 for(
unsigned i = 0; i < N_COMPONENTS; i++){
174 x_i[i] = molarRhoField[i]/molarRho;
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);
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 {
190 for(
unsigned i = 0; i < N_COMPONENTS; i++){
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));
196 molarRhoField[i] = rhoField[i]/_M[i];
197 molarRho += molarRhoField[i];
199 for(
unsigned i = 0; i < N_COMPONENTS; i++){
200 x_i[i] = molarRhoField[i]/molarRho;
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);
216template <
typename T,
typename S>
227template <
typename T,
typename S>
232 bool operator() (T psi[],
const S rho[])
override;
237template <
typename T,
typename S>
244 bool operator() (T psi[],
const S rho[])
override;
247template <
typename T,
typename S>
259 PengRobinson(T G, T acentricFactor=0.334, T a=2./49., T b=2./21., T tr=.8);
262 bool operator() (T psi[],
const S rho[],
const S t[]);
265template <
typename T,
typename S>
275 bool operator() (T psi[],
const S rho[])
override;
277 bool operator() (T psi[],
const S rho[],
const S t[]);
286template <
typename T,
typename S>
292 Krause(T rhoZero=1., T psiZero=1.9);
296template <
typename T,
typename S>
306template <
typename T,
typename S>
312 Normal(T sigma=1., T mu=1.);
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
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)
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
ADf< T, DIM > exp(const ADf< T, DIM > &a)
Top level namespace for all of OpenLB.
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 ¶ms) 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 ¶ms) any_platform
V compute(V rho, PARAMETERS ¶ms) any_platform