Cantera  2.0
L_matrix.h
Go to the documentation of this file.
1 /**
2  * @file L_matrix.h
3  * Functions to evaluate portions of the L matrix needed for
4  * multicomponent transport properties.
5  */
6 
7 #ifndef CT_LMATRIX_H
8 #define CT_LMATRIX_H
9 
10 #include "cantera/numerics/DenseMatrix.h"
11 #include "cantera/base/ct_defs.h"
12 
13 #include <vector>
14 
15 /////////////////////////////////////////////////////////////////////
16 
17 namespace Cantera
18 {
19 //====================================================================================================================
20 // #define CHEMKIN_COMPATIBILITY_MODE
21 
22 //! Constant to compare dimensionless heat capacities against zero
23 const doublereal Min_C_Internal = 0.001;
24 //====================================================================================================================
25 bool MultiTransport::hasInternalModes(size_t j)
26 {
27 #ifdef CHEMKIN_COMPATIBILITY_MODE
28  return (m_crot[j] > Min_C_Internal);
29 #else
30  return (m_cinternal[j] > Min_C_Internal);
31 #endif
32 }
33 
34 //====================================================================================================================
35 /*
36  * Evaluate the upper-left block of the L matrix.
37  */
38 void MultiTransport::eval_L0000(const doublereal* const x)
39 {
40 
41  doublereal prefactor = 16.0*m_temp/25.0;
42  doublereal sum;
43  for (size_t i = 0; i < m_nsp; i++) {
44  // subtract-off the k=i term to account for the first delta
45  // function in Eq. (12.121)
46 
47  sum = -x[i]/m_bdiff(i,i);
48  for (size_t k = 0; k < m_nsp; k++) {
49  sum += x[k]/m_bdiff(i,k);
50  }
51 
52  sum /= m_mw[i];
53  for (size_t j = 0; j != m_nsp; ++j) {
54  m_Lmatrix(i,j) = prefactor * x[j]
55  * (m_mw[j] * sum + x[i]/m_bdiff(i,j));
56  }
57  // diagonal term is zero
58  m_Lmatrix(i,i) = 0.0;
59  }
60 }
61 //====================================================================================================================
62 void MultiTransport::eval_L0010(const doublereal* const x)
63 {
64 
65  doublereal prefactor = 1.6*m_temp;
66 
67  doublereal sum, wj, xj;
68  for (size_t j = 0; j < m_nsp; j++) {
69  //constant = prefactor * x[j];
70  xj = x[j];
71  wj = m_mw[j];
72  sum = 0.0;
73  for (size_t i = 0; i < m_nsp; i++) {
74  m_Lmatrix(i,j + m_nsp) = - prefactor * x[i] * xj * m_mw[i] *
75  (1.2 * m_cstar(j,i) - 1.0) /
76  ((wj + m_mw[i]) * m_bdiff(j,i));
77 
78  // the next term is independent of "j";
79  // need to do it for the "j,j" term
80  sum -= m_Lmatrix(i,j+m_nsp);
81  }
82  m_Lmatrix(j,j+m_nsp) += sum;
83  }
84 }
85 //====================================================================================================================
86 void MultiTransport::eval_L1000()
87 {
88  for (size_t j = 0; j < m_nsp; j++) {
89  for (size_t i = 0; i < m_nsp; i++) {
90  m_Lmatrix(i+m_nsp,j) = m_Lmatrix(j,i+m_nsp);
91  }
92  }
93 }
94 //====================================================================================================================
95 void MultiTransport::eval_L1010(const doublereal* x)
96 {
97 
98  const doublereal fiveover3pi = 5.0/(3.0*Pi);
99  doublereal prefactor = (16.0*m_temp)/25.0;
100 
101  doublereal constant1, wjsq, constant2, constant3, constant4,
102  fourmj, threemjsq, sum, sumwij;;
103  doublereal term1, term2;
104 
105  for (size_t j = 0; j < m_nsp; j++) {
106 
107  // get constant terms that depend on just species "j"
108 
109  constant1 = prefactor*x[j];
110  wjsq = m_mw[j]*m_mw[j];
111  constant2 = 13.75*wjsq;
112  constant3 = m_crot[j]/m_rotrelax[j];
113  constant4 = 7.5*wjsq;
114  fourmj = 4.0*m_mw[j];
115  threemjsq = 3.0*m_mw[j]*m_mw[j];
116  sum = 0.0;
117  for (size_t i = 0; i < m_nsp; i++) {
118 
119  sumwij = m_mw[i] + m_mw[j];
120  term1 = m_bdiff(i,j) * sumwij*sumwij;
121  term2 = fourmj*m_astar(i,j)*(1.0 + fiveover3pi*
122  (constant3 +
123  (m_crot[i]/m_rotrelax[i]))); // see Eq. (12.125)
124 
125  m_Lmatrix(i+m_nsp,j+m_nsp) = constant1*x[i]*m_mw[i] /(m_mw[j]*term1) *
126  (constant2 - threemjsq*m_bstar(i,j)
127  - term2*m_mw[j]);
128 
129  sum += x[i] /(term1) *
130  (constant4 + m_mw[i]*m_mw[i]*
131  (6.25 - 3.0*m_bstar(i,j)) + term2*m_mw[i]);
132  }
133 
134  m_Lmatrix(j+m_nsp,j+m_nsp) -= sum*constant1;
135  }
136 }
137 //====================================================================================================================
138 void MultiTransport::eval_L1001(const doublereal* x)
139 {
140 
141  doublereal prefactor = 32.00*m_temp/(5.00*Pi);
142  doublereal constant, sum;
143  size_t n2 = 2*m_nsp;
144  int npoly = 0;
145  for (size_t j = 0; j < m_nsp; j++) {
146  // collect terms that depend only on "j"
147  if (hasInternalModes(j)) {
148  constant = prefactor*m_mw[j]*x[j]*m_crot[j]/(m_cinternal[j]*m_rotrelax[j]);
149  sum = 0.0;
150  for (size_t i = 0; i < m_nsp; i++) {
151  // see Eq. (12.127)
152  m_Lmatrix(i+m_nsp,j+n2) = constant * m_astar(j,i) * x[i] /
153  ((m_mw[j] + m_mw[i]) * m_bdiff(j,i));
154  sum += m_Lmatrix(i+m_nsp,j+n2);
155  }
156  npoly++;
157  m_Lmatrix(j+m_nsp,j+n2) += sum;
158  } else {
159  for (size_t i = 0; i < m_nsp; i++) {
160  m_Lmatrix(i+m_nsp,j+n2) = 0.0;
161  }
162  }
163  }
164 }
165 //====================================================================================================================
166 
167 void MultiTransport::eval_L0001()
168 {
169  size_t n2 = 2*m_nsp;
170  for (size_t j = 0; j < m_nsp; j++) {
171  for (size_t i = 0; i < m_nsp; i++) {
172  m_Lmatrix(i,j+n2) = 0.0;
173  }
174  }
175 }
176 //====================================================================================================================
177 
178 void MultiTransport::eval_L0100()
179 {
180  size_t n2 = 2*m_nsp;
181  for (size_t j = 0; j < m_nsp; j++)
182  for (size_t i = 0; i < m_nsp; i++) {
183  m_Lmatrix(i+n2,j) = 0.0; // see Eq. (12.123)
184  }
185 }
186 //====================================================================================================================
187 
188 void MultiTransport::eval_L0110()
189 {
190  size_t n2 = 2*m_nsp;
191  for (size_t j = 0; j < m_nsp; j++)
192  for (size_t i = 0; i < m_nsp; i++) {
193  m_Lmatrix(i+n2,j+m_nsp) = m_Lmatrix(j+m_nsp,i+n2); // see Eq. (12.123)
194  }
195 }
196 //====================================================================================================================
197 void MultiTransport::eval_L0101(const doublereal* x)
198 {
199 
200  const doublereal fivepi = 5.00*Pi;
201  const doublereal eightoverpi = 8.0 / Pi;
202 
203  doublereal prefactor = 4.00*m_temp;
204  size_t n2 = 2*m_nsp;
205  doublereal constant1, constant2, diff_int, sum;
206  for (size_t i = 0; i < m_nsp; i++) {
207  if (hasInternalModes(i)) {
208  // collect terms that depend only on "i"
209  constant1 = prefactor*x[i]/m_cinternal[i];
210  constant2 = 12.00*m_mw[i]*m_crot[i] /
211  (fivepi*m_cinternal[i]*m_rotrelax[i]);
212  sum = 0.0;
213  for (size_t k = 0; k < m_nsp; k++) {
214  // see Eq. (12.131)
215  diff_int = m_bdiff(i,k);
216  m_Lmatrix(k+n2,i+n2) = 0.0;
217  sum += x[k]/diff_int;
218  if (k != i) sum += x[k]*m_astar(i,k)*constant2 /
219  (m_mw[k]*diff_int);
220  }
221  // see Eq. (12.130)
222  m_Lmatrix(i+n2,i+n2) =
223  - eightoverpi*m_mw[i]*x[i]*x[i]*m_crot[i] /
224  (m_cinternal[i]*m_cinternal[i]*GasConstant*m_visc[i]*m_rotrelax[i])
225  - constant1*sum;
226  } else {
227  for (size_t k = 0; k < m_nsp; k++) {
228  m_Lmatrix(i+n2,i+n2) = 1.0;
229  }
230  }
231  }
232 }
233 }
234 //======================================================================================================================
235 #endif
Generated by   doxygen 1.8.2