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