diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp index 4ffb881a3..5bfa039c3 100644 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -761,47 +761,45 @@ namespace Cantera { tr.xml->XML_comment(logfile, "*** polynomial coefficients not printed (log_level < 3) ***"); } #endif - for (i = 0; i < nsp; i++) - { - for (j = i; j < nsp; j++) - { - // Chemkin fits only delta* = 0 - if (mode != CK_Mode) - dstar = tr.delta(i,j); - else - dstar = 0.0; + for (i = 0; i < nsp; i++) { + for (j = i; j < nsp; j++) { + // Chemkin fits only delta* = 0 + if (mode != CK_Mode) { + dstar = tr.delta(i,j); + } else { + dstar = 0.0; + } - // if a fit has already been generated for - // delta* = tr.delta(i,j), then use it. Otherwise, - // make a new fit, and add tr.delta(i,j) to the list - // of delta* values for which fits have been done. + // if a fit has already been generated for + // delta* = tr.delta(i,j), then use it. Otherwise, + // make a new fit, and add tr.delta(i,j) to the list + // of delta* values for which fits have been done. - // 'find' returns a pointer to end() if not found - if (dptr = find(tr.fitlist.begin(), tr.fitlist.end(), - dstar), dptr == tr.fitlist.end()) - { - vector_fp ca(degree+1), cb(degree+1), cc(degree+1); - vector_fp co22(degree+1); - m_integrals->fit(logfile, degree, dstar, - DATA_PTR(ca), DATA_PTR(cb), DATA_PTR(cc)); - m_integrals->fit_omega22(logfile, degree, dstar, - DATA_PTR(co22)); - tr.omega22_poly.push_back(co22); - tr.astar_poly.push_back(ca); - tr.bstar_poly.push_back(cb); - tr.cstar_poly.push_back(cc); - tr.poly[i][j] = static_cast(tr.astar_poly.size()) - 1; - tr.fitlist.push_back(dstar); - } + // 'find' returns a pointer to end() if not found + dptr = find(tr.fitlist.begin(), tr.fitlist.end(), dstar); + if (dptr == tr.fitlist.end()) { + vector_fp ca(degree+1), cb(degree+1), cc(degree+1); + vector_fp co22(degree+1); + m_integrals->fit(logfile, degree, dstar, + DATA_PTR(ca), DATA_PTR(cb), DATA_PTR(cc)); + m_integrals->fit_omega22(logfile, degree, dstar, + DATA_PTR(co22)); + tr.omega22_poly.push_back(co22); + tr.astar_poly.push_back(ca); + tr.bstar_poly.push_back(cb); + tr.cstar_poly.push_back(cc); + tr.poly[i][j] = static_cast(tr.astar_poly.size()) - 1; + tr.fitlist.push_back(dstar); + } - // delta* found in fitlist, so just point to this - // polynomial - else { - tr.poly[i][j] = static_cast((dptr - tr.fitlist.begin())); - } - tr.poly[j][i] = tr.poly[i][j]; - } - } + // delta* found in fitlist, so just point to this + // polynomial + else { + tr.poly[i][j] = static_cast((dptr - tr.fitlist.begin())); + } + tr.poly[j][i] = tr.poly[i][j]; + } + } #ifdef DEBUG_MODE if (m_verbose) { tr.xml->XML_close(logfile, "tstar_fits"); @@ -1462,71 +1460,69 @@ namespace Cantera { mxerr = 0.0, mxrelerr = 0.0; vector_fp diff(np + 1); doublereal eps, sigma; - for (k = 0; k < tr.nsp_; k++) - { - for (j = k; j < tr.nsp_; j++) { + for (k = 0; k < tr.nsp_; k++) { + for (j = k; j < tr.nsp_; j++) { - ipoly = tr.poly[k][j]; - for (n = 0; n < np; n++) { + ipoly = tr.poly[k][j]; + for (n = 0; n < np; n++) { - t = tr.tmin + dt*n; + t = tr.tmin + dt*n; - eps = tr.epsilon(j,k); - tstar = Boltzmann * t/eps; - sigma = tr.diam(j,k); - om11 = m_integrals->omega11(tstar, tr.delta(j,k)); + eps = tr.epsilon(j,k); + tstar = Boltzmann * t/eps; + sigma = tr.diam(j,k); + om11 = m_integrals->omega11(tstar, tr.delta(j,k)); - diffcoeff = ThreeSixteenths * - sqrt( 2.0 * Pi/tr.reducedMass(k,j) ) * - pow((Boltzmann * t), 1.5)/ - (Pi * sigma * sigma * om11); + diffcoeff = ThreeSixteenths * + sqrt( 2.0 * Pi/tr.reducedMass(k,j) ) * + pow((Boltzmann * t), 1.5)/ + (Pi * sigma * sigma * om11); - // 2nd order correction - // NOTE: THIS CORRECTION IS NOT APPLIED - doublereal fkj, fjk; - getBinDiffCorrection(t, tr, k, j, 1.0, 1.0, fkj, fjk); - //diffcoeff *= fkj; + // 2nd order correction + // NOTE: THIS CORRECTION IS NOT APPLIED + doublereal fkj, fjk; + getBinDiffCorrection(t, tr, k, j, 1.0, 1.0, fkj, fjk); + //diffcoeff *= fkj; - if (mode == CK_Mode) { - diff[n] = log(diffcoeff); - w[n] = -1.0; - } - else { - diff[n] = diffcoeff/pow(t, 1.5); - w[n] = 1.0/(diff[n]*diff[n]); - } + if (mode == CK_Mode) { + diff[n] = log(diffcoeff); + w[n] = -1.0; } - polyfit(np, DATA_PTR(tlog), DATA_PTR(diff), - DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c)); - - doublereal pre; - for (n = 0; n < np; n++) { - if (mode == CK_Mode) { - val = exp(diff[n]); - fit = exp(poly3(tlog[n], DATA_PTR(c))); - } - else { - t = exp(tlog[n]); - pre = pow(t, 1.5); - val = pre * diff[n]; - fit = pre * poly4(tlog[n], DATA_PTR(c)); - } - err = fit - val; - relerr = err/val; - if (fabs(err) > mxerr) mxerr = fabs(err); - if (fabs(relerr) > mxrelerr) mxrelerr = fabs(relerr); + else { + diff[n] = diffcoeff/pow(t, 1.5); + w[n] = 1.0/(diff[n]*diff[n]); } - tr.diffcoeffs.push_back(c); -#ifdef DEBUG_MODE - if (tr.log_level >= 2 && m_verbose) { - tr.xml->XML_writeVector(logfile, " ", tr.thermo->speciesName(k) - + "__"+tr.thermo->speciesName(j), c.size(), DATA_PTR(c)); - } -#endif } + polyfit(np, DATA_PTR(tlog), DATA_PTR(diff), + DATA_PTR(w), degree, ndeg, 0.0, DATA_PTR(c)); + + doublereal pre; + for (n = 0; n < np; n++) { + if (mode == CK_Mode) { + val = exp(diff[n]); + fit = exp(poly3(tlog[n], DATA_PTR(c))); + } else { + t = exp(tlog[n]); + pre = pow(t, 1.5); + val = pre * diff[n]; + fit = pre * poly4(tlog[n], DATA_PTR(c)); + } + err = fit - val; + relerr = err/val; + if (fabs(err) > mxerr) mxerr = fabs(err); + if (fabs(relerr) > mxrelerr) mxrelerr = fabs(relerr); + } + tr.diffcoeffs.push_back(c); +#ifdef DEBUG_MODE + if (tr.log_level >= 2 && m_verbose) { + tr.xml->XML_writeVector(logfile, " ", tr.thermo->speciesName(k) + + "__"+tr.thermo->speciesName(j), c.size(), DATA_PTR(c)); + } +#endif } + } #ifdef DEBUG_MODE if (m_verbose) { sprintf(s,"Maximum binary diffusion coefficient absolute error:" diff --git a/Cantera/src/transport/TransportParams.h b/Cantera/src/transport/TransportParams.h index dd42f2c86..de03f32bd 100644 --- a/Cantera/src/transport/TransportParams.h +++ b/Cantera/src/transport/TransportParams.h @@ -86,12 +86,12 @@ namespace Cantera { public: + //! Constructor GasTransportParams() : TransportParams(), visccoeffs(0), condcoeffs(0), diffcoeffs(0), - polytempvec(0), poly(0), omega22_poly(0), astar_poly(0), @@ -112,32 +112,89 @@ namespace Cantera { { } + //! Destructor virtual ~GasTransportParams() {} // polynomial fits //! temperature-fit of the viscosity + /*! + * The outer loop the number of species, nsp + * The inner loop is over degree + 1, which is the polynomial order of the collision integral fit. + */ std::vector visccoeffs; - //temperature-fit heat conduction - std::vector condcoeffs; + //! temperature-fits of the heat conduction + /*! + * The outer loop the number of species, nsp + * The inner loop is over degree + 1, which is the polynomial order of the collision integral fit. + */ + std::vector condcoeffs; - //temperature-fit diffusivity - std::vector diffcoeffs; - vector_fp polytempvec; + //! temperature-fits of the diffusivity + /*! + * The outer loop the number of species, nsp + * The inner loop is over degree + 1, which is the polynomial order of the collision integral fit. + */ + std::vector diffcoeffs; + + //! This is vector of vectors containing the integer ookup value for the (i,j) interaction + /*! + * The outer loop is over a flat (i,j) index that is parameterized on the tr.delta(i,j) value. + * Unique values of delta get their own spot in the array. The values of delta are storred in + * the fitlist vector. + * + * The inner loop is over degree + 1, which is the polynomial order of the collision integral fit. + */ std::vector > poly; - std::vector omega22_poly; - std::vector astar_poly; - std::vector bstar_poly; - std::vector cstar_poly; + + //! This is vector of vectors containing the astar fit. + /*! + * The outer loop is over a flat (i,j) index that is parameterized on the tr.delta(i,j) value. + * Unique values of delta get their own spot in the array. The values of delta are storred in + * the fitlist vector. + * + * The inner loop is over degree + 1, which is the polynomial order of the collision integral fit. + */ + std::vector omega22_poly; + + //! This is vector of vectors containing the astar fit. + /*! + * The outer loop is over a flat (i,j) index that is parameterized on the tr.delta(i,j) value. + * Unique values of delta get their own spot in the array. The values of delta are storred in + * the fitlist vector. + * + * The inner loop is over degree + 1, which is the polynomial order of the collision integral fit. + */ + std::vector astar_poly; + + //! This is vector of vectors containing the astar fit. + /*! + * The outer loop is over a flat (i,j) index that is parameterized on the tr.delta(i,j) value. + * Unique values of delta get their own spot in the array. The values of delta are storred in + * the fitlist vector. + * + * The inner loop is over degree + 1, which is the polynomial order of the collision integral fit. + */ + std::vector bstar_poly; + + //! This is vector of vectors containing the astar fit. + /*! + * The outer loop is over a flat (i,j) index that is parameterized on the tr.delta(i,j) value. + * Unique values of delta get their own spot in the array. The values of delta are storred in + * the fitlist vector. + * + * The inner loop is over degree + 1, which is the polynomial order of the collision integral fit. + */ + std::vector cstar_poly; //! Rotational relaxation number for the species in the current phase /*! * length is the number of species in the phase * units are dimensionless */ - vector_fp zrot; + vector_fp zrot; //! Dimensionless rotational heat capacity of the species in the current phase /*! @@ -145,31 +202,93 @@ namespace Cantera { * length is the number of species in the pahse * units are dimensionless (Cr / R) */ - vector_fp crot; + vector_fp crot; + //! Vector of booleans indicating whether a species is a polar molecule + /*! + * Length is nsp + */ std::vector polar; + + //! Polarizability of each species in the phase + /*! + * Length = nsp + * Units = m^3 + */ vector_fp alpha; - vector_fp fitlist; + + //! This is vector containing the values of delta(i,j) that are used in the collision integral fits. + /*! + * This is used in astar_poly, bstar_poly, cstar_poly, and omega22_poly. + * The outer loop is over a flat (i,j) index that is parameterized on the tr.delta(i,j) value. + * Unique values of delta get their own spot in the array. The values of delta are storred in + * the fitlist vector. + * + */ + vector_fp fitlist; //! Lennard-Jones well-depth of the species in the current phase /*! * length is the number of species in the phase - * Units are Joules (Note this is not Joules/kmol) + * Units are Joules (Note this is not Joules/kmol) (note, no kmol -> this is a per molecule amount) */ - vector_fp eps; + vector_fp eps; //! Lennard-Jones diameter of the species in the current phase /*! * length is the number of species in the phase * units are in meters. */ - vector_fp sigma; + vector_fp sigma; + //! This is the reduced mass of the interaction between species i and j + /*! + * tr.reducedMass(i,j) = tr.mw[i] * tr.mw[j] / (Avogadro * (tr.mw[i] + tr.mw[j])); + * + * Units are kg (note, no kmol -> this is a per molecule amount) + * + * Length nsp * nsp. This is a symmetric matrix + */ DenseMatrix reducedMass; - DenseMatrix diam; - DenseMatrix epsilon; - DenseMatrix dipole; - DenseMatrix delta; + + //! hard-sphere diameter for (i,j) collision + /*! + * diam(i,j) = 0.5*(tr.sigma[i] + tr.sigma[j]); + * Units are m (note, no kmol -> this is a per molecule amount) + * + * Length nsp * nsp. This is a symmetric matrix. + */ + DenseMatrix diam; + + //! The effective well depth for (i,j) collisions + /*! + * epsilon(i,j) = sqrt(tr.eps[i]*tr.eps[j]); + * Units are Joules (note, no kmol -> this is a per molecule amount) + * + * Length nsp * nsp. This is a symmetric matrix. + */ + DenseMatrix epsilon; + + //! The effective dipole moment for (i,j) collisions + /*! + * tr.dipoleMoment has units of Debye's. A Debye is 10-18 cm3/2 erg1/2 + * + * tr.dipole(i,i) = 1.e-25 * SqrtTen * trdat.dipoleMoment; + * tr.dipole(i,j) = sqrt(tr.dipole(i,i)*tr.dipole(j,j)); + * Units are in Debye (note, no kmol -> this is a per molecule amount) + * + * Length nsp * nsp. This is a symmetric matrix. + */ + DenseMatrix dipole; + + //! Matrix containing the reduced dipole moment of the interaction between two species + /*! + * This is the reduced dipole moment of the interaction between two species + * 0.5 * tr.dipole(i,j)*tr.dipole(i,j) (epsilon(i,j) * d * d * d); + * + * Length nsp * nsp .This is a symmetric matrix + */ + DenseMatrix delta; };