diff --git a/src/transport/DustyGasTransport.cpp b/src/transport/DustyGasTransport.cpp index 8d29c9497..2ed6cd013 100644 --- a/src/transport/DustyGasTransport.cpp +++ b/src/transport/DustyGasTransport.cpp @@ -136,9 +136,8 @@ void DustyGasTransport::updateKnudsenDiffCoeffs() return; } doublereal K_g = m_pore_radius * m_porosity / m_tortuosity; - const doublereal TwoThirds = 2.0/3.0; for (size_t k = 0; k < m_nsp; k++) { - m_dk[k] = TwoThirds * K_g * sqrt((8.0 * GasConstant * m_temp)/ + m_dk[k] = 2.0/3.0 * K_g * sqrt((8.0 * GasConstant * m_temp)/ (Pi * m_mw[k])); } m_knudsen_ok = true; @@ -148,7 +147,6 @@ void DustyGasTransport::eval_H_matrix() { updateBinaryDiffCoeffs(); updateKnudsenDiffCoeffs(); - doublereal sum; for (size_t k = 0; k < m_nsp; k++) { // evaluate off-diagonal terms for (size_t j = 0; j < m_nsp; j++) { @@ -156,7 +154,7 @@ void DustyGasTransport::eval_H_matrix() } // evaluate diagonal term - sum = 0.0; + double sum = 0.0; for (size_t j = 0; j < m_nsp; j++) { if (j != k) { sum += m_x[j]/m_d(k,j); @@ -171,7 +169,6 @@ void DustyGasTransport::getMolarFluxes(const doublereal* const state1, const doublereal delta, doublereal* const fluxes) { - doublereal conc1, conc2; // cbar will be the average concentration between the two points doublereal* const cbar = m_spwork.data(); doublereal* const gradc = m_spwork2.data(); @@ -184,8 +181,8 @@ void DustyGasTransport::getMolarFluxes(const doublereal* const state1, doublereal c1sum = 0.0, c2sum = 0.0; for (size_t k = 0; k < m_nsp; k++) { - conc1 = rho1 * y1[k] / m_mw[k]; - conc2 = rho2 * y2[k] / m_mw[k]; + double conc1 = rho1 * y1[k] / m_mw[k]; + double conc2 = rho2 * y2[k] / m_mw[k]; cbar[k] = 0.5*(conc1 + conc2); gradc[k] = (conc2 - conc1) / delta; c1sum += conc1; diff --git a/src/transport/GasTransport.cpp b/src/transport/GasTransport.cpp index c91027488..5078dd308 100644 --- a/src/transport/GasTransport.cpp +++ b/src/transport/GasTransport.cpp @@ -154,7 +154,6 @@ doublereal GasTransport::viscosity() void GasTransport::updateViscosity_T() { - doublereal vratiokj, wratiojk, factor1; if (!m_spvisc_ok) { updateSpeciesViscosities(); } @@ -162,11 +161,11 @@ void GasTransport::updateViscosity_T() // see Eq. (9-5.15) of Reid, Prausnitz, and Poling for (size_t j = 0; j < m_nsp; j++) { for (size_t k = j; k < m_nsp; k++) { - vratiokj = m_visc[k]/m_visc[j]; - wratiojk = m_mw[j]/m_mw[k]; + double vratiokj = m_visc[k]/m_visc[j]; + double wratiojk = m_mw[j]/m_mw[k]; // Note that m_wratjk(k,j) holds the square root of m_wratjk(j,k)! - factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j); + double factor1 = 1.0 + (m_sqvisc[k]/m_sqvisc[j]) * m_wratjk(k,j); m_phi(k,j) = factor1*factor1 / (sqrt(8.0) * m_wratkj1(j,k)); m_phi(j,k) = m_phi(k,j)/(vratiokj * wratiojk); } @@ -507,8 +506,6 @@ void GasTransport::makePolarCorrections(size_t i, size_t j, void GasTransport::fitCollisionIntegrals(MMCollisionInt& integrals) { - double dstar; - // Chemkin fits to sixth order polynomials int degree = (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE); if (m_log_level) { @@ -523,11 +520,7 @@ void GasTransport::fitCollisionIntegrals(MMCollisionInt& integrals) for (size_t i = 0; i < m_nsp; i++) { for (size_t j = i; j < m_nsp; j++) { // Chemkin fits only delta* = 0 - if (m_mode != CK_Mode) { - dstar = m_delta(i,j); - } else { - dstar = 0.0; - } + double dstar = (m_mode != CK_Mode) ? m_delta(i,j) : 0.0; // if a fit has already been generated for delta* = m_delta(i,j), // then use it. Otherwise, make a new fit, and add m_delta(i,j) to @@ -577,7 +570,7 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) if (m_log_level && m_log_level < 2) { writelog("*** polynomial coefficients not printed (log_level < 2) ***\n"); } - double sqrt_T, visc, err, relerr, + double visc, err, relerr, mxerr = 0.0, mxrelerr = 0.0, mxerr_cond = 0.0, mxrelerr_cond = 0.0; if (m_log_level) { @@ -590,9 +583,6 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) } } - double cp_R, cond, w_RT, f_int, A_factor, B_factor, c1, cv_rot, cv_int, - f_rot, f_trans, om11, diffcoeff; - const vector_fp& mw = m_thermo->molecularWeights(); for (size_t k = 0; k < m_nsp; k++) { for (size_t n = 0; n < np; n++) { @@ -600,32 +590,31 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) m_thermo->setTemperature(t); vector_fp cp_R_all(m_thermo->nSpecies()); m_thermo->getCp_R_ref(&cp_R_all[0]); - cp_R = cp_R_all[k]; + double cp_R = cp_R_all[k]; double tstar = Boltzmann * t/ m_eps[k]; - sqrt_T = sqrt(t); + double sqrt_T = sqrt(t); double om22 = integrals.omega22(tstar, m_delta(k,k)); - om11 = integrals.omega11(tstar, m_delta(k,k)); + double om11 = integrals.omega11(tstar, m_delta(k,k)); // self-diffusion coefficient, without polar corrections - diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,k)) * - pow((Boltzmann * t), 1.5)/ - (Pi * m_sigma[k] * m_sigma[k] * om11); + double diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,k)) * + pow((Boltzmann * t), 1.5)/ + (Pi * m_sigma[k] * m_sigma[k] * om11); // viscosity visc = 5.0/16.0 * sqrt(Pi * mw[k] * Boltzmann * t / Avogadro) / (om22 * Pi * m_sigma[k]*m_sigma[k]); // thermal conductivity - w_RT = mw[k]/(GasConstant * t); - f_int = w_RT * diffcoeff/visc; - cv_rot = m_crot[k]; - A_factor = 2.5 - f_int; - B_factor = m_zrot[k] + 2.0/Pi * (5.0/3.0 * cv_rot + f_int); - c1 = 2.0/Pi * A_factor/B_factor; - cv_int = cp_R - 2.5 - cv_rot; - f_rot = f_int * (1.0 + c1); - f_trans = 2.5 * (1.0 - c1 * cv_rot/1.5); - cond = (visc/mw[k])*GasConstant*(f_trans * 1.5 + double f_int = mw[k]/(GasConstant * t) * diffcoeff/visc; + double cv_rot = m_crot[k]; + double A_factor = 2.5 - f_int; + double B_factor = m_zrot[k] + 2.0/Pi * (5.0/3.0 * cv_rot + f_int); + double c1 = 2.0/Pi * A_factor/B_factor; + double cv_int = cp_R - 2.5 - cv_rot; + double f_rot = f_int * (1.0 + c1); + double f_trans = 2.5 * (1.0 - c1 * cv_rot/1.5); + double cond = (visc/mw[k])*GasConstant*(f_trans * 1.5 + f_rot * cv_rot + f_int * cv_int); if (m_mode == CK_Mode) { @@ -660,7 +649,7 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) val = exp(spvisc[n]); fit = exp(poly3(tlog[n], c.data())); } else { - sqrt_T = exp(0.5*tlog[n]); + double sqrt_T = exp(0.5*tlog[n]); val = sqrt_T * pow(spvisc[n],2); fit = sqrt_T * pow(poly4(tlog[n], c.data()),2); } @@ -677,7 +666,7 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) val = exp(spcond[n]); fit = exp(poly3(tlog[n], c2.data())); } else { - sqrt_T = exp(0.5*tlog[n]); + double sqrt_T = exp(0.5*tlog[n]); val = sqrt_T * spcond[n]; fit = sqrt_T * poly4(tlog[n], c2.data()); } @@ -723,18 +712,16 @@ void GasTransport::fitProperties(MMCollisionInt& integrals) mxerr = 0.0, mxrelerr = 0.0; vector_fp diff(np + 1); - double eps, sigma; for (size_t k = 0; k < m_nsp; k++) { for (size_t j = k; j < m_nsp; j++) { for (size_t n = 0; n < np; n++) { double t = m_thermo->minTemp() + dt*n; - eps = m_epsilon(j,k); + double eps = m_epsilon(j,k); double tstar = Boltzmann * t/eps; - sigma = m_diam(j,k); - om11 = integrals.omega11(tstar, m_delta(j,k)); - diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,j)) * - pow(Boltzmann * t, 1.5) / - (Pi * sigma * sigma * om11); + double sigma = m_diam(j,k); + double om11 = integrals.omega11(tstar, m_delta(j,k)); + double diffcoeff = 3.0/16.0 * sqrt(2.0 * Pi/m_reducedMass(k,j)) + * pow(Boltzmann * t, 1.5) / (Pi * sigma * sigma * om11); // 2nd order correction // NOTE: THIS CORRECTION IS NOT APPLIED diff --git a/src/transport/HighPressureGasTransport.cpp b/src/transport/HighPressureGasTransport.cpp index 5e78e6589..776b561c8 100644 --- a/src/transport/HighPressureGasTransport.cpp +++ b/src/transport/HighPressureGasTransport.cpp @@ -132,7 +132,6 @@ void HighPressureGasTransport::getThermalDiffCoeffs(doublereal* const dt) void HighPressureGasTransport::getBinaryDiffCoeffs(const size_t ld, doublereal* const d) { - doublereal P_corr_ij, Tr_ij, Pr_ij; vector_fp PcP(5); size_t nsp = m_thermo->nSpecies(); vector_fp molefracs(nsp); @@ -161,9 +160,10 @@ void HighPressureGasTransport::getBinaryDiffCoeffs(const size_t ld, doublereal* x_j = x_j/(x_i + x_j); //Calculate Tr and Pr based on mole-fraction-weighted crit constants: - Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j)); - Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j)); + double Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j)); + double Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j)); + double P_corr_ij; if (Pr_ij < 0.1) { // If pressure is low enough, no correction is needed: P_corr_ij = 1; @@ -204,7 +204,6 @@ void HighPressureGasTransport::getMultiDiffCoeffs(const size_t ld, doublereal* c // Correct the binary diffusion coefficients for high-pressure effects; this // is basically the same routine used in 'getBinaryDiffCoeffs,' above: - doublereal P_corr_ij, Tr_ij, Pr_ij; size_t nsp = m_thermo->nSpecies(); vector_fp molefracs(nsp); m_thermo->getMoleFractions(&molefracs[0]); @@ -225,9 +224,10 @@ void HighPressureGasTransport::getMultiDiffCoeffs(const size_t ld, doublereal* c doublereal x_j = std::max(Tiny, molefracs[j]); x_i = x_i/(x_i+x_j); x_j = x_j/(x_i+x_j); - Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j)); - Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j)); + double Tr_ij = m_temp/(x_i*Tcrit_i(i) + x_j*Tcrit_i(j)); + double Pr_ij = m_thermo->pressure()/(x_i*Pcrit_i(i) + x_j*Pcrit_i(j)); + double P_corr_ij; if (Pr_ij < 0.1) { P_corr_ij = 1; }else { @@ -260,14 +260,12 @@ void HighPressureGasTransport::getMultiDiffCoeffs(const size_t ld, doublereal* c m_l0000_ok = false; // matrix is overwritten by inverse m_lmatrix_soln_ok = false; - doublereal pres = m_thermo->pressure(); doublereal prefactor = 16.0 * m_temp - *m_thermo->meanMolecularWeight()/(25.0*pres); - doublereal c; + *m_thermo->meanMolecularWeight()/(25.0*m_thermo->pressure()); for (size_t i = 0; i < m_nsp; i++) { for (size_t j = 0; j < m_nsp; j++) { - c = prefactor/m_mw[j]; + double c = prefactor/m_mw[j]; d[ld*j + i] = c*molefracs[i]*(m_Lmatrix(i,j) - m_Lmatrix(i,i)); } } @@ -280,7 +278,6 @@ doublereal HighPressureGasTransport::viscosity() double Pc_mix_n = 0.; double Pc_mix_d = 0.; double MW_mix = m_thermo->meanMolecularWeight(); - doublereal x_H, Tc, Zc, Tr, Afac, Z1m, Z2m; double MW_H = m_mw[0]; double MW_L = m_mw[0]; doublereal FP_mix_o = 0; @@ -291,13 +288,13 @@ doublereal HighPressureGasTransport::viscosity() vector_fp molefracs(nsp); m_thermo->getMoleFractions(&molefracs[0]); - x_H = molefracs[0]; + double x_H = molefracs[0]; for (size_t i = 0; i < m_nsp; i++) { // Calculate pure-species critical constants and add their contribution // to the mole-fraction-weighted mixture averages: - Tc = Tcrit_i(i); - Tr = tKelvin/Tc; - Zc = Zcrit_i(i); + double Tc = Tcrit_i(i); + double Tr = tKelvin/Tc; + double Zc = Zcrit_i(i); Tc_mix += Tc*molefracs[i]; Pc_mix_n += molefracs[i]*Zc; //numerator Pc_mix_d += molefracs[i]*Vcrit_i(i); //denominator @@ -345,17 +342,15 @@ doublereal HighPressureGasTransport::viscosity() *pow(Pc_mix,4)),1.0/6.0); if (ratio > 9 && x_H > 0.05 && x_H < 0.7) { - Afac = 1 - 0.01*pow(ratio,0.87); - } else { - Afac = 1; + FQ_mix_o *= 1 - 0.01*pow(ratio,0.87); } - FQ_mix_o *= Afac; // Calculate Z1m - Z1m = (0.807*pow(Tr_mix,0.618) - 0.357*exp(-0.449*Tr_mix) - + 0.340*exp(-4.058*Tr_mix)+0.018)*FP_mix_o*FQ_mix_o; + double Z1m = (0.807*pow(Tr_mix,0.618) - 0.357*exp(-0.449*Tr_mix) + + 0.340*exp(-4.058*Tr_mix)+0.018)*FP_mix_o*FQ_mix_o; // Calculate Z2m: + double Z2m; if (Tr_mix <= 1.0) { if (Pr_mix < Pvp_mix/Pc_mix) { doublereal alpha = 3.262 + 14.98*pow(Pr_mix,5.508); @@ -396,9 +391,8 @@ doublereal HighPressureGasTransport::viscosity() // Pure species critical properties - Tc, Pc, Vc, Zc: doublereal HighPressureGasTransport::Tcrit_i(size_t i) { - size_t nsp = m_thermo->nSpecies(); // Store current molefracs and set temp molefrac of species i to 1.0: - vector_fp molefracs = store(i,nsp); + vector_fp molefracs = store(i, m_thermo->nSpecies()); double tc = m_thermo->critTemperature(); // Restore actual molefracs: @@ -408,9 +402,8 @@ doublereal HighPressureGasTransport::Tcrit_i(size_t i) doublereal HighPressureGasTransport::Pcrit_i(size_t i) { - size_t nsp = m_thermo->nSpecies(); // Store current molefracs and set temp molefrac of species i to 1.0: - vector_fp molefracs = store(i,nsp); + vector_fp molefracs = store(i, m_thermo->nSpecies()); double pc = m_thermo->critPressure(); // Restore actual molefracs: @@ -420,9 +413,8 @@ doublereal HighPressureGasTransport::Pcrit_i(size_t i) doublereal HighPressureGasTransport::Vcrit_i(size_t i) { - size_t nsp = m_thermo->nSpecies(); // Store current molefracs and set temp molefrac of species i to 1.0: - vector_fp molefracs = store(i,nsp); + vector_fp molefracs = store(i, m_thermo->nSpecies()); double vc = m_thermo->critVolume(); // Restore actual molefracs: @@ -432,9 +424,8 @@ doublereal HighPressureGasTransport::Vcrit_i(size_t i) doublereal HighPressureGasTransport::Zcrit_i(size_t i) { - size_t nsp = m_thermo->nSpecies(); // Store current molefracs and set temp molefrac of species i to 1.0: - vector_fp molefracs = store(i,nsp); + vector_fp molefracs = store(i, m_thermo->nSpecies()); double zc = m_thermo->critCompressibility(); // Restore actual molefracs: @@ -446,12 +437,8 @@ vector_fp HighPressureGasTransport::store(size_t i, size_t nsp) { vector_fp molefracs(nsp); m_thermo->getMoleFractions(&molefracs[0]); - vector_fp mf_temp(nsp); - for (size_t j = 0; j < nsp; j++) { - if (j == i) { - mf_temp[j] = 1; - } else {mf_temp[j] = 0;} - } + vector_fp mf_temp(nsp, 0.0); + mf_temp[i] = 1; m_thermo->setMoleFractions(&mf_temp[0]); return molefracs; } diff --git a/src/transport/LTPspecies.cpp b/src/transport/LTPspecies.cpp index bad007607..acccc3adf 100644 --- a/src/transport/LTPspecies.cpp +++ b/src/transport/LTPspecies.cpp @@ -16,7 +16,6 @@ namespace Cantera class LTPError : public CanteraError { public: - //! Constructor is a wrapper around CanteraError /*! * @param msg Informative message @@ -44,8 +43,7 @@ static void getArrhenius(const XML_Node& node, // parse the children for the A, b, and E components. A = getFloat(node, "A", "toSI"); b = getFloat(node, "b"); - E = getFloat(node, "E", "actEnergy"); - E /= GasConstant; + E = getFloat(node, "E", "actEnergy") / GasConstant; } LTPspecies::LTPspecies(const XML_Node* const propNode, const std::string name, diff --git a/src/transport/LiquidTranInteraction.cpp b/src/transport/LiquidTranInteraction.cpp index 593910f25..60c0d2694 100644 --- a/src/transport/LiquidTranInteraction.cpp +++ b/src/transport/LiquidTranInteraction.cpp @@ -35,20 +35,16 @@ LiquidTranInteraction::LiquidTranInteraction(TransportPropertyType tp_ind) : LiquidTranInteraction::~LiquidTranInteraction() { - size_t kmax = m_Aij.size(); - for (size_t k = 0; k < kmax; k++) { + for (size_t k = 0; k < m_Aij.size(); k++) { delete m_Aij[k]; } - kmax = m_Bij.size(); - for (size_t k = 0; k < kmax; k++) { + for (size_t k = 0; k < m_Bij.size(); k++) { delete m_Bij[k]; } - kmax = m_Hij.size(); - for (size_t k = 0; k < kmax; k++) { + for (size_t k = 0; k < m_Hij.size(); k++) { delete m_Hij[k]; } - kmax = m_Sij.size(); - for (size_t k = 0; k < kmax; k++) { + for (size_t k = 0; k < m_Sij.size(); k++) { delete m_Sij[k]; } } @@ -60,19 +56,16 @@ void LiquidTranInteraction::init(const XML_Node& compModelNode, size_t nsp = thermo->nSpecies(); m_Dij.resize(nsp, nsp, 0.0); m_Eij.resize(nsp, nsp, 0.0); - std::string speciesA; - std::string speciesB; - size_t num = compModelNode.nChildren(); - for (size_t iChild = 0; iChild < num; iChild++) { + for (size_t iChild = 0; iChild < compModelNode.nChildren(); iChild++) { XML_Node& xmlChild = compModelNode.child(iChild); std::string nodeName = lowercase(xmlChild.name()); if (nodeName != "interaction") { throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", "expected element and got <" + nodeName + ">"); } - speciesA = xmlChild.attrib("speciesA"); - speciesB = xmlChild.attrib("speciesB"); + string speciesA = xmlChild.attrib("speciesA"); + string speciesB = xmlChild.attrib("speciesB"); size_t iSpecies = m_thermo->speciesIndex(speciesA); if (iSpecies == npos) { throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", @@ -417,8 +410,7 @@ doublereal LTI_Log_MoleFracs::getMixTransProp(std::vector LTPptrs) } } } - value = exp(value); - return value; + return exp(value); } void LTI_Pairwise_Interaction::setParameters(LiquidTransportParams& trParam) @@ -436,22 +428,12 @@ void LTI_Pairwise_Interaction::setParameters(LiquidTransportParams& trParam) doublereal LTI_Pairwise_Interaction::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight) { - size_t nsp = m_thermo->nSpecies(); - vector_fp molefracs(nsp); - m_thermo->getMoleFractions(&molefracs[0]); - doublereal value = 0; throw LTPmodelError("Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense."); - return value; } doublereal LTI_Pairwise_Interaction::getMixTransProp(std::vector LTPptrs) { - size_t nsp = m_thermo->nSpecies(); - vector_fp molefracs(nsp); - m_thermo->getMoleFractions(&molefracs[0]); - doublereal value = 0; throw LTPmodelError("Calling LTI_Pairwise_Interaction::getMixTransProp does not make sense."); - return value; } void LTI_Pairwise_Interaction::getMatrixTransProp(DenseMatrix& mat, doublereal* speciesValues) @@ -478,18 +460,17 @@ void LTI_Pairwise_Interaction::getMatrixTransProp(DenseMatrix& mat, doublereal* void LTI_StefanMaxwell_PPN::setParameters(LiquidTransportParams& trParam) { size_t nsp = m_thermo->nSpecies(); - size_t nsp2 = nsp*nsp; m_ionCondMix = 0; m_ionCondMixModel = trParam.ionConductivity; m_ionCondSpecies.resize(nsp,0); m_mobRatMix.resize(nsp,nsp,0.0); - m_mobRatMixModel.resize(nsp2); - m_mobRatSpecies.resize(nsp2); + m_mobRatMixModel.resize(nsp*nsp); + m_mobRatSpecies.resize(nsp*nsp); m_selfDiffMix.resize(nsp,0.0); m_selfDiffMixModel.resize(nsp); m_selfDiffSpecies.resize(nsp); - for (size_t k = 0; k < nsp2; k++) { + for (size_t k = 0; k < nsp*nsp; k++) { m_mobRatMixModel[k] = trParam.mobilityRatio[k]; m_mobRatSpecies[k].resize(nsp,0); } @@ -501,7 +482,7 @@ void LTI_StefanMaxwell_PPN::setParameters(LiquidTransportParams& trParam) for (size_t k = 0; k < nsp; k++) { LiquidTransportData& ltd = trParam.LTData[k]; m_ionCondSpecies[k] = ltd.ionConductivity; - for (size_t j = 0; j < nsp2; j++) { + for (size_t j = 0; j < nsp*nsp; j++) { m_mobRatSpecies[j][k] = ltd.mobilityRatio[j]; } for (size_t j = 0; j < nsp; j++) { @@ -512,22 +493,12 @@ void LTI_StefanMaxwell_PPN::setParameters(LiquidTransportParams& trParam) doublereal LTI_StefanMaxwell_PPN::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight) { - size_t nsp = m_thermo->nSpecies(); - vector_fp molefracs(nsp); - m_thermo->getMoleFractions(&molefracs[0]); - doublereal value = 0; throw LTPmodelError("Calling LTI_StefanMaxwell_PPN::getMixTransProp does not make sense."); - return value; } doublereal LTI_StefanMaxwell_PPN::getMixTransProp(std::vector LTPptrs) { - size_t nsp = m_thermo->nSpecies(); - vector_fp molefracs(nsp); - m_thermo->getMoleFractions(&molefracs[0]); - doublereal value = 0; throw LTPmodelError("Calling LTI_StefanMaxwell_PPN::getMixTransProp does not make sense."); - return value; } void LTI_StefanMaxwell_PPN::getMatrixTransProp(DenseMatrix& mat, doublereal* speciesValues) @@ -588,16 +559,13 @@ void LTI_StefanMaxwell_PPN::getMatrixTransProp(DenseMatrix& mat, doublereal* spe double vM = viS[anion[0]]; double zP = charges[cation[0]]; double zM = charges[anion[0]]; - doublereal xA, xB, eps; - doublereal inv_vP_vM_MutualDiff; - vector_fp dlnActCoeffdlnN_diag; - dlnActCoeffdlnN_diag.resize(neut_molefracs.size(),0.0); + vector_fp dlnActCoeffdlnN_diag(neut_molefracs.size(),0.0); marg_thermo->getdlnActCoeffdlnN_diag(&dlnActCoeffdlnN_diag[0]); - xA = neut_molefracs[neutMolIndex[cation[0]]]; - xB = neut_molefracs[neutMolIndex[cation[1]]]; - eps = (1-m_mobRatMix(cation[1],cation[0]))/(xA+xB*m_mobRatMix(cation[1],cation[0])); - inv_vP_vM_MutualDiff = (xA*(1-xB+dlnActCoeffdlnN_diag[neutMolIndex[cation[1]]])/m_selfDiffMix[cation[1]]+xB*(1-xA+dlnActCoeffdlnN_diag[neutMolIndex[cation[0]]])/m_selfDiffMix[cation[0]]); + double xA = neut_molefracs[neutMolIndex[cation[0]]]; + double xB = neut_molefracs[neutMolIndex[cation[1]]]; + double eps = (1-m_mobRatMix(cation[1],cation[0]))/(xA+xB*m_mobRatMix(cation[1],cation[0])); + double inv_vP_vM_MutualDiff = (xA*(1-xB+dlnActCoeffdlnN_diag[neutMolIndex[cation[1]]])/m_selfDiffMix[cation[1]]+xB*(1-xA+dlnActCoeffdlnN_diag[neutMolIndex[cation[0]]])/m_selfDiffMix[cation[0]]); mat.resize(nsp, nsp, 0.0); mat(cation[0],cation[1]) = mat(cation[1],cation[0]) = (1+vM/vP)*(1+eps*xB)*(1-eps*xA)*inv_vP_vM_MutualDiff-zP*zP*Faraday*Faraday/GasConstant/temp/m_ionCondMix/vol; @@ -607,22 +575,12 @@ void LTI_StefanMaxwell_PPN::getMatrixTransProp(DenseMatrix& mat, doublereal* spe doublereal LTI_StokesEinstein::getMixTransProp(doublereal* speciesValues, doublereal* speciesWeight) { - size_t nsp = m_thermo->nSpecies(); - vector_fp molefracs(nsp); - m_thermo->getMoleFractions(&molefracs[0]); - doublereal value = 0; throw LTPmodelError("Calling LTI_StokesEinstein::getMixTransProp does not make sense."); - return value; } doublereal LTI_StokesEinstein::getMixTransProp(std::vector LTPptrs) { - size_t nsp = m_thermo->nSpecies(); - vector_fp molefracs(nsp); - m_thermo->getMoleFractions(&molefracs[0]); - doublereal value = 0; throw LTPmodelError("Calling LTI_StokesEinstein::getMixTransProp does not make sense."); - return value; } void LTI_StokesEinstein::setParameters(LiquidTransportParams& trParam) diff --git a/src/transport/LiquidTransport.cpp b/src/transport/LiquidTransport.cpp index 87643b9e1..2ca59a84e 100644 --- a/src/transport/LiquidTransport.cpp +++ b/src/transport/LiquidTransport.cpp @@ -563,14 +563,8 @@ void LiquidTransport::set_Grad_X(const doublereal* grad_X) doublereal LiquidTransport::getElectricConduct() { vector_fp gradT(m_nDim,0.0); - vector_fp gradX(m_nDim * m_nsp); - vector_fp gradV(m_nDim); - for (size_t i = 0; i < m_nDim; i++) { - for (size_t k = 0; k < m_nsp; k++) { - gradX[ i*m_nDim + k] = 0.0; - } - gradV[i] = 1.0; - } + vector_fp gradX(m_nDim * m_nsp, 0.0); + vector_fp gradV(m_nDim, 1.0); set_Grad_T(&gradT[0]); set_Grad_X(&gradX[0]); @@ -887,9 +881,8 @@ void LiquidTransport::updateHydrodynamicRadius_T() void LiquidTransport::update_Grad_lnAC() { - doublereal grad_T; for (size_t k = 0; k < m_nDim; k++) { - grad_T = m_Grad_T[k]; + double grad_T = m_Grad_T[k]; size_t start = m_nsp*k; m_thermo->getdlnActCoeffds(grad_T, &m_Grad_X[start], &m_Grad_lnAC[start]); for (size_t i = 0; i < m_nsp; i++) { @@ -900,12 +893,10 @@ void LiquidTransport::update_Grad_lnAC() } } } - return; } void LiquidTransport::stefan_maxwell_solve() { - doublereal tmp; m_B.resize(m_nsp, m_nDim, 0.0); m_A.resize(m_nsp, m_nsp, 0.0); @@ -1002,7 +993,7 @@ void LiquidTransport::stefan_maxwell_solve() m_A(i,i) = 0.0; for (size_t j = 0; j < m_nsp; j++) { if (j != i) { - tmp = m_molefracs_tran[j] * m_bdiff(i,j); + double tmp = m_molefracs_tran[j] * m_bdiff(i,j); m_A(i,i) -= tmp; m_A(i,j) = tmp; } @@ -1044,7 +1035,7 @@ void LiquidTransport::stefan_maxwell_solve() m_A(i,i) = 0.0; for (size_t j = 0; j < m_nsp; j++) { if (j != i) { - tmp = m_molefracs_tran[j] * m_bdiff(i,j); + double tmp = m_molefracs_tran[j] * m_bdiff(i,j); m_A(i,i) -= tmp; m_A(i,j) = tmp; } @@ -1084,7 +1075,7 @@ void LiquidTransport::stefan_maxwell_solve() m_A(i,i) = 0.0; for (size_t j = 0; j < m_nsp; j++) { if (j != i) { - tmp = m_molefracs_tran[j] * m_bdiff(i,j); + double tmp = m_molefracs_tran[j] * m_bdiff(i,j); m_A(i,i) -= tmp; m_A(i,j) = tmp; } diff --git a/src/transport/LiquidTransportData.cpp b/src/transport/LiquidTransportData.cpp index 85000a484..97793db51 100644 --- a/src/transport/LiquidTransportData.cpp +++ b/src/transport/LiquidTransportData.cpp @@ -80,14 +80,10 @@ LiquidTransportData::~LiquidTransportData() delete ionConductivity; for (size_t k = 0; k < mobilityRatio.size(); k++) { - if (mobilityRatio[k]) { - delete mobilityRatio[k]; - } + delete mobilityRatio[k]; } for (size_t k = 0; k < selfDiffusion.size(); k++) { - if (selfDiffusion[k]) { - delete selfDiffusion[k]; - } + delete selfDiffusion[k]; } delete thermalCond; diff --git a/src/transport/MMCollisionInt.cpp b/src/transport/MMCollisionInt.cpp index ce9b79428..5d233ec04 100644 --- a/src/transport/MMCollisionInt.cpp +++ b/src/transport/MMCollisionInt.cpp @@ -20,13 +20,12 @@ double MMCollisionInt::delta[8] = {0.0, 0.25, 0.50, 0.75, 1.0, doublereal quadInterp(doublereal x0, doublereal* x, doublereal* y) { - doublereal dx21, dx32, dx31, dy32, dy21, a; - dx21 = x[1] - x[0]; - dx32 = x[2] - x[1]; - dx31 = dx21 + dx32; - dy32 = y[2] - y[1]; - dy21 = y[1] - y[0]; - a = (dx21*dy32 - dy21*dx32)/(dx21*dx31*dx32); + double dx21 = x[1] - x[0]; + double dx32 = x[2] - x[1]; + double dx31 = dx21 + dx32; + double dy32 = y[2] - y[1]; + double dy21 = y[1] - y[0]; + double a = (dx21*dy32 - dy21*dx32)/(dx21*dx31*dx32); return a*(x0 - x[0])*(x0 - x[1]) + (dy21/dx21)*(x0 - x[1]) + y[1]; } @@ -237,7 +236,7 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level) writelogf("T*_max = %g\n", tstar[m_nmax + 1]); } m_logTemp.resize(37); - doublereal rmserr, e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0; + doublereal e22 = 0.0, ea = 0.0, eb = 0.0, ec = 0.0; if (m_loglevel > 0) { writelog("Collision integral fits at each tabulated T* vs. delta*.\n" @@ -253,7 +252,7 @@ void MMCollisionInt::init(doublereal tsmin, doublereal tsmax, int log_level) m_logTemp[i] = log(tstar[i+1]); vector_fp c(DeltaDegree+1); - rmserr = fitDelta(0, i, DeltaDegree, c.data()); + double rmserr = fitDelta(0, i, DeltaDegree, c.data()); if (log_level > 3) { writelogf("\ndelta* fit at T* = %.6g\n", tstar[i+1]); writelog("omega22 = [" + vec2str(c) + "]\n"); @@ -319,12 +318,13 @@ doublereal MMCollisionInt::fitDelta(int table, int ntstar, int degree, doublerea doublereal MMCollisionInt::omega22(double ts, double deltastar) { int i; - for (i = 0; i < 37; i++) if (ts < tstar22[i]) { + for (i = 0; i < 37; i++) { + if (ts < tstar22[i]) { break; } - int i1, i2; - i1 = std::max(i - 1, 0); - i2 = i1+3; + } + int i1 = std::max(i - 1, 0); + int i2 = i1+3; if (i2 > 36) { i2 = 36; i1 = i2 - 3; @@ -346,9 +346,8 @@ doublereal MMCollisionInt::astar(double ts, double deltastar) for (i = 0; i < 37; i++) if (ts < tstar22[i]) { break; } - int i1, i2; - i1 = std::max(i - 1, 0); - i2 = i1+3; + int i1 = std::max(i - 1, 0); + int i2 = i1+3; if (i2 > 36) { i2 = 36; i1 = i2 - 3; @@ -370,9 +369,8 @@ doublereal MMCollisionInt::bstar(double ts, double deltastar) for (i = 0; i < 37; i++) if (ts < tstar22[i]) { break; } - int i1, i2; - i1 = std::max(i - 1, 0); - i2 = i1+3; + int i1 = std::max(i - 1, 0); + int i2 = i1+3; if (i2 > 36) { i2 = 36; i1 = i2 - 3; @@ -394,9 +392,8 @@ doublereal MMCollisionInt::cstar(double ts, double deltastar) for (i = 0; i < 37; i++) if (ts < tstar22[i]) { break; } - int i1, i2; - i1 = std::max(i - 1,0); - i2 = i1+3; + int i1 = std::max(i - 1,0); + int i2 = i1+3; if (i2 > 36) { i2 = 36; i1 = i2 - 3; @@ -415,12 +412,11 @@ doublereal MMCollisionInt::cstar(double ts, double deltastar) void MMCollisionInt::fit_omega22(int degree, doublereal deltastar, doublereal* o22) { - int i, n = m_nmax - m_nmin + 1; + int n = m_nmax - m_nmin + 1; vector_fp values(n); - doublereal rmserr; vector_fp w(n); doublereal* logT = &m_logTemp[m_nmin]; - for (i = 0; i < n; i++) { + for (int i = 0; i < n; i++) { if (deltastar == 0.0) { values[i] = omega22_table[8*(i + m_nmin)]; } else { @@ -428,7 +424,7 @@ void MMCollisionInt::fit_omega22(int degree, doublereal deltastar, } } w[0]= -1.0; - rmserr = polyfit(n, degree, logT, values.data(), w.data(), o22); + double rmserr = polyfit(n, degree, logT, values.data(), w.data(), o22); if (m_loglevel > 0 && rmserr > 0.01) { writelogf("Warning: RMS error = %12.6g in omega_22 fit" "with delta* = %12.6g\n", rmserr, deltastar); @@ -438,12 +434,11 @@ void MMCollisionInt::fit_omega22(int degree, doublereal deltastar, void MMCollisionInt::fit(int degree, doublereal deltastar, doublereal* a, doublereal* b, doublereal* c) { - int i, n = m_nmax - m_nmin + 1; + int n = m_nmax - m_nmin + 1; vector_fp values(n); - doublereal rmserr; vector_fp w(n); doublereal* logT = &m_logTemp[m_nmin]; - for (i = 0; i < n; i++) { + for (int i = 0; i < n; i++) { if (deltastar == 0.0) { values[i] = astar_table[8*(i + m_nmin + 1)]; } else { @@ -451,9 +446,9 @@ void MMCollisionInt::fit(int degree, doublereal deltastar, } } w[0]= -1.0; - rmserr = polyfit(n, degree, logT, values.data(), w.data(), a); + double rmserr = polyfit(n, degree, logT, values.data(), w.data(), a); - for (i = 0; i < n; i++) { + for (int i = 0; i < n; i++) { if (deltastar == 0.0) { values[i] = bstar_table[8*(i + m_nmin + 1)]; } else { @@ -463,7 +458,7 @@ void MMCollisionInt::fit(int degree, doublereal deltastar, w[0]= -1.0; rmserr = polyfit(n, degree, logT, values.data(), w.data(), b); - for (i = 0; i < n; i++) { + for (int i = 0; i < n; i++) { if (deltastar == 0.0) { values[i] = cstar_table[8*(i + m_nmin + 1)]; } else { diff --git a/src/transport/MultiTransport.cpp b/src/transport/MultiTransport.cpp index 604cbf046..bbc922342 100644 --- a/src/transport/MultiTransport.cpp +++ b/src/transport/MultiTransport.cpp @@ -219,23 +219,20 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_ } // copy grad_X to fluxes - const doublereal* gx; for (size_t n = 0; n < ndim; n++) { - gx = grad_X + ldx*n; + const double* gx = grad_X + ldx*n; copy(gx, gx + m_nsp, fluxes + ldf*n); fluxes[jmax + n*ldf] = 0.0; } // solve the equations solve(m_aa, fluxes, ndim, ldf); - - size_t offset; doublereal pp = pressure_ig(); // multiply diffusion velocities by rho * V to create mass fluxes, and // restore the gradx elements that were modified for (size_t n = 0; n < ndim; n++) { - offset = n*ldf; + size_t offset = n*ldf; for (size_t i = 0; i < m_nsp; i++) { fluxes[i + offset] *= rho * y[i] / pp; } @@ -244,7 +241,7 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_ // thermal diffusion if (addThermalDiffusion) { for (size_t n = 0; n < ndim; n++) { - offset = n*ldf; + size_t offset = n*ldf; doublereal grad_logt = grad_T[n]/m_temp; for (size_t i = 0; i < m_nsp; i++) { fluxes[i + offset] -= m_spwork[i]*grad_logt; @@ -259,7 +256,7 @@ void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* s double* x1 = m_spwork1.data(); double* x2 = m_spwork2.data(); double* x3 = m_spwork3.data(); - size_t n, nsp = m_thermo->nSpecies(); + size_t nsp = m_thermo->nSpecies(); m_thermo->restoreState(nsp+2, state1); double p1 = m_thermo->pressure(); double t1 = state1[0]; @@ -273,7 +270,7 @@ void MultiTransport::getMassFluxes(const doublereal* state1, const doublereal* s double p = 0.5*(p1 + p2); double t = 0.5*(state1[0] + state2[0]); - for (n = 0; n < nsp; n++) { + for (size_t n = 0; n < nsp; n++) { x3[n] = 0.5*(x1[n] + x2[n]); } m_thermo->setState_TPX(t, p, x3); @@ -380,10 +377,9 @@ void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d) doublereal prefactor = 16.0 * m_temp * m_thermo->meanMolecularWeight()/(25.0 * p); - doublereal c; for (size_t i = 0; i < m_nsp; i++) { for (size_t j = 0; j < m_nsp; j++) { - c = prefactor/m_mw[j]; + double c = prefactor/m_mw[j]; d[ld*j + i] = c*m_molefracs[i]* (m_Lmatrix(i,j) - m_Lmatrix(i,i)); } @@ -430,12 +426,10 @@ void MultiTransport::updateThermal_T() updateDiff_T(); // evaluate polynomial fits for A*, B*, C* - doublereal z; - int ipoly; for (size_t i = 0; i < m_nsp; i++) { for (size_t j = i; j < m_nsp; j++) { - z = m_logt - m_log_eps_k(i,j); - ipoly = m_poly[i][j]; + double z = m_logt - m_log_eps_k(i,j); + int ipoly = m_poly[i][j]; if (m_mode == CK_Mode) { m_om22(i,j) = poly6(z, m_omega22_poly[ipoly].data()); m_astar(i,j) = poly6(z, m_astar_poly[ipoly].data()); @@ -456,18 +450,15 @@ void MultiTransport::updateThermal_T() m_abc_ok = true; // evaluate the temperature-dependent rotational relaxation rate - doublereal tr, sqtr; for (size_t k = 0; k < m_nsp; k++) { - tr = m_eps[k]/ m_kbt; - sqtr = m_sqrt_eps_k[k] / m_sqrt_t; + double tr = m_eps[k]/ m_kbt; + double sqtr = m_sqrt_eps_k[k] / m_sqrt_t; m_rotrelax[k] = std::max(1.0,m_zrot[k]) * m_frot_298[k]/Frot(tr, sqtr); } - doublereal d; doublereal c = 1.2*GasConstant*m_temp; for (size_t k = 0; k < m_nsp; k++) { - d = c * m_visc[k] * m_astar(k,k)/m_mw[k]; - m_bdiff(k,k) = d; + m_bdiff(k,k) = c * m_visc[k] * m_astar(k,k)/m_mw[k]; } // Calculate the internal heat capacities by subtracting off the translational contributions @@ -520,11 +511,10 @@ void MultiTransport::eval_L0000(const doublereal* const x) void MultiTransport::eval_L0010(const doublereal* const x) { doublereal prefactor = 1.6*m_temp; - doublereal sum, wj, xj; for (size_t j = 0; j < m_nsp; j++) { - xj = x[j]; - wj = m_mw[j]; - sum = 0.0; + double xj = x[j]; + double wj = m_mw[j]; + double sum = 0.0; for (size_t 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) / @@ -551,26 +541,22 @@ void MultiTransport::eval_L1010(const doublereal* x) { const doublereal fiveover3pi = 5.0/(3.0*Pi); doublereal prefactor = (16.0*m_temp)/25.0; - doublereal constant1, wjsq, constant2, constant3, constant4, - fourmj, threemjsq, sum, sumwij;; - doublereal term1, term2; for (size_t j = 0; j < m_nsp; 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; + double constant1 = prefactor*x[j]; + double wjsq = m_mw[j]*m_mw[j]; + double constant2 = 13.75*wjsq; + double constant3 = m_crot[j]/m_rotrelax[j]; + double constant4 = 7.5*wjsq; + double fourmj = 4.0*m_mw[j]; + double threemjsq = 3.0*m_mw[j]*m_mw[j]; + double sum = 0.0; for (size_t 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) + double sumwij = m_mw[i] + m_mw[j]; + double term1 = m_bdiff(i,j) * sumwij*sumwij; + double 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) @@ -588,25 +574,21 @@ void MultiTransport::eval_L1010(const doublereal* x) void MultiTransport::eval_L1001(const doublereal* x) { doublereal prefactor = 32.00*m_temp/(5.00*Pi); - doublereal constant, sum; - size_t n2 = 2*m_nsp; - int npoly = 0; for (size_t 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; + double constant = prefactor*m_mw[j]*x[j]*m_crot[j]/(m_cinternal[j]*m_rotrelax[j]); + double sum = 0.0; for (size_t 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_Lmatrix(i+m_nsp,j+2*m_nsp) = 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); + sum += m_Lmatrix(i+m_nsp,j+2*m_nsp); } - npoly++; - m_Lmatrix(j+m_nsp,j+n2) += sum; + m_Lmatrix(j+m_nsp,j+2*m_nsp) += sum; } else { for (size_t i = 0; i < m_nsp; i++) { - m_Lmatrix(i+m_nsp,j+n2) = 0.0; + m_Lmatrix(i+m_nsp,j+2*m_nsp) = 0.0; } } } @@ -614,64 +596,57 @@ void MultiTransport::eval_L1001(const doublereal* x) void MultiTransport::eval_L0001() { - size_t n2 = 2*m_nsp; for (size_t j = 0; j < m_nsp; j++) { for (size_t i = 0; i < m_nsp; i++) { - m_Lmatrix(i,j+n2) = 0.0; + m_Lmatrix(i,j+2*m_nsp) = 0.0; } } } void MultiTransport::eval_L0100() { - size_t n2 = 2*m_nsp; for (size_t j = 0; j < m_nsp; j++) { for (size_t i = 0; i < m_nsp; i++) { - m_Lmatrix(i+n2,j) = 0.0; // see Eq. (12.123) + m_Lmatrix(i+2*m_nsp,j) = 0.0; // see Eq. (12.123) } } } void MultiTransport::eval_L0110() { - size_t n2 = 2*m_nsp; for (size_t j = 0; j < m_nsp; j++) { for (size_t i = 0; i < m_nsp; i++) { - m_Lmatrix(i+n2,j+m_nsp) = m_Lmatrix(j+m_nsp,i+n2); // see Eq. (12.123) + m_Lmatrix(i+2*m_nsp,j+m_nsp) = m_Lmatrix(j+m_nsp,i+2*m_nsp); // 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; - size_t n2 = 2*m_nsp; - doublereal constant1, constant2, diff_int, sum; for (size_t 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; + double constant1 = 4*m_temp*x[i]/m_cinternal[i]; + double constant2 = 12*m_mw[i]*m_crot[i] / + (5*Pi*m_cinternal[i]*m_rotrelax[i]); + double sum = 0.0; for (size_t k = 0; k < m_nsp; k++) { // see Eq. (12.131) - diff_int = m_bdiff(i,k); - m_Lmatrix(k+n2,i+n2) = 0.0; + double diff_int = m_bdiff(i,k); + m_Lmatrix(k+2*m_nsp,i+2*m_nsp) = 0.0; sum += x[k]/diff_int; - if (k != i) sum += x[k]*m_astar(i,k)*constant2 / - (m_mw[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_Lmatrix(i+2*m_nsp,i+2*m_nsp) = + - 8/Pi*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 (size_t k = 0; k < m_nsp; k++) { - m_Lmatrix(i+n2,i+n2) = 1.0; + m_Lmatrix(i+2*m_nsp,i+2*m_nsp) = 1.0; } } } diff --git a/src/transport/SimpleTransport.cpp b/src/transport/SimpleTransport.cpp index 7ca1a1691..c766486ba 100644 --- a/src/transport/SimpleTransport.cpp +++ b/src/transport/SimpleTransport.cpp @@ -176,9 +176,7 @@ bool SimpleTransport::initLiquid(LiquidTransportParams& tr) m_viscSpecies.resize(m_nsp); m_coeffVisc_Ns.clear(); m_coeffVisc_Ns.resize(m_nsp); - std::string spName = m_thermo->speciesName(0); for (size_t k = 0; k < m_nsp; k++) { - spName = m_thermo->speciesName(k); LiquidTransportData& ltd = tr.LTData[k]; m_coeffVisc_Ns[k] = ltd.viscosity; ltd.viscosity = 0; @@ -189,7 +187,6 @@ bool SimpleTransport::initLiquid(LiquidTransportParams& tr) m_coeffLambda_Ns.clear(); m_coeffLambda_Ns.resize(m_nsp); for (size_t k = 0; k < m_nsp; k++) { - spName = m_thermo->speciesName(k); LiquidTransportData& ltd = tr.LTData[k]; m_coeffLambda_Ns[k] = ltd.thermalCond; ltd.thermalCond = 0; @@ -201,7 +198,7 @@ bool SimpleTransport::initLiquid(LiquidTransportParams& tr) m_coeffDiff_Ns.clear(); m_coeffDiff_Ns.resize(m_nsp); for (size_t k = 0; k < m_nsp; k++) { - spName = m_thermo->speciesName(k); + string spName = m_thermo->speciesName(k); LiquidTransportData& ltd = tr.LTData[k]; m_coeffDiff_Ns[k] = ltd.speciesDiffusivity; ltd.speciesDiffusivity = 0; @@ -280,7 +277,6 @@ void SimpleTransport::getSpeciesViscosities(doublereal* const visc) void SimpleTransport::getBinaryDiffCoeffs(size_t ld, doublereal* d) { - double bdiff; update_T(); // if necessary, evaluate the species diffusion coefficients @@ -291,8 +287,7 @@ void SimpleTransport::getBinaryDiffCoeffs(size_t ld, doublereal* d) for (size_t i = 0; i < m_nsp; i++) { for (size_t j = 0; j < m_nsp; j++) { - bdiff = 0.5 * (m_diffSpecies[i] + m_diffSpecies[j]); - d[i*m_nsp+j] = bdiff; + d[i*m_nsp+j] = 0.5 * (m_diffSpecies[i] + m_diffSpecies[j]); } } } diff --git a/src/transport/SolidTransport.cpp b/src/transport/SolidTransport.cpp index 26d670838..371fb17d8 100644 --- a/src/transport/SolidTransport.cpp +++ b/src/transport/SolidTransport.cpp @@ -160,8 +160,7 @@ void SolidTransport::getMobilities(doublereal* const mobil) void SolidTransport::getMixDiffCoeffs(doublereal* const d) { - size_t nsp = m_thermo->nSpecies(); - for (size_t k = 0; k < nsp; k++) { + for (size_t k = 0; k < m_thermo->nSpecies(); k++) { d[k] = 0.0; } } diff --git a/src/transport/TransportFactory.cpp b/src/transport/TransportFactory.cpp index df4bf7900..066b51cb6 100644 --- a/src/transport/TransportFactory.cpp +++ b/src/transport/TransportFactory.cpp @@ -100,26 +100,19 @@ std::string TransportFactory::modelName(int model) LTPspecies* TransportFactory::newLTP(const XML_Node& trNode, const std::string& name, TransportPropertyType tp_ind, thermo_t* thermo) { - LTPspecies* ltps = 0; std::string model = lowercase(trNode["model"]); switch (m_LTRmodelMap[model]) { case LTP_TD_CONSTANT: - ltps = new LTPspecies_Const(trNode, name, tp_ind, thermo); - break; + return new LTPspecies_Const(trNode, name, tp_ind, thermo); case LTP_TD_ARRHENIUS: - ltps = new LTPspecies_Arrhenius(trNode, name, tp_ind, thermo); - break; + return new LTPspecies_Arrhenius(trNode, name, tp_ind, thermo); case LTP_TD_POLY: - ltps = new LTPspecies_Poly(trNode, name, tp_ind, thermo); - break; + return new LTPspecies_Poly(trNode, name, tp_ind, thermo); case LTP_TD_EXPT: - ltps = new LTPspecies_ExpT(trNode, name, tp_ind, thermo); - break; + return new LTPspecies_ExpT(trNode, name, tp_ind, thermo); default: throw CanteraError("TransportFactory::newLTP","unknown transport model: " + model); - ltps = new LTPspecies(&trNode, name, tp_ind, thermo); } - return ltps; } LiquidTranInteraction* TransportFactory::newLTI(const XML_Node& trNode, @@ -127,54 +120,52 @@ LiquidTranInteraction* TransportFactory::newLTI(const XML_Node& trNode, LiquidTransportParams& trParam) { LiquidTranInteraction* lti = 0; - thermo_t* thermo = trParam.thermo; - std::string model = trNode["model"]; - switch (m_LTImodelMap[model]) { + switch (m_LTImodelMap[trNode["model"]]) { case LTI_MODEL_SOLVENT: lti = new LTI_Solvent(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); break; case LTI_MODEL_MOLEFRACS: lti = new LTI_MoleFracs(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); break; case LTI_MODEL_MASSFRACS: lti = new LTI_MassFracs(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); break; case LTI_MODEL_LOG_MOLEFRACS: lti = new LTI_Log_MoleFracs(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); break; case LTI_MODEL_PAIRWISE_INTERACTION: lti = new LTI_Pairwise_Interaction(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); lti->setParameters(trParam); break; case LTI_MODEL_STEFANMAXWELL_PPN: lti = new LTI_StefanMaxwell_PPN(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); lti->setParameters(trParam); break; case LTI_MODEL_STOKES_EINSTEIN: lti = new LTI_StokesEinstein(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); lti->setParameters(trParam); break; case LTI_MODEL_MOLEFRACS_EXPT: lti = new LTI_MoleFracs_ExpT(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); break; case LTI_MODEL_NOTSET: case LTI_MODEL_NONE: case LTI_MODEL_MULTIPLE: lti = new LiquidTranInteraction(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); break; default: // @TODO make sure we can throw an error here with existing datasets and tests before changing code lti = new LiquidTranInteraction(tp_ind); - lti->init(trNode, thermo); + lti->init(trNode, trParam.thermo); } return lti; } @@ -269,9 +260,7 @@ void TransportFactory::setupLiquidTransport(thermo_t* thermo, int log_level, trParam.log_level = log_level; // Get the molecular weights and load them into trParam - trParam.mw.resize(nsp); - copy(trParam.thermo->molecularWeights().begin(), - trParam.thermo->molecularWeights().end(), trParam.mw.begin()); + trParam.mw = trParam.thermo->molecularWeights(); // Resize all other vectors in trParam trParam.LTData.resize(nsp); @@ -282,7 +271,7 @@ void TransportFactory::setupLiquidTransport(thermo_t* thermo, int log_level, trParam.diff_Dij.resize(nsp,nsp); trParam.radius_Aij.resize(nsp,nsp); - XML_Node root, log; + XML_Node log; // Note that getLiquidSpeciesTransportData just populates the pure species transport data. getLiquidSpeciesTransportData(species_database, log, trParam.thermo->speciesNames(), trParam); @@ -302,21 +291,17 @@ void TransportFactory::setupSolidTransport(thermo_t* thermo, int log_level, // constant mixture attributes trParam.thermo = thermo; trParam.nsp_ = trParam.thermo->nSpecies(); - size_t nsp = trParam.nsp_; trParam.tmin = thermo->minTemp(); trParam.tmax = thermo->maxTemp(); trParam.log_level = log_level; // Get the molecular weights and load them into trParam - trParam.mw.resize(nsp); - copy(trParam.thermo->molecularWeights().begin(), - trParam.thermo->molecularWeights().end(), trParam.mw.begin()); - - XML_Node root, log; + trParam.mw = trParam.thermo->molecularWeights(); // getSolidTransportData() populates the phase transport models like // electronic conductivity thermal conductivity, interstitial diffusion if (phase_database->hasChild("transport")) { + XML_Node log; XML_Node& transportNode = phase_database->child("transport"); getSolidTransportData(transportNode, log, thermo->name(), trParam); } @@ -347,8 +332,6 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vector &names, LiquidTransportParams& trParam) { - std::string name; - // Create a map of species names versus liquid transport data parameters std::map datatable; @@ -363,8 +346,7 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vectorspeciesIndex(firstSpec)+nsp*temp_thermo->speciesIndex(secondSpec); - data.mobilityRatio[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], temp_thermo); + size_t index = trParam.thermo->speciesIndex(firstSpec)+nsp*trParam.thermo->speciesIndex(secondSpec); + data.mobilityRatio[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], trParam.thermo); }; }; break; @@ -406,8 +387,8 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vectorspeciesIndex(specName); - data.selfDiffusion[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], temp_thermo); + size_t index = trParam.thermo->speciesIndex(specName); + data.selfDiffusion[index] = newLTP(propSpecNode, name, m_tranPropMap[nodeName], trParam.thermo); }; }; break; @@ -415,25 +396,25 @@ void TransportFactory::getLiquidSpeciesTransportData(const std::vectorsecond; // Now, transfer these objects into LTData in the correct phase index // order by calling the default copy constructor for // LiquidTransportData. - trParam.LTData.push_back(trdat); + trParam.LTData.push_back(it->second); } } @@ -475,14 +455,12 @@ void TransportFactory::getLiquidInteractionsTransportData(const XML_Node& transp try { size_t nsp = trParam.nsp_; size_t nBinInt = nsp*(nsp-1)/2; - size_t num = transportNode.nChildren(); - for (size_t iChild = 0; iChild < num; iChild++) { + for (size_t iChild = 0; iChild < transportNode.nChildren(); iChild++) { //tranTypeNode is a type of transport property like viscosity XML_Node& tranTypeNode = transportNode.child(iChild); std::string nodeName = tranTypeNode.name(); trParam.mobilityRatio.resize(nsp*nsp,0); trParam.selfDiffusion.resize(nsp,0); - ThermoPhase* temp_thermo = trParam.thermo; if (tranTypeNode.name() == "compositionDependence") { std::string modelName = tranTypeNode.attrib("model"); @@ -514,7 +492,7 @@ void TransportFactory::getLiquidInteractionsTransportData(const XML_Node& transp size_t loc = specName.find(":"); string firstSpec = specName.substr(0,loc); string secondSpec = specName.substr(loc+1); - size_t index = temp_thermo->speciesIndex(firstSpec)+nsp*temp_thermo->speciesIndex(secondSpec); + size_t index = trParam.thermo->speciesIndex(firstSpec)+nsp*trParam.thermo->speciesIndex(secondSpec); trParam.mobilityRatio[index] = newLTI(propSpecNode, m_tranPropMap[nodeName], trParam); @@ -525,7 +503,7 @@ void TransportFactory::getLiquidInteractionsTransportData(const XML_Node& transp for (size_t iSpec = 0; iSpec< nsp; iSpec++) { XML_Node& propSpecNode = compDepNode.child(iSpec); string specName = propSpecNode.name(); - size_t index = temp_thermo->speciesIndex(specName); + size_t index = trParam.thermo->speciesIndex(specName); trParam.selfDiffusion[index] = newLTI(propSpecNode, m_tranPropMap[nodeName], trParam); @@ -591,39 +569,37 @@ void TransportFactory::getSolidTransportData(const XML_Node& transportNode, const std::string phaseName, SolidTransportData& trParam) { - size_t num = transportNode.nChildren(); - for (size_t iChild = 0; iChild < num; iChild++) { + for (size_t iChild = 0; iChild < transportNode.nChildren(); iChild++) { //tranTypeNode is a type of transport property like viscosity XML_Node& tranTypeNode = transportNode.child(iChild); std::string nodeName = tranTypeNode.name(); - ThermoPhase* temp_thermo = trParam.thermo; //tranTypeNode contains the interaction model switch (m_tranPropMap[nodeName]) { case TP_IONCONDUCTIVITY: trParam.ionConductivity = newLTP(tranTypeNode, phaseName, m_tranPropMap[nodeName], - temp_thermo); + trParam.thermo); break; case TP_THERMALCOND: trParam.thermalConductivity = newLTP(tranTypeNode, phaseName, m_tranPropMap[nodeName], - temp_thermo); + trParam.thermo); break; case TP_DEFECTDIFF: trParam.defectDiffusivity = newLTP(tranTypeNode, phaseName, m_tranPropMap[nodeName], - temp_thermo); + trParam.thermo); break; case TP_DEFECTCONC: trParam.defectActivity = newLTP(tranTypeNode, phaseName, m_tranPropMap[nodeName], - temp_thermo); + trParam.thermo); break; case TP_ELECTCOND: trParam.electConductivity = newLTP(tranTypeNode, phaseName, m_tranPropMap[nodeName], - temp_thermo); + trParam.thermo); break; default: throw CanteraError("getSolidTransportData","unknown transport property: " + nodeName);