1 /**
2 * @file multiGasTransport.cpp
3 * test problem for multi transport
4 */
5
6 // This file is part of Cantera. See License.txt in the top-level directory or
7 // at https://cantera.org/license.txt for license and copyright information.
8
9 // Example
10 //
11 // Test case for mixture transport in a gas
12 // The basic idea is to set up a gradient of some kind.
13 // Then the resulting transport coefficients out.
14 // Essentially all of the interface routines should be
15 // exercised and the results dumped out.
16 //
17 // A blessed solution test will make sure that the actual
18 // solution doesn't change as a function of time or
19 // further development.
20
21 // perhaps, later, an analytical solution could be added
22
23 #include "cantera/thermo/IdealGasPhase.h"
24 #include "cantera/transport/MultiTransport.h"
25 #include "cantera/base/Solution.h"
26
27 #include <iostream>
28
29 using namespace std;
30 using namespace Cantera;
31
cutoff(double val,double atol=1.0E-15)32 static double cutoff(double val, double atol=1.0E-15)
33 {
34 if (fabs(val) < atol) {
35 return 0.0;
36 }
37 return val;
38 }
39
main(int argc,char ** argv)40 int main(int argc, char** argv)
41 {
42 #if defined(_MSC_VER) && _MSC_VER < 1900
43 _set_output_format(_TWO_DIGIT_EXPONENT);
44 #endif
45 try {
46 auto sol = newSolution("gri30.yaml", "gri30", "Multi");
47 auto gas = sol->thermo();
48 size_t nsp = gas->nSpecies();
49 double pres = 1.0E5;
50 vector_fp Xset(nsp, 0.0);
51 Xset[0] = 0.269205 ;
52 Xset[1] = 0.000107082;
53 Xset[2] = 1.36377e-09 ;
54 Xset[3] = 4.35475e-10;
55 Xset[4] = 4.34036e-06 ;
56 Xset[5] = 0.192249;
57 Xset[6] = 3.59356e-13;
58 Xset[7] = 2.78061e-12 ;
59 Xset[8] = 4.7406e-18 ;
60 Xset[9] = 4.12955e-17 ;
61 Xset[10] = 2.58549e-14 ;
62 Xset[11] = 8.96502e-16 ;
63 Xset[12] = 6.09056e-11 ;
64 Xset[13] = 7.56752e-09 ;
65 Xset[14] = 0.192253;
66 Xset[15] = 0.0385036;
67 Xset[16] = 1.49596e-08 ;
68 Xset[17] = 2.22378e-08 ;
69 Xset[18] = 1.43096e-13 ;
70 Xset[19] = 1.45312e-15 ;
71 Xset[20] = 1.96948e-12 ;
72 Xset[21] = 8.41937e-19;
73 Xset[22] = 3.18852e-13 ;
74 Xset[23] = 7.93625e-18 ;
75 Xset[24] = 3.20653e-15 ;
76 Xset[25] = 1.15149e-19 ;
77 Xset[26] = 1.61189e-18 ;
78 Xset[27] = 1.4719e-15 ;
79 Xset[28] = 5.24728e-13 ;
80 Xset[29] = 6.90582e-17 ;
81 Xset[30] = 6.37248e-12 ;
82 Xset[31] = 5.93728e-11 ;
83 Xset[32] = 2.71219e-09 ;
84 Xset[33] = 2.66645e-06 ;
85 Xset[34] = 6.57142e-11 ;
86 Xset[35] = 9.52453e-08 ;
87 Xset[36] = 1.26006e-14;
88 Xset[37] = 3.49802e-12;
89 Xset[38] = 1.19232e-11 ;
90 Xset[39] = 7.17782e-13;
91 Xset[40] = 1.85347e-07 ;
92 Xset[41] = 8.25325e-14 ;
93 Xset[42] = 5.00914e-20 ;
94 Xset[43] = 1.54407e-16 ;
95 Xset[44] =3.07176e-11 ;
96 Xset[45] =4.93198e-08 ;
97 Xset[46] =4.84792e-12 ;
98 Xset[47] = 0.307675 ;
99 Xset[48] =0;
100 Xset[49] =6.21649e-29;
101 Xset[50] = 8.42393e-28 ;
102 Xset[51] = 6.77865e-18;
103 Xset[52] = 2.19225e-16;
104 double T1 = 1500.;
105
106 double sum = 0.0;
107 for (size_t k = 0; k < nsp; k++) {
108 sum += Xset[k];
109 }
110 for (size_t k = 0; k < nsp; k++) {
111 Xset[k] /= sum;
112 }
113
114 vector_fp X2set(nsp, 0.0);
115 X2set[0] = 0.25 ;
116 X2set[5] = 0.17;
117 X2set[14] = 0.15;
118 X2set[15] = 0.05;
119 X2set[47] = 0.38 ;
120 double T2 = 1200.;
121
122 double dist = 0.1;
123
124 vector_fp X3set(nsp, 0.0);
125 X3set[0] = 0.27 ;
126 X3set[5] = 0.15;
127 X3set[14] = 0.18;
128 X3set[15] = 0.06;
129 X3set[47] = 0.36 ;
130 double T3 = 1400.;
131
132 vector_fp grad_T(3, 0.0);
133 Array2D grad_X(nsp, 2, 0.0);
134
135 for (size_t k = 0; k < nsp; k++) {
136 grad_X(k,0) = (X2set[k] - Xset[k])/dist;
137 grad_X(k,1) = (X3set[k] - Xset[k])/dist;
138 }
139
140 grad_T[0] = (T2 - T1) / dist;
141 grad_T[1] = (T3 - T1) / dist;
142
143 auto tran = sol->transport();
144 gas->setState_TPX(1500.0, pres, Xset.data());
145 vector_fp mixDiffs(nsp, 0.0);
146
147 tran->getMixDiffCoeffs(mixDiffs.data());
148 printf(" Dump of the mixture Diffusivities:\n");
149 for (size_t k = 0; k < nsp; k++) {
150 string sss = gas->speciesName(k);
151 printf(" %15s %13.2g\n", sss.c_str(), mixDiffs[k]);
152 }
153
154 vector_fp specVisc(nsp, 0.0);
155 tran->getSpeciesViscosities(specVisc.data());
156 printf(" Dump of the species viscosities:\n");
157 for (size_t k = 0; k < nsp; k++) {
158 string sss = gas->speciesName(k);
159 printf(" %15s %13.4g\n", sss.c_str(), specVisc[k]);
160 }
161
162 vector_fp thermDiff(nsp, 0.0);
163 tran->getThermalDiffCoeffs(thermDiff.data());
164 printf(" Dump of the Thermal Diffusivities :\n");
165 for (size_t k = 0; k < nsp; k++) {
166 string sss = gas->speciesName(k);
167 double ddd = cutoff(thermDiff[k]);
168 printf(" %15s %13.4g\n", sss.c_str(), ddd);
169 }
170
171 printf("Viscosity and thermal Cond vs. T\n");
172 for (size_t k = 0; k < 10; k++) {
173 T1 = 400. + 100. * k;
174 gas->setState_TPX(T1, pres, Xset.data());
175 double visc = tran->viscosity();
176 double cond = tran->thermalConductivity();
177 printf(" %13g %13.4g %13.4g\n", T1, visc, cond);
178 }
179
180 gas->setState_TPX(T1, pres, Xset.data());
181
182 Array2D Bdiff(nsp, nsp, 0.0);
183 printf("Binary Diffusion Coefficients H2 vs species\n");
184
185 tran->getBinaryDiffCoeffs(nsp, Bdiff.ptrColumn(0));
186 for (size_t k = 0; k < nsp; k++) {
187 string sss = gas->speciesName(k);
188 printf(" H2 - %15s %13.4g %13.4g\n", sss.c_str(), Bdiff(0,k), Bdiff(k,0));
189 }
190
191 vector_fp specMob(nsp, 0.0);
192
193 printf(" Dump of the species mobilities:\n");
194 for (size_t k = 0; k < nsp; k++) {
195 string sss = gas->speciesName(k);
196 printf(" %15s %13.4g\n", sss.c_str(), specMob[k]);
197 }
198
199 Array2D fluxes(nsp, 2, 0.0);
200 tran->getSpeciesFluxes(2, grad_T.data(), nsp,
201 grad_X.ptrColumn(0), nsp, fluxes.ptrColumn(0));
202 printf(" Dump of the species fluxes:\n");
203 double sum1 = 0.0;
204 double sum2 = 0.0;
205 double max1 = 0.0;
206 double max2 = 0.0;
207 for (size_t k = 0; k < nsp; k++) {
208 string sss = gas->speciesName(k);
209 double ddd = cutoff(fluxes(k,0));
210 double eee = cutoff(fluxes(k,1));
211 printf(" %15s %13.4g %13.4g\n", sss.c_str(), ddd, eee);
212 sum1 += fluxes(k,0);
213 if (fabs(fluxes(k,0)) > max1) {
214 max1 = fabs(fluxes(k,0));
215 }
216 sum2 += fluxes(k,1);
217 if (fabs(fluxes(k,1)) > max2) {
218 max2 = fabs(fluxes(k,0));
219 }
220 }
221
222 // Make sure roundoff error doesn't interfere with the printout.
223 // these should be zero.
224 if (fabs(sum1) * 1.0E14 > max1) {
225 printf("sum in x direction = %13.4g\n", sum1);
226 } else {
227 printf("sum in x direction = 0\n");
228 }
229 if (fabs(sum2) * 1.0E14 > max2) {
230 printf("sum in y direction = %13.4g\n", sum1);
231 } else {
232 printf("sum in y direction = 0\n");
233 }
234
235 Array2D MDdiff(nsp, nsp, 0.0);
236 printf("Multicomponent Diffusion Coefficients H2 vs species\n");
237
238 tran->getMultiDiffCoeffs(nsp, MDdiff.ptrColumn(0));
239 for (size_t k = 0; k < nsp; k++) {
240 string sss = gas->speciesName(k);
241 printf(" H2 - %15s %13.4g %13.4g\n", sss.c_str(), MDdiff(0,k), MDdiff(k,0));
242 }
243
244 } catch (CanteraError& err) {
245 std::cout << err.what() << std::endl;
246 return -1;
247 }
248
249 return 0;
250 }
251 /***********************************************************/
252