Line data Source code
1 : #ifndef MODEL_GAS_TRANSPORT_MIXTURE_AVERAGED_H_
2 : #define MODEL_GAS_TRANSPORT_MIXTURE_AVERAGED_H_
3 :
4 : #include <vector>
5 : #include <memory>
6 : #include "IO/ParmParse.H"
7 : #include "Model/Gas/Thermo/CpConstant.H"
8 : #include "Model/Gas/Transport/Transport.H"
9 :
10 : // Type: 0,1 determines if per species viscosity and thermal conductivity are
11 : // constant (type=0) or computed using kinetic theory (type=1)
12 : // Either way, mixture averaging is used to compute the mixture values
13 : //
14 : // For kinetic theory: Chapman-Enskog is used to compute viscosity
15 : // and Eucken relation used to compute conductivity
16 :
17 :
18 : // dynamic viscosity, mu, computed using Chapman-Enskog theory
19 : //
20 : // 5E20 sqrt(pi * m * kb * T)
21 : // mu = ----------------------------
22 : // 16 pi σ^2 Ω
23 : //
24 : // MW
25 : // where m = -- (mass per particle), MW: molecular weight, NA: Avogadro's number, kb: Boltzman constant
26 : // NA
27 : //
28 : // σ: Lennard-Jones collision diameter in Angstroms, Ω: collision integral
29 :
30 :
31 : // thermal conductivity coefficient, k, computed using Eucken relation
32 : //
33 : // k = mu * (cp + 1.25R), where cp and R are both per kg
34 : //
35 :
36 : namespace Model {
37 : namespace Gas {
38 : namespace Transport {
39 :
40 : class Mixture_Averaged // public Transport
41 : {
42 : public:
43 : static constexpr const char* name = "mixture_averaged";
44 :
45 : private:
46 : int nspecies;
47 : double Rg;
48 : std::vector<double> MW;
49 : Thermo::Thermo<Thermo::CpConstant> thermo;
50 : std::string type;
51 : std::vector<double> val1; // constant viscosity mu or Lennard-Jones collision diameter depending on type
52 : std::vector<double> val2; // constant thermal conductivity or Lennard-Jones potential well depth depending on type
53 :
54 : public:
55 8 : Mixture_Averaged()
56 8 : {
57 8 : Rg = Set::Constant::Rg;
58 8 : }
59 :
60 8 : static void Parse( Mixture_Averaged & value, IO::ParmParse & pp,
61 : int a_nspecies, std::vector<double> &a_MW, Thermo::Thermo<Thermo::CpConstant> a_thermo)
62 : {
63 8 : value.nspecies = a_nspecies;
64 8 : value.MW = a_MW;
65 8 : value.thermo = a_thermo;
66 :
67 : // Compute individual species transport properties as "constant" or from "LJ" (Lennard-Jones) theory
68 32 : pp.query_validate("type", value.type, {"constant", "LJ"});
69 8 : if (value.type == "constant")
70 : {
71 : // Dynamic viscosity if type=constant
72 16 : pp.queryarr_required("mu", value.val1, Unit::Pressure() * Unit::Time());
73 : // Thermal conductivity if type=constant
74 24 : pp.queryarr_required("k", value.val2, Unit::Power() / Unit::Length() / Unit::Temperature());
75 : }
76 0 : else if (value.type == "LJ")
77 : {
78 : // Lennard-Jones diameter if type=LJ
79 0 : pp.queryarr_required("LJdiameter", value.val1, Unit::Length());
80 : // Lennard-Jones well depth if type=LJ
81 0 : pp.queryarr_required("LJwelldepth", value.val2, Unit::Temperature());
82 : }
83 : else
84 : {
85 0 : Util::Abort(INFO, "Unrecognized type for transport model ", value.type);
86 : }
87 : //pp.query_default("P_1atm", value.P_1atm, "1_atm", Unit::Pressure());
88 :
89 56 : Util::Assert(INFO, TEST(value.val1.size() == (size_t)value.nspecies));
90 56 : Util::Assert(INFO, TEST(value.val2.size() == (size_t)value.nspecies));
91 8 : }
92 :
93 : private:
94 0 : double collision_integral(double Tstar) const {
95 0 : double S = 1.16145 * pow(Tstar, -0.14874)
96 0 : + 0.52487 / exp(0.77320 * Tstar)
97 0 : + 2.16178 / exp(2.43787 * Tstar);
98 0 : return S;
99 : }
100 0 : double mu_enskog(double T, int i) const {
101 0 : return 2.669570e-6 * sqrt(MW[i] * T) / (val1[i] * val1[i] * collision_integral(T / val2[i]));
102 : }
103 0 : double k_eucken(double T, Set::Patch<const Set::Scalar>& X, int i, int j, int k, int n) const
104 : {
105 0 : double cp = thermo.cp_mol(T, X, i, j, k);
106 0 : return mu_enskog(T, n) * (X(i, j, k, n) * cp + 1.25 * Rg) / MW[n];
107 : }
108 :
109 : public:
110 13235200 : double dynamic_viscosity(double T, Set::Patch<const Set::Scalar>& X, int i, int j, int k) const
111 : {
112 : // Wilke's Average: Dynamic viscosity, Pa-s
113 13235200 : double mu_mix = 0.0;
114 13235200 : if ( type == "constant" )
115 : {
116 : // dynamic viscosity, mu, is constant per species
117 26470400 : for (int a = 0; a < nspecies; ++a) {
118 13235200 : double mua = val1[a];
119 13235200 : double MWa = MW[a];
120 13235200 : double phi = 0.0;
121 26470400 : for (int b = 0; b < nspecies; ++b) {
122 13235200 : double mub = val1[b];
123 13235200 : double MWb = MW[b];
124 13235200 : double ratio = pow(MWb / MWa, 0.25);
125 13235200 : double mu_ratio = sqrt(mua / mub);
126 13235200 : double term = sqrt(1.0 + MWa / MWb);
127 13235200 : phi += X(i,j,k,b) * pow(1.0 + mu_ratio * ratio, 2.0) / (sqrt(8.0) * term);
128 : }
129 13235200 : mu_mix += X(i, j, k, a) * mua / phi;
130 : }
131 : }
132 0 : else if ( type == "LJ" )
133 : {
134 : // dynamic viscosity, mu, computed by Chapman-Enskog per species
135 0 : for (int a = 0; a < nspecies; ++a) {
136 0 : double mua = mu_enskog(T, a);
137 0 : double MWa = MW[a];
138 0 : double phi = 0.0;
139 0 : for (int b = 0; b < nspecies; ++b) {
140 0 : double mub = mu_enskog(T, b);
141 0 : double MWb = MW[b];
142 0 : double ratio = pow(MWb / MWa, 0.25);
143 0 : double mu_ratio = sqrt(mua / mub);
144 0 : double term = sqrt(1.0 + MWa / MWb);
145 0 : phi += X(i, j, k, b) * pow(1.0 + mu_ratio * ratio, 2.0) / (sqrt(8.0) * term);
146 : }
147 0 : mu_mix += X(i, j, k, a) * mua / phi;
148 : }
149 : }
150 : else
151 : {
152 0 : Util::Abort(INFO, "[Transport::Mixture_Averaged] Unknown type.");
153 : }
154 :
155 13235200 : if (!(mu_mix == mu_mix) ) mu_mix = 0.0;
156 13235200 : return mu_mix;
157 : }
158 :
159 0 : double thermal_conductivity(double T, Set::Patch<const Set::Scalar>& X, int i, int j, int k) const
160 : {
161 : // Wilke's Average: Thermal conductivity coefficient, W/(m-K)
162 0 : double k_mix = 0.0;
163 0 : if ( type == "constant" )
164 : {
165 : // thermal conductivity coefficient, k, is constant per species
166 0 : for (int a = 0; a < nspecies; ++a) {
167 0 : double mua = val1[a];
168 0 : double MWa = MW[a];
169 0 : double phi = 0.0;
170 0 : for (int b = 0; b < nspecies; ++b) {
171 0 : double mub = val1[b];
172 0 : double MWb = MW[b];
173 0 : double ratio = pow(MWb / MWa, 0.25);
174 0 : double mu_ratio = sqrt(mua / mub);
175 0 : double term = sqrt(1.0 + MWa / MWb);
176 0 : phi += X(i, j, k, b) * pow(1.0 + mu_ratio * ratio, 2.0) / (sqrt(8.0) * term);
177 : }
178 0 : k_mix += X(i, j, k, a) * val2[a] / phi;
179 : }
180 : }
181 0 : else if ( type == "LJ" )
182 : {
183 : // thermal conductivity coefficient, k, computed by Eucken relation per species
184 0 : for (int a = 0; a < nspecies; ++a) {
185 0 : double mua = mu_enskog(T, a);
186 0 : double MWa = MW[a];
187 0 : double phi = 0.0;
188 0 : for (int b = 0; b < nspecies; ++b) {
189 0 : double mub = mu_enskog(T, b);
190 0 : double MWb = MW[b];
191 0 : double ratio = pow(MWb / MWa, 0.25);
192 0 : double mu_ratio = sqrt(mua / mub);
193 0 : double term = sqrt(1.0 + MWa / MWb);
194 0 : phi += X(i, j, k, b) * pow(1.0 + mu_ratio * ratio, 2.0) / (sqrt(8.0) * term);
195 : }
196 0 : k_mix += X(i, j, k, a) * k_eucken(T, X, i, j, k, a) / phi;
197 : }
198 : }
199 : else
200 : {
201 0 : Util::Abort(INFO, "[Transport::Mixture_Averaged] Unknown type.");
202 : }
203 :
204 0 : return k_mix;
205 : }
206 :
207 0 : void diffusion_coeffs(
208 : Set::Patch<Set::Scalar>& DKM, double T, double P,
209 : Set::Patch<const Set::Scalar>& X, int i, int j, int k) const
210 : {
211 : // Mixture diffusion coefficients (Chapman–Enskog), m^2/s
212 0 : for (int a = 0; a < nspecies; ++a) {
213 0 : double sumD = 0.0;
214 0 : for (int b = 0; b < nspecies; ++b) {
215 0 : if (a == b) continue;
216 :
217 0 : double sigmaAB = 0.5 * (val1[a] + val1[b]);
218 0 : double epsAB = sqrt(val2[a] * val2[b]);
219 0 : double Tstar = T / epsAB;
220 :
221 : double omegaAB =
222 0 : 1.06036 / pow(Tstar, 0.15610) +
223 0 : 0.19300 / exp(0.47635 * Tstar) +
224 0 : 1.03587 / exp(1.52996 * Tstar) +
225 0 : 1.76474 / exp(3.89411 * Tstar);
226 :
227 : double DAB =
228 0 : 0.0018583 * sqrt(T*T*T * (1.0 / MW[a] + 1.0 / MW[b])) /
229 0 : ( (P / 101325.0) * sigmaAB*sigmaAB * omegaAB );
230 0 : DAB /= 1.0e4;
231 :
232 0 : sumD += X(i, j, k, b) / DAB;
233 : }
234 :
235 0 : DKM(i, j, k, a) = (1.0 - X(i, j, k ,a)) / sumD;
236 0 : if (!(DKM(i, j, k, a) == DKM(i, j, k, a))) DKM(i, j, k, a) = 0.0; // NaN guard
237 : }
238 0 : }
239 :
240 : }; // class Mixture_Averaged
241 :
242 : } // namespace Transport
243 : } // namespace Gas
244 : } // namespace Model
245 :
246 : #endif
|
