From a20321e5c9211cf9e3e7e38f81cbe6ad05b7c918 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Fri, 19 Mar 2010 16:44:17 +0000 Subject: [PATCH] Worked on documenting using doxygen. --- Cantera/src/transport/L_matrix.h | 372 ++++++++++++----------- Cantera/src/transport/MultiTransport.cpp | 14 +- 2 files changed, 199 insertions(+), 187 deletions(-) diff --git a/Cantera/src/transport/L_matrix.h b/Cantera/src/transport/L_matrix.h index 238d68370..ecf5cfe59 100644 --- a/Cantera/src/transport/L_matrix.h +++ b/Cantera/src/transport/L_matrix.h @@ -25,229 +25,231 @@ namespace Cantera { - // #define CHEMKIN_COMPATIBILITY_MODE + // #define CHEMKIN_COMPATIBILITY_MODE - const doublereal Min_C_Internal = 0.001; + //! Constant to compare dimensionless heat capacities against zero + const doublereal Min_C_Internal = 0.001; - bool MultiTransport::hasInternalModes(int j) { + bool MultiTransport::hasInternalModes(int j) { #ifdef CHEMKIN_COMPATIBILITY_MODE - return (m_crot[j] > Min_C_Internal); + return (m_crot[j] > Min_C_Internal); #else - return (m_cinternal[j] > Min_C_Internal); + return (m_cinternal[j] > Min_C_Internal); #endif + } + + + /** + * Evaluate the upper-left block of the L matrix. + */ + void MultiTransport::eval_L0000(const doublereal* x) { + + doublereal prefactor = 16.0*m_temp/25.0; + doublereal sum; + int i, j, k; + for (i = 0; i < m_nsp; i++) { + // subtract-off the k=i term to account for the first delta + // function in Eq. (12.121) + + sum = -x[i]/m_bdiff(i,i); + for (k = 0; k < m_nsp; k++) sum += x[k]/m_bdiff(i,k); + + sum /= m_mw[i]; + for (j = 0; j != m_nsp; ++j) { + m_Lmatrix(i,j) = prefactor * x[j] + * ( m_mw[j] * sum + x[i]/m_bdiff(i,j) ); + } + // diagonal term is zero + m_Lmatrix(i,i) = 0.0; } + } - /** - * Evaluate the upper-left block of the L matrix. - */ - void MultiTransport::eval_L0000(const doublereal* x) { + //////////////////////////////////////////////////////////////////////////// - doublereal prefactor = 16.0*m_temp/25.0; - doublereal sum; - int i, j, k; - for (i = 0; i < m_nsp; i++) - { - // subtract-off the k=i term to account for the first delta - // function in Eq. (12.121) - sum = -x[i]/m_bdiff(i,i); - for (k = 0; k < m_nsp; k++) sum += x[k]/m_bdiff(i,k); + void MultiTransport::eval_L0010(const doublereal* x) { - sum /= m_mw[i]; - for (j = 0; j != m_nsp; ++j) { - m_Lmatrix(i,j) = prefactor * x[j] - * ( m_mw[j] * sum + x[i]/m_bdiff(i,j) ); - } - // diagonal term is zero - m_Lmatrix(i,i) = 0.0; - } + doublereal prefactor = 1.6*m_temp; + + doublereal sum, wj, xj; + int i, j; + for (j = 0; j < m_nsp; j++) { + //constant = prefactor * x[j]; + xj = x[j]; + wj = m_mw[j]; + sum = 0.0; + for (i = 0; i < m_nsp; i++) { + m_Lmatrix(i,j + m_nsp) = - prefactor * x[i] * xj * m_mw[i] * + (1.2 * m_cstar(j,i) - 1.0) / + ( (wj + m_mw[i]) * m_bdiff(j,i) ); + + // the next term is independent of "j"; + // need to do it for the "j,j" term + sum -= m_Lmatrix(i,j+m_nsp); + } + m_Lmatrix(j,j+m_nsp) += sum; } + } - //////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////// - void MultiTransport::eval_L0010(const doublereal* x) { - - doublereal prefactor = 1.6*m_temp; - - doublereal sum, wj, xj; - int i, j; - for (j = 0; j < m_nsp; j++) { - //constant = prefactor * x[j]; - xj = x[j]; - wj = m_mw[j]; - sum = 0.0; - for (i = 0; i < m_nsp; i++) { - m_Lmatrix(i,j + m_nsp) = - prefactor * x[i] * xj * m_mw[i] * - (1.2 * m_cstar(j,i) - 1.0) / - ( (wj + m_mw[i]) * m_bdiff(j,i) ); - - // the next term is independent of "j"; - // need to do it for the "j,j" term - sum -= m_Lmatrix(i,j+m_nsp); - } - m_Lmatrix(j,j+m_nsp) += sum; - } + void MultiTransport::eval_L1000() { + int i, j; + for (j = 0; j < m_nsp; j++) { + for (i = 0; i < m_nsp; i++) { + m_Lmatrix(i+m_nsp,j) = m_Lmatrix(j,i+m_nsp); + } } + } + ////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////// + void MultiTransport::eval_L1010(const doublereal* x) { - void MultiTransport::eval_L1000() { - int i, j; - for (j = 0; j < m_nsp; j++) - for (i = 0; i < m_nsp; i++) - m_Lmatrix(i+m_nsp,j) = m_Lmatrix(j,i+m_nsp); - } + const doublereal fiveover3pi = 5.0/(3.0*Pi); + doublereal prefactor = (16.0*m_temp)/25.0; + int i, j; + doublereal constant1, wjsq, constant2, constant3, constant4, + fourmj, threemjsq, sum, sumwij;; + doublereal term1, term2; - ////////////////////////////////////////////////////////////////////// + for (j = 0; j < m_nsp; j++) { - void MultiTransport::eval_L1010(const doublereal* x) { + // get constant terms that depend on just species "j" - const doublereal fiveover3pi = 5.0/(3.0*Pi); - doublereal prefactor = (16.0*m_temp)/25.0; + constant1 = prefactor*x[j]; + wjsq = m_mw[j]*m_mw[j]; + constant2 = 13.75*wjsq; + constant3 = m_crot[j]/m_rotrelax[j]; + constant4 = 7.5*wjsq; + fourmj = 4.0*m_mw[j]; + threemjsq = 3.0*m_mw[j]*m_mw[j]; + sum = 0.0; + for (i = 0; i < m_nsp; i++) { - int i, j; - doublereal constant1, wjsq, constant2, constant3, constant4, - fourmj, threemjsq, sum, sumwij;; - doublereal term1, term2; + sumwij = m_mw[i] + m_mw[j]; + term1 = m_bdiff(i,j) * sumwij*sumwij; + term2 = fourmj*m_astar(i,j)*(1.0 + fiveover3pi* + (constant3 + + (m_crot[i]/m_rotrelax[i]))); // see Eq. (12.125) - for (j = 0; j < m_nsp; j++) { + m_Lmatrix(i+m_nsp,j+m_nsp) = constant1*x[i]*m_mw[i] /(m_mw[j]*term1) * + (constant2 - threemjsq*m_bstar(i,j) + - term2*m_mw[j]); - // get constant terms that depend on just species "j" - - constant1 = prefactor*x[j]; - wjsq = m_mw[j]*m_mw[j]; - constant2 = 13.75*wjsq; - constant3 = m_crot[j]/m_rotrelax[j]; - constant4 = 7.5*wjsq; - fourmj = 4.0*m_mw[j]; - threemjsq = 3.0*m_mw[j]*m_mw[j]; - sum = 0.0; - for (i = 0; i < m_nsp; i++) { - - sumwij = m_mw[i] + m_mw[j]; - term1 = m_bdiff(i,j) * sumwij*sumwij; - term2 = fourmj*m_astar(i,j)*(1.0 + fiveover3pi* - (constant3 + - (m_crot[i]/m_rotrelax[i]))); // see Eq. (12.125) - - m_Lmatrix(i+m_nsp,j+m_nsp) = constant1*x[i]*m_mw[i] /(m_mw[j]*term1) * - (constant2 - threemjsq*m_bstar(i,j) - - term2*m_mw[j]); - - sum += x[i] /(term1) * - (constant4 + m_mw[i]*m_mw[i]* - (6.25 - 3.0*m_bstar(i,j)) + term2*m_mw[i]); - } + sum += x[i] /(term1) * + (constant4 + m_mw[i]*m_mw[i]* + (6.25 - 3.0*m_bstar(i,j)) + term2*m_mw[i]); + } - m_Lmatrix(j+m_nsp,j+m_nsp) -= sum*constant1; - } + m_Lmatrix(j+m_nsp,j+m_nsp) -= sum*constant1; } + } - ////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////////////////////// - void MultiTransport::eval_L1001(const doublereal* x) { + void MultiTransport::eval_L1001(const doublereal* x) { - doublereal prefactor = 32.00*m_temp/(5.00*Pi); - int i,j; - doublereal constant, sum; - int n2 = 2*m_nsp; - int npoly = 0; - for (j = 0; j < m_nsp; j++) { - // collect terms that depend only on "j" - if (hasInternalModes(j)) { - constant = prefactor*m_mw[j]*x[j]*m_crot[j]/(m_cinternal[j]*m_rotrelax[j]); - sum = 0.0; - for (i = 0; i < m_nsp; i++) { - // see Eq. (12.127) - m_Lmatrix(i+m_nsp,j+n2) = constant * m_astar(j,i) * x[i] / - ( (m_mw[j] + m_mw[i] ) * m_bdiff(j,i)); - sum += m_Lmatrix(i+m_nsp,j+n2); - } - npoly++; - m_Lmatrix(j+m_nsp,j+n2) += sum; - } - else { - for (i = 0; i < m_nsp; i++) m_Lmatrix(i+m_nsp,j+n2) = 0.0; - } - } + doublereal prefactor = 32.00*m_temp/(5.00*Pi); + int i,j; + doublereal constant, sum; + int n2 = 2*m_nsp; + int npoly = 0; + for (j = 0; j < m_nsp; j++) { + // collect terms that depend only on "j" + if (hasInternalModes(j)) { + constant = prefactor*m_mw[j]*x[j]*m_crot[j]/(m_cinternal[j]*m_rotrelax[j]); + sum = 0.0; + for (i = 0; i < m_nsp; i++) { + // see Eq. (12.127) + m_Lmatrix(i+m_nsp,j+n2) = constant * m_astar(j,i) * x[i] / + ( (m_mw[j] + m_mw[i] ) * m_bdiff(j,i)); + sum += m_Lmatrix(i+m_nsp,j+n2); + } + npoly++; + m_Lmatrix(j+m_nsp,j+n2) += sum; + } + else { + for (i = 0; i < m_nsp; i++) m_Lmatrix(i+m_nsp,j+n2) = 0.0; + } } + } - //////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////// - void MultiTransport::eval_L0001() { - int i, j; - int n2 = 2*m_nsp; - for (j = 0; j < m_nsp; j++) - for (i = 0; i < m_nsp; i++) - m_Lmatrix(i,j+n2) = 0.0; - } - - //////////////////////////////////////////////////////////////////////// - - void MultiTransport::eval_L0100() { - int i, j; - int n2 = 2*m_nsp; - for (j = 0; j < m_nsp; j++) - for (i = 0; i < m_nsp; i++) - m_Lmatrix(i+n2,j) = 0.0; // see Eq. (12.123) - } - - //////////////////////////////////////////////////////////////////////// - - void MultiTransport::eval_L0110() { - int i, j; - int n2 = 2*m_nsp; - for (j = 0; j < m_nsp; j++) - for (i = 0; i < m_nsp; i++) - m_Lmatrix(i+n2,j+m_nsp) = m_Lmatrix(j+m_nsp,i+n2); // see Eq. (12.123) - } - - //////////////////////////////////////////////////////////////////////// - - - void MultiTransport::eval_L0101(const doublereal* x) { - - const doublereal fivepi = 5.00*Pi; - const doublereal eightoverpi = 8.0 / Pi; - - doublereal prefactor = 4.00*m_temp; - int n2 = 2*m_nsp; - int i,k; - doublereal constant1, constant2, diff_int, sum; - for (i = 0; i < m_nsp; i++) { - if (hasInternalModes(i)) { - // collect terms that depend only on "i" - constant1 = prefactor*x[i]/m_cinternal[i]; - constant2 = 12.00*m_mw[i]*m_crot[i] / - (fivepi*m_cinternal[i]*m_rotrelax[i]); - sum = 0.0; - for (k = 0; k < m_nsp; k++) { - // see Eq. (12.131) - diff_int = m_bdiff(i,k); - m_Lmatrix(k+n2,i+n2) = 0.0; - sum += x[k]/diff_int; - if (k != i) sum += x[k]*m_astar(i,k)*constant2 / - (m_mw[k]*diff_int); - } - // see Eq. (12.130) - m_Lmatrix(i+n2,i+n2) = - - eightoverpi*m_mw[i]*x[i]*x[i]*m_crot[i] / - (m_cinternal[i]*m_cinternal[i]*GasConstant*m_visc[i]*m_rotrelax[i]) - - constant1*sum; - } - else { - for (k = 0; k < m_nsp; k++) - m_Lmatrix(i+n2,i+n2) = 1.0; - } - } + void MultiTransport::eval_L0001() { + int i, j; + int n2 = 2*m_nsp; + for (j = 0; j < m_nsp; j++) + for (i = 0; i < m_nsp; i++) + m_Lmatrix(i,j+n2) = 0.0; + } + + //////////////////////////////////////////////////////////////////////// + + void MultiTransport::eval_L0100() { + int i, j; + int n2 = 2*m_nsp; + for (j = 0; j < m_nsp; j++) + for (i = 0; i < m_nsp; i++) + m_Lmatrix(i+n2,j) = 0.0; // see Eq. (12.123) + } + + //////////////////////////////////////////////////////////////////////// + + void MultiTransport::eval_L0110() { + int i, j; + int n2 = 2*m_nsp; + for (j = 0; j < m_nsp; j++) + for (i = 0; i < m_nsp; i++) + m_Lmatrix(i+n2,j+m_nsp) = m_Lmatrix(j+m_nsp,i+n2); // see Eq. (12.123) + } + + //////////////////////////////////////////////////////////////////////// + + + void MultiTransport::eval_L0101(const doublereal* x) { + + const doublereal fivepi = 5.00*Pi; + const doublereal eightoverpi = 8.0 / Pi; + + doublereal prefactor = 4.00*m_temp; + int n2 = 2*m_nsp; + int i,k; + doublereal constant1, constant2, diff_int, sum; + for (i = 0; i < m_nsp; i++) { + if (hasInternalModes(i)) { + // collect terms that depend only on "i" + constant1 = prefactor*x[i]/m_cinternal[i]; + constant2 = 12.00*m_mw[i]*m_crot[i] / + (fivepi*m_cinternal[i]*m_rotrelax[i]); + sum = 0.0; + for (k = 0; k < m_nsp; k++) { + // see Eq. (12.131) + diff_int = m_bdiff(i,k); + m_Lmatrix(k+n2,i+n2) = 0.0; + sum += x[k]/diff_int; + if (k != i) sum += x[k]*m_astar(i,k)*constant2 / + (m_mw[k]*diff_int); + } + // see Eq. (12.130) + m_Lmatrix(i+n2,i+n2) = + - eightoverpi*m_mw[i]*x[i]*x[i]*m_crot[i] / + (m_cinternal[i]*m_cinternal[i]*GasConstant*m_visc[i]*m_rotrelax[i]) + - constant1*sum; + } + else { + for (k = 0; k < m_nsp; k++) + m_Lmatrix(i+n2,i+n2) = 1.0; + } } + } } #endif diff --git a/Cantera/src/transport/MultiTransport.cpp b/Cantera/src/transport/MultiTransport.cpp index 4bc72954c..ebf8b047e 100644 --- a/Cantera/src/transport/MultiTransport.cpp +++ b/Cantera/src/transport/MultiTransport.cpp @@ -1010,9 +1010,19 @@ namespace Cantera { m_bdiff(k,k) = d; } - // internal heat capacities + // Calculate the internal heat capacities by subtracting off the translational contributions + /* + * HKM Exploratory comment: + * The translational component is 1.5 + * The rotational component is 1.0 for a linear molecule and 1.5 for a nonlinear molecule + * and zero for a monotomic. + * Chemkin has traditionally subtracted 1.5 here (SAND86-8246). + * The original Dixon-Lewis paper subtracted 1.5 here. + */ const array_fp& cp = ((IdealGasPhase*)m_thermo)->cp_R_ref(); - for (k = 0; k < m_nsp; k++) m_cinternal[k] = cp[k] - 2.5; + for (k = 0; k < m_nsp; k++) { + m_cinternal[k] = cp[k] - 2.5; + } } /**