From 01f915d01a7ad993beb1d4e40bf319f049b75aae Mon Sep 17 00:00:00 2001 From: Christopher Lueth Date: Thu, 4 Mar 2010 18:02:58 +0000 Subject: [PATCH] Changed storage of stefan maxwell coefficients from Dij to 1/Dij...to allow for negative and zero values of the friction between species. Corrected some memory issues (LTP.cpp was creating null pointers for ionConductivity that were being used in LT.cpp) --- Cantera/src/transport/LiquidTransport.cpp | 48 ++++++++++--------- .../src/transport/LiquidTransportParams.cpp | 22 ++++----- 2 files changed, 37 insertions(+), 33 deletions(-) diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp index d6303dc28..499672a64 100644 --- a/Cantera/src/transport/LiquidTransport.cpp +++ b/Cantera/src/transport/LiquidTransport.cpp @@ -802,11 +802,15 @@ namespace Cantera { // from the polynomial fits if (!m_diff_temp_ok) updateDiff_T(); - for (i = 0; i < m_nsp; i++) - for (j = 0; j < m_nsp; j++) { - d[ld*j + i] = m_bdiff(i,j); - + for (i = 0; i < m_nsp; i++) { + for (j = 0; j < m_nsp; j++){ + //if (!( ( m_bdiff(i,j) > 0.0 ) | ( m_bdiff(i,j) < 0.0 ))){ + // throw CanteraError("LiquidTransport::getBinaryDiffCoeffs ", + // "m_bdiff has zero entry in non-diagonal.");} + d[ld*j + i] = 1.0 / m_bdiff(i,j); + } + } } @@ -1759,10 +1763,10 @@ namespace Cantera { m_A(i,i) = 0.0; for (j = 0; j < m_nsp; j++){ if (j != i) { - if ( !( m_bdiff(i,j) > 0.0 ) ) - throw CanteraError("LiquidTransport::stefan_maxwell_solve", - "m_bdiff has zero entry in non-diagonal."); - tmp = m_molefracs_tran[j] / m_bdiff(i,j); + //if ( !( m_bdiff(i,j) > 0.0 ) ) + //throw CanteraError("LiquidTransport::stefan_maxwell_solve", + // "m_bdiff has zero entry in non-diagonal."); + tmp = m_molefracs_tran[j] * m_bdiff(i,j); m_A(i,i) -= tmp; m_A(i,j) = + tmp; } @@ -1773,12 +1777,12 @@ namespace Cantera { solve(m_A, m_B); /* - condSum2 = m_chargeSpecies[1]*m_chargeSpecies[1]*m_molefracs_tran[1]/m_bdiff(2,3) + - m_chargeSpecies[2]*m_chargeSpecies[2]*m_molefracs_tran[2]/m_bdiff(1,3) + - m_chargeSpecies[3]*m_chargeSpecies[3]*m_molefracs_tran[3]/m_bdiff(1,2); - condSum1 = m_molefracs_tran[1]/m_bdiff(1,2)/m_bdiff(1,3) + - m_molefracs_tran[2]/m_bdiff(2,3)/m_bdiff(1,2) + - m_molefracs_tran[3]/m_bdiff(1,3)/m_bdiff(2,3); + condSum2 = m_chargeSpecies[1]*m_chargeSpecies[1]*m_molefracs_tran[1]*m_bdiff(2,3) + + m_chargeSpecies[2]*m_chargeSpecies[2]*m_molefracs_tran[2]*m_bdiff(1,3) + + m_chargeSpecies[3]*m_chargeSpecies[3]*m_molefracs_tran[3]*m_bdiff(1,2); + condSum1 = m_molefracs_tran[1]*m_bdiff(1,2)*m_bdiff(1,3) + + m_molefracs_tran[2]*m_bdiff(2,3)*m_bdiff(1,2) + + m_molefracs_tran[3]*m_bdiff(1,3)*m_bdiff(2,3); condSum2 = condSum2/condSum1*Faraday*Faraday/GasConstant/T/vol; */ @@ -1820,10 +1824,10 @@ namespace Cantera { m_A(i,i) = 0.0; for (j = 0; j < m_nsp; j++) { if (j != i) { - if ( !( m_bdiff(i,j) > 0.0 ) ) - throw CanteraError("LiquidTransport::stefan_maxwell_solve", - "m_bdiff has zero entry in non-diagonal."); - tmp = m_molefracs_tran[j] / m_bdiff(i,j); + //if ( !( m_bdiff(i,j) > 0.0 ) ) + //throw CanteraError("LiquidTransport::stefan_maxwell_solve", + // "m_bdiff has zero entry in non-diagonal."); + tmp = m_molefracs_tran[j] * m_bdiff(i,j); m_A(i,i) -= tmp; m_A(i,j) = + tmp; } @@ -1862,10 +1866,10 @@ namespace Cantera { m_A(i,i) = 0.0; for (j = 0; j < m_nsp; j++) { if (j != i) { - if ( !( m_bdiff(i,j) > 0.0 ) ) - throw CanteraError("LiquidTransport::stefan_maxwell_solve", - "m_bdiff has zero entry in non-diagonal."); - tmp = m_molefracs_tran[j] / m_bdiff(i,j); + //if ( !( m_bdiff(i,j) > 0.0 ) ) + //throw CanteraError("LiquidTransport::stefan_maxwell_solve", + // "m_bdiff has zero entry in non-diagonal."); + tmp = m_molefracs_tran[j] * m_bdiff(i,j); m_A(i,i) -= tmp; m_A(i,j) = + tmp; } diff --git a/Cantera/src/transport/LiquidTransportParams.cpp b/Cantera/src/transport/LiquidTransportParams.cpp index 446639def..c92e913f5 100644 --- a/Cantera/src/transport/LiquidTransportParams.cpp +++ b/Cantera/src/transport/LiquidTransportParams.cpp @@ -602,11 +602,11 @@ namespace Cantera { mat.resize( nsp, nsp, 0.0 ); for ( int i = 0; i < nsp; i++ ) for ( int j = 0; j < i; j++ ) - mat(i,j) = mat(j,i) = m_Dij(i,j) * exp( - m_Eij(i,j) / temp ); + mat(i,j) = mat(j,i) = exp( m_Eij(i,j) / temp ) / m_Dij(i,j); for ( int i = 0; i < nsp; i++ ) if ( mat(i,i) == 0.0 && m_diagonals[i] ) - mat(i,i) = m_diagonals[i]->getSpeciesTransProp() ; + mat(i,i) = 1.0 / m_diagonals[i]->getSpeciesTransProp() ; } @@ -617,7 +617,7 @@ namespace Cantera { m_ionCondMix = 0; m_ionCondMixModel = trParam.ionConductivity; - trParam.ionConductivity = 0; + //trParam.ionConductivity = 0; m_ionCondSpecies.resize(nsp,0); m_mobRatMix.resize(nsp,nsp,0.0); m_mobRatMixModel.resize(nBinInt); @@ -630,13 +630,13 @@ namespace Cantera { for ( int k = 0; k < nBinInt; k++ ) { m_mobRatMixModel[k] = trParam.mobilityRatio[k]; - trParam.mobilityRatio[k] = 0; + //trParam.mobilityRatio[k] = 0; m_mobRatSpecies[k].resize(nsp,0); m_mobRatIndex[k] = trParam.mobRatIndex[k]; } for ( int k = 0; k < nsp; k++ ) { m_selfDiffMixModel[k] = trParam.selfDiffusion[k]; - trParam.selfDiffusion[k] = 0; + //trParam.selfDiffusion[k] = 0; m_selfDiffSpecies[k].resize(nsp,0); m_selfDiffIndex[k] = trParam.selfDiffIndex[k]; } @@ -644,14 +644,14 @@ namespace Cantera { for (int k = 0; k < nsp; k++) { Cantera::LiquidTransportData <d = trParam.LTData[k]; m_ionCondSpecies[k] = ltd.ionConductivity; - ltd.ionConductivity = 0; + //ltd.ionConductivity = 0; for ( int j = 0; j < nBinInt; j++ ){ m_mobRatSpecies[j][k] = ltd.mobilityRatio[j]; - ltd.mobilityRatio[j] = 0; + //ltd.mobilityRatio[j] = 0; } for ( int j = 0; j < nsp; j++ ){ m_selfDiffSpecies[j][k] = ltd.selfDiffusion[j]; - ltd.selfDiffusion[j] = 0; + //ltd.selfDiffusion[j] = 0; } } } @@ -805,13 +805,14 @@ namespace Cantera { mat(cation[0],anion[0]) = mat(anion[0],cation[0]) = (1+vP/vM)*(-eps*xB*(1-eps*xA)*inv_vP_vM_MutualDiff)-zP*zM*Faraday*Faraday/GasConstant/temp/kappa/vol; mat(cation[1],anion[0]) = mat(anion[0],cation[1]) = (1+vP/vM)*(eps*xA*(1+eps*xB)*inv_vP_vM_MutualDiff)-zP*zM*Faraday*Faraday/GasConstant/temp/kappa/vol; - +/* for ( i = 0; i < nsp; i++ ) { for ( j = 0; j < nsp; j++ ) { mat(i,j) = 1.0/mat(i,j); //cout << "D" << i << j << " = " << mat(i,j) << endl; } } +*/ } @@ -871,8 +872,7 @@ mat(cation[1],anion[0]) = mat(anion[0],cation[1]) = (1+vP/vM)*(eps*xA*(1+eps*xB) mat.resize(nsp,nsp, 0.0); for (int i = 0; i < nsp; i++) for (int j = 0; j < nsp; j++) { - mat(i,j) = GasConstant * temp - / ( 6.0 * Pi * radiusSpec[i] * viscSpec[j] ) ; + mat(i,j) = ( 6.0 * Pi * radiusSpec[i] * viscSpec[j] ) / GasConstant / temp; } delete radiusSpec; delete viscSpec;