diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt index 47e573b6f13..28ff5208229 100644 --- a/src/coreComponents/constitutive/CMakeLists.txt +++ b/src/coreComponents/constitutive/CMakeLists.txt @@ -43,21 +43,23 @@ set( constitutive_headers fluid/multifluid/blackOil/PVTOData.hpp fluid/multifluid/CO2Brine/CO2BrineFluid.hpp fluid/multifluid/CO2Brine/PhaseModel.hpp - fluid/multifluid/CO2Brine/functions/PhillipsBrineDensity.hpp - fluid/multifluid/CO2Brine/functions/PhillipsBrineViscosity.hpp + fluid/multifluid/CO2Brine/functions/BrineEnthalpy.hpp + fluid/multifluid/CO2Brine/functions/CO2Enthalpy.hpp + fluid/multifluid/CO2Brine/functions/CO2EOSSolver.hpp + fluid/multifluid/CO2Brine/functions/CO2Solubility.hpp + fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.hpp + fluid/multifluid/CO2Brine/functions/CO2SolubilitySpycherPruess.hpp fluid/multifluid/CO2Brine/functions/EzrokhiBrineDensity.hpp fluid/multifluid/CO2Brine/functions/EzrokhiBrineViscosity.hpp - fluid/multifluid/CO2Brine/functions/CO2Solubility.hpp fluid/multifluid/CO2Brine/functions/FenghourCO2Viscosity.hpp fluid/multifluid/CO2Brine/functions/FlashModelBase.hpp - fluid/multifluid/CO2Brine/functions/PVTFunctionBase.hpp fluid/multifluid/CO2Brine/functions/NoOpPVTFunction.hpp + fluid/multifluid/CO2Brine/functions/PhillipsBrineDensity.hpp + fluid/multifluid/CO2Brine/functions/PhillipsBrineViscosity.hpp + fluid/multifluid/CO2Brine/functions/PureWaterProperties.hpp + fluid/multifluid/CO2Brine/functions/PVTFunctionBase.hpp fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp fluid/multifluid/CO2Brine/functions/SpanWagnerCO2Density.hpp - fluid/multifluid/CO2Brine/functions/BrineEnthalpy.hpp - fluid/multifluid/CO2Brine/functions/CO2Enthalpy.hpp - fluid/multifluid/CO2Brine/functions/CO2EOSSolver.hpp - fluid/multifluid/CO2Brine/functions/PureWaterProperties.hpp fluid/multifluid/CO2Brine/functions/WaterDensity.hpp fluid/multifluid/compositional/functions/CompositionalProperties.hpp fluid/multifluid/compositional/functions/CubicEOSPhaseModel.hpp @@ -202,6 +204,8 @@ set( constitutive_sources fluid/multifluid/CO2Brine/functions/EzrokhiBrineDensity.cpp fluid/multifluid/CO2Brine/functions/EzrokhiBrineViscosity.cpp fluid/multifluid/CO2Brine/functions/CO2Solubility.cpp + fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.cpp + fluid/multifluid/CO2Brine/functions/CO2SolubilitySpycherPruess.cpp fluid/multifluid/CO2Brine/functions/FenghourCO2Viscosity.cpp fluid/multifluid/CO2Brine/functions/SpanWagnerCO2Density.cpp fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.cpp diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp index 0a02e19ad51..b30d367dd4c 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.cpp @@ -377,9 +377,10 @@ void CO2BrineFluid< PHASE1, PHASE2, FLASH >::createPVTModels() // If 1 table is provided, it is the CO2 solubility table and water vapourisation is zero // If 2 tables are provided, they are the CO2 solubility and water vapourisation tables depending // on how phaseNames is arranged + string const solubilityModel = EnumStrings< CO2Solubility::SolubilityModel >::toString( CO2Solubility::SolubilityModel::Tables ); string_array strs; strs.emplace_back( "FlashModel" ); - strs.emplace_back( "Tables" ); // Marker to indicate that tables are provided + strs.emplace_back( solubilityModel ); // Marker to indicate that tables are provided strs.emplace_back( "" ); // 2 empty strings for the 2 phase tables gas first, then water strs.emplace_back( "" ); if( m_solubilityTables.size() == 2 ) diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.hpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.hpp index b68851478f2..835a475fec1 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/CO2BrineFluid.hpp @@ -194,7 +194,6 @@ class CO2BrineFluid : public MultiFluidBase // Flash model std::unique_ptr< FLASH > m_flash; - }; // these aliases are useful in constitutive dispatch diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.cpp index ce4041bce68..f43009e5204 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.cpp @@ -17,9 +17,8 @@ */ #include "constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.hpp" - -#include "constitutive/fluid/multifluid/CO2Brine/functions/CO2EOSSolver.hpp" -#include "constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp" +#include "constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilitySpycherPruess.hpp" +#include "constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.hpp" #include "functions/FunctionManager.hpp" #include "common/Units.hpp" @@ -29,337 +28,185 @@ namespace geos using namespace stringutilities; -namespace constitutive -{ - -namespace PVTProps -{ - namespace { -constexpr real64 P_Pa_f = 1e+5; -constexpr real64 P_c = 73.773 * P_Pa_f; -constexpr real64 T_c = 304.1282; -constexpr real64 Rgas = constants::gasConstant; -constexpr real64 V_c = Rgas*T_c/P_c; - -// these coefficients are in Table (A1) of Duan and Sun (2003) -constexpr real64 acoef[] = -{ 8.99288497e-2, -4.94783127e-1, 4.77922245e-2, 1.03808883e-2, -2.82516861e-2, 9.49887563e-2, 5.20600880e-4, - -2.93540971e-4, -1.77265112e-3, -2.51101973e-5, 8.93353441e-5, 7.88998563e-5, -1.66727022e-2, 1.398, 2.96e-2 }; - -real64 co2EOS( real64 const & T, real64 const & P, real64 const & V_r ) -{ - // reduced pressure - real64 const P_r = P*P_Pa_f/P_c; - // reduced temperature - real64 const T_r = units::convertCToK( T )/T_c; - - // CO2 equation of state - // see equation (A1) in Duan and Sun (2003) - real64 const f_Z = 1.0 - + ( acoef[0] + acoef[1]/(T_r * T_r) + acoef[2]/(T_r * T_r * T_r) )/V_r - + ( acoef[3] + acoef[4]/(T_r * T_r) + acoef[5]/(T_r * T_r * T_r) )/(V_r*V_r) - + ( acoef[6] + acoef[7]/(T_r * T_r) + acoef[8]/(T_r * T_r * T_r) )/(V_r*V_r*V_r*V_r) - + ( acoef[9] + acoef[10]/(T_r * T_r) + acoef[11]/(T_r * T_r * T_r) )/(V_r*V_r*V_r*V_r*V_r) - + acoef[12]/(T_r * T_r * T_r)/(V_r * V_r) * (acoef[13] + acoef[14]/(V_r * V_r)) * exp( -acoef[14]/(V_r * V_r)) - P_r * V_r / T_r; - - return f_Z; -} - -real64 PWater( real64 const & T ) +TableFunction const * makeTable( string const & tableName, + constitutive::PVTProps::PTTableCoordinates const & tableCoords, + array1d< real64 > && values, + FunctionManager & functionManager ) { - // these coefficients are defined in Table (B1) of Duan and Sun (2003) - constexpr real64 ccoef[] = { -38.640844, 5.8948420, 59.876516, 26.654627, 10.637097 }; - - // H2O critical pressure (bars) - real64 const P_c_w = 220.85; - // H2O critical temperature (K) - real64 const T_c_w = 647.29; - real64 const tt = ( units::convertCToK( T )-T_c_w )/T_c_w; - // Empirical model for water pressure of equation (B1) of Duan and Sun (2003) - real64 const x = ( P_c_w*units::convertCToK( T )/T_c_w ) - * (1 - + ccoef[0]*pow( -tt, 1.9 ) - + ccoef[1]*tt - + ccoef[2]*tt*tt - + ccoef[3]*tt*tt*tt - + ccoef[4]*tt*tt*tt*tt); - - return x; -} - -real64 logF( real64 const & T, real64 const & P, real64 const & V_r ) -{ - // reduced pressure - real64 const P_r = P*P_Pa_f/P_c; - // reduced temperature - real64 const T_r = units::convertCToK( T ) / T_c; - real64 const Z = P_r * V_r/T_r; - - // fugacity coefficient of CO2, equation (A6) of Duan and Sun (2003) - real64 const log_f = Z - 1 - log( Z ) + - ( acoef[0] + acoef[1]/T_r/T_r + acoef[2]/T_r/T_r/T_r )/V_r - + ( acoef[3] + acoef[4]/T_r/T_r + acoef[5]/T_r/T_r/T_r )/2.0/V_r/V_r - + ( acoef[6] + acoef[7]/T_r/T_r + acoef[8]/T_r/T_r/T_r )/4.0/V_r/V_r/V_r/V_r - + ( acoef[9] + acoef[10]/T_r/T_r + acoef[11]/T_r/T_r/T_r )/5.0/V_r/V_r/V_r/V_r/V_r - + acoef[12]/2.0/T_r/T_r/T_r/acoef[14] * ( acoef[13] + 1.0 - (acoef[13] + 1.0 + acoef[14]/V_r/V_r) * exp( -acoef[14]/V_r/V_r ) ); - - return log_f; -} - -real64 Par( real64 const & T, real64 const & P, real64 const * cc ) -{ - // "equation for the parameters", see equation (7) of Duan and Sun (2003) - real64 x = cc[0] - + cc[1]*T - + cc[2]/T - + cc[3]*T*T - + cc[4]/(630.0-T) - + cc[5]*P - + cc[6]*P *log( T ) - + cc[7]*P/T - + cc[8]*P/(630.0-T) - + cc[9]*P*P/(630.0-T)/(630.0-T) - + cc[10]*T *log( P ); - - return x; -} - -real64 CO2SolubilityFunction( string const & name, - real64 const & tolerance, - real64 const & T, - real64 const & P, - real64 (* f)( real64 const & x1, real64 const & x2, real64 const & x3 ) ) -{ - // compute the initial guess for Newton's method - real64 const initialReducedVolume = 0.75*Rgas*units::convertCToK( T )/(P*P_Pa_f)*(1/V_c); - - // define the local solver parameters - // for now, this is hard-coded, but we may want to let the user access the parameters at some point - integer const maxNumNewtonIter = 500; - integer const maxNumBacktrackIter = 8; - real64 const maxAbsUpdate = 1e12; - real64 const minAbsDeriv = 0; - real64 const allowedMinValue = 0.05; // value chosen to match previous implementation - real64 const presMultiplierForReporting = 1e5; // this is because P is in hectopascal in this function - - // solve the CO2 equation of state for this pair of (pres, temp) - // return the reduced volume - return CO2EOSSolver::solve( name, - maxNumNewtonIter, - maxNumBacktrackIter, - tolerance, - minAbsDeriv, - maxAbsUpdate, - allowedMinValue, - initialReducedVolume, - T, - P, - presMultiplierForReporting, - f ); + TableFunction * tableFunction = nullptr; + if( functionManager.hasGroup< TableFunction >( tableName ) ) + { + tableFunction = functionManager.getGroupPointer< TableFunction >( tableName ); + } + else + { + tableFunction = dynamicCast< TableFunction * >( functionManager.createChild( "TableFunction", tableName ) ); + tableFunction->setTableCoordinates( tableCoords.getCoords(), { units::Pressure, units::TemperatureInC } ); + tableFunction->setTableValues( values, units::Solubility ); + tableFunction->setInterpolationMethod( TableFunction::InterpolationType::Linear ); + } + tableFunction->initializeFunction(); + return tableFunction; } -void calculateCO2Solubility( string const & functionName, - real64 const & tolerance, - PTTableCoordinates const & tableCoords, - real64 const & salinity, - array1d< real64 > const & values ) +std::pair< TableFunction const *, TableFunction const * > +makeSolubilityTables( string const & functionName, + string_array const & inputParams, + constitutive::PVTProps::CO2Solubility::SolubilityModel const & solubilityModel ) { - // Interaction parameters, see Table 2 of Duan and Sun (2003) - constexpr real64 mu[] = - { 28.9447706, -0.0354581768, -4770.67077, 1.02782768e-5, 33.8126098, 9.04037140e-3, - -1.14934031e-3, -0.307405726, -0.0907301486, 9.32713393e-4, 0 }; - constexpr real64 lambda[] = { -0.411370585, 6.07632013e-4, 97.5347708, 0, 0, 0, 0, -0.0237622469, 0.0170656236, 0, 1.41335834e-5 }; - constexpr real64 zeta[] = { 3.36389723e-4, -1.98298980e-5, 0, 0, 0, 0, 0, 2.12220830e-3, -5.24873303e-3, 0, 0 }; + FunctionManager & functionManager = FunctionManager::getInstance(); + constitutive::PVTProps::PTTableCoordinates tableCoords; - localIndex const nPressures = tableCoords.nPressures(); - localIndex const nTemperatures = tableCoords.nTemperatures(); - - for( localIndex i = 0; i < nPressures; ++i ) + // Check solubility model for explicit table input + if( solubilityModel == constitutive::PVTProps::CO2Solubility::SolubilityModel::Tables ) { - real64 const P = tableCoords.getPressure( i ) / P_Pa_f; - - for( localIndex j = 0; j < nTemperatures; ++j ) + // The default table is a table with all zeros unless the name is explicitly provided + // The pressure and temperature values below will be used only to create the zero table so they + // simply give a range large enough to cover most values. + tableCoords.appendPressure( 1.0e5 ).appendPressure( 1.0e8 ) + .appendTemperature( 0.0 ).appendTemperature( 800.0 ); + + TableFunction const * tables[2] = { nullptr, nullptr }; + for( integer tableIndex : { 0, 1 } ) { - real64 const T = tableCoords.getTemperature( j ); + array1d< real64 > values( 4 ); + values.zero(); - // compute reduced volume by solving the CO2 equation of state - real64 const V_r = CO2SolubilityFunction( functionName, tolerance, T, P, &co2EOS ); - - // compute equation (6) of Duan and Sun (2003) - real64 const logK = Par( units::convertCToK( T ), P, mu ) - - logF( T, P, V_r ) - + 2*Par( units::convertCToK( T ), P, lambda ) * salinity - + Par( units::convertCToK( T ), P, zeta ) * salinity * salinity; - real64 const expLogK = exp( logK ); - - // mole fraction of CO2 in vapor phase, equation (4) of Duan and Sun (2003) - real64 const Pw = PWater( T ); - real64 const y_CO2 = (P - Pw)/P; - values[j*nPressures+i] = y_CO2 * P / expLogK; - - GEOS_WARNING_IF( expLogK <= 1e-10, - GEOS_FMT( "CO2Solubility: exp(logK) = {} is too small (logK = {}, P = {}, T = {}, V_r = {}), resulting solubility value is {}", - expLogK, logK, P, T, V_r, values[j*nPressures+i] )); - - if( values[j*nPressures+i] < 0 ) + string inputTableName = inputParams[2 + tableIndex]; + if( inputTableName.empty() ) + { + inputTableName = GEOS_FMT( "{}_zeroDissolution_table", constitutive::PVTProps::CO2Solubility::catalogName() ); + } + else { - GEOS_LOG_RANK_0( GEOS_FMT( "CO2Solubility: negative solubility value = {}, y_CO2 = {}, P = {}, PWater(T) = {}; corrected to 0", - values[j * nPressures + i], y_CO2, P, Pw ) ); - values[j*nPressures+i] = 0.0; + // If a name is explicitly given, then check that it exists + GEOS_THROW_IF( !functionManager.hasGroup< TableFunction >( inputTableName ), + GEOS_FMT( "{}: Could not find TableFunction with name {}", functionName, inputTableName ), + InputError ); } + tables[tableIndex] = makeTable( inputTableName, tableCoords, std::move ( values ), functionManager ); } + return { tables[0], tables[1] }; } -} -TableFunction const * getSolubilityTable( string const & tableName, - FunctionManager & functionManager ) -{ - TableFunction * const table = functionManager.getGroupPointer< TableFunction >( tableName ); - table->initializeFunction(); - table->setDimUnits( { units::Pressure, units::TemperatureInC } ); - table->setValueUnits( units::Solubility ); - return table; -} + // Initialize the (p,T) coordinates + constitutive::PVTProps::PVTFunctionHelpers::initializePropertyTable( inputParams, tableCoords ); -TableFunction const * makeSolubilityTable( array1d< real64_array > const & coords, - array1d< real64 > const & values, - string const & tableName, - FunctionManager & functionManager ) -{ - TableFunction * const table = dynamicCast< TableFunction * >( functionManager.createChild( "TableFunction", tableName ) ); - table->setTableCoordinates( coords, { units::Pressure, units::TemperatureInC } ); - table->setTableValues( values, units::Solubility ); - table->setInterpolationMethod( TableFunction::InterpolationType::Linear ); - return table; -} + // Initialize salinity and tolerance + GEOS_THROW_IF_LT_MSG( inputParams.size(), 9, + GEOS_FMT( "{}: insufficient number of model parameters", functionName ), + InputError ); -TableFunction const * makeZeroTable( string const & tableName, - FunctionManager & functionManager ) -{ - if( functionManager.hasGroup< TableFunction >( tableName ) ) - { - return getSolubilityTable( tableName, functionManager ); - } - else + real64 tolerance = 1e-9; + real64 salinity = 0.0; + try { - array1d< array1d< real64 > > coords( 2 ); - for( integer dim = 0; dim < 2; ++dim ) + salinity = stod( inputParams[8] ); + if( inputParams.size() >= 10 ) { - coords[dim].emplace_back( -1.0e10 ); - coords[dim].emplace_back( 1.0e10 ); + tolerance = stod( inputParams[9] ); } - array1d< real64 > values( 4 ); - values.zero(); - - return makeSolubilityTable( coords, values, tableName, functionManager ); } -} - -TableFunction const * makeSolubilityTable( string_array const & inputParams, - string const & functionName, - FunctionManager & functionManager ) -{ - // Check the second argument - if( inputParams[1] == "Tables" ) + catch( const std::invalid_argument & e ) { - string const inputTableName = inputParams[2]; - if( inputTableName.empty()) - { - return makeZeroTable( GEOS_FMT( "{}_zeroDissolution_table", CO2Solubility::catalogName() ), functionManager ); - } - else - { - GEOS_THROW_IF( !functionManager.hasGroup< TableFunction >( inputTableName ), - GEOS_FMT( "{}: Could not find TableFunction with name {}", functionName, inputTableName ), - InputError ); - return getSolubilityTable( inputTableName, functionManager ); - } + GEOS_THROW( GEOS_FMT( "{}: invalid model parameter value: {}", functionName, e.what() ), InputError ); } - string const tableName = functionName + "_co2Dissolution_table"; + integer const nPressures = tableCoords.nPressures(); + integer const nTemperatures = tableCoords.nTemperatures(); - if( functionManager.hasGroup< TableFunction >( tableName ) ) + array1d< real64 > co2Solubility( nPressures * nTemperatures ); + array1d< real64 > h2oSolubility( nPressures * nTemperatures ); + + if( solubilityModel == constitutive::PVTProps::CO2Solubility::SolubilityModel::DuanSun ) { - return getSolubilityTable( tableName, functionManager ); + constitutive::PVTProps::CO2SolubilityDuanSun::populateSolubilityTables( + functionName, + tableCoords, + salinity, + tolerance, + co2Solubility, + h2oSolubility ); } - else + else if( solubilityModel == constitutive::PVTProps::CO2Solubility::SolubilityModel::SpycherPruess ) { - // initialize the (p,T) coordinates - PTTableCoordinates tableCoords; - PVTFunctionHelpers::initializePropertyTable( inputParams, tableCoords ); - - // initialize salinity and tolerance - GEOS_THROW_IF_LT_MSG( inputParams.size(), 9, - GEOS_FMT( "{}: insufficient number of model parameters", functionName ), - InputError ); - - real64 tolerance = 1e-9; - real64 salinity = 0.0; - try + constitutive::PVTProps::CO2SolubilitySpycherPruess::populateSolubilityTables( + functionName, + tableCoords, + salinity, + tolerance, + co2Solubility, + h2oSolubility ); + } + + // Truncate negative solubility and warn + integer constexpr maxBad = 5; // Maximum number of bad values to report + stackArray2d< real64, maxBad *4 > badValues( maxBad, 4 ); + integer badCount = 0; + for( localIndex i = 0; i < nPressures; ++i ) + { + real64 const P = tableCoords.getPressure( i ); + for( localIndex j = 0; j < nTemperatures; ++j ) { - salinity = stod( inputParams[8] ); - if( inputParams.size() >= 10 ) + real64 const T = tableCoords.getTemperature( j ); + if( co2Solubility[j*nPressures+i] < 0.0 || h2oSolubility[j*nPressures+i] < 0.0 ) { - tolerance = stod( inputParams[9] ); + badValues( badCount % maxBad, 0 ) = P; + badValues( badCount % maxBad, 1 ) = T; + badValues( badCount % maxBad, 2 ) = co2Solubility[j*nPressures+i]; + badValues( badCount % maxBad, 3 ) = h2oSolubility[j*nPressures+i]; + ++badCount; + } + if( co2Solubility[j*nPressures+i] < 0.0 ) + { + co2Solubility[j*nPressures+i] = 0.0; + } + if( h2oSolubility[j*nPressures+i] < 0.0 ) + { + h2oSolubility[j*nPressures+i] = 0.0; } } - catch( const std::invalid_argument & e ) - { - GEOS_THROW( GEOS_FMT( "{}: invalid model parameter value: {}", functionName, e.what() ), InputError ); - } - - array1d< real64 > values( tableCoords.nPressures() * tableCoords.nTemperatures() ); - calculateCO2Solubility( functionName, tolerance, tableCoords, salinity, values ); - - return makeSolubilityTable( tableCoords.getCoords(), values, tableName, functionManager ); } -} -TableFunction const * makeVapourisationTable( string_array const & inputParams, - string const & functionName, - FunctionManager & functionManager ) -{ - if( inputParams[1] == "Tables" ) + if( 0 < badCount ) { - string const inputTableName = inputParams[3]; - if( inputTableName.empty()) - { - return makeZeroTable( GEOS_FMT( "{}_zeroDissolution_table", CO2Solubility::catalogName() ), functionManager ); - } - else + std::ostringstream badValueTable; + badValueTable + << std::setw( 15 ) << "Pressure (Pa)" << " " + << std::setw( 15 ) << "Temperature (C)" << " " + << std::setw( 23 ) << "CO2 solubility (mol/kg)" << " " + << std::setw( 15 ) << "H2O solubility (mol/kg)" << " " + << "\n"; + for( integer row = 0; row < LvArray::math::min( maxBad, badCount ); ++row ) { - GEOS_THROW_IF( !functionManager.hasGroup< TableFunction >( inputTableName ), - GEOS_FMT( "{}: Could not find TableFunction with name {}", functionName, inputTableName ), - InputError ); - return getSolubilityTable( inputTableName, functionManager ); + badValueTable + << std::setw( 15 ) << badValues( row, 0 ) << " " + << std::setw( 15 ) << badValues( row, 1 ) << " " + << std::setw( 23 ) << badValues( row, 2 ) << " " + << std::setw( 15 ) << badValues( row, 3 ) << " " + << "\n"; } + GEOS_LOG_RANK_0( GEOS_FMT( "CO2Solubility: {} negative solubility values encountered. These will be truncated to zero.\n Check out report table {} with max {} values", + badCount, badValueTable.str(), maxBad ) ); } - string const tableName = functionName + "_waterVaporization_table"; - - if( functionManager.hasGroup< TableFunction >( tableName ) ) - { - return getSolubilityTable( tableName, functionManager ); - } - else - { - // initialize the (p,T) coordinates - PTTableCoordinates tableCoords; - PVTFunctionHelpers::initializePropertyTable( inputParams, tableCoords ); - - // Currently initialise to all zeros - array1d< real64 > values( tableCoords.nPressures() * tableCoords.nTemperatures() ); - values.zero(); - - return makeSolubilityTable( tableCoords.getCoords(), values, tableName, functionManager ); - } + string const co2TableName = functionName + "_co2Dissolution_table"; + TableFunction const * co2SolubilityTable = makeTable( co2TableName, tableCoords, std::move( co2Solubility ), functionManager ); + string const h2oTableName = functionName + "_waterVaporization_table"; + TableFunction const * h2oSolubilityTable = makeTable( h2oTableName, tableCoords, std::move( h2oSolubility ), functionManager ); + return {co2SolubilityTable, h2oSolubilityTable}; } } // namespace +namespace constitutive +{ +namespace PVTProps +{ + CO2Solubility::CO2Solubility( string const & name, string_array const & inputParams, string_array const & phaseNames, @@ -389,8 +236,18 @@ CO2Solubility::CO2Solubility( string const & name, string const expectedWaterPhaseNames[] = { "Water", "water", "Liquid", "liquid" }; m_phaseLiquidIndex = PVTFunctionHelpers::findName( phaseNames, expectedWaterPhaseNames, "phaseNames" ); - m_CO2SolubilityTable = makeSolubilityTable( inputParams, m_modelName, FunctionManager::getInstance() ); - m_WaterVapourisationTable = makeVapourisationTable( inputParams, m_modelName, FunctionManager::getInstance() ); + SolubilityModel solubilityModel = SolubilityModel::DuanSun; // Default solubility model + if( inputParams[1] == EnumStrings< SolubilityModel >::toString( SolubilityModel::Tables ) ) + { + solubilityModel = SolubilityModel::Tables; + } + else if( 11 <= inputParams.size() ) + { + solubilityModel = EnumStrings< SolubilityModel >::fromString( inputParams[10] ); + } + + std::tie( m_CO2SolubilityTable, m_WaterVapourisationTable ) = makeSolubilityTables( m_modelName, inputParams, solubilityModel ); + if( printTable ) { m_CO2SolubilityTable->print( m_CO2SolubilityTable->getName() ); @@ -418,8 +275,6 @@ CO2Solubility::KernelWrapper CO2Solubility::createKernelWrapper() const m_phaseLiquidIndex ); } -REGISTER_CATALOG_ENTRY( FlashModelBase, CO2Solubility, string const &, string_array const &, string_array const &, string_array const &, array1d< real64 > const &, bool const ) - } // end namespace PVTProps } // namespace constitutive diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.hpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.hpp index bc0492008fa..9221c48503b 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.hpp @@ -103,7 +103,14 @@ class CO2SolubilityUpdate final : public FlashModelBaseUpdate class CO2Solubility : public FlashModelBase { public: + enum class SolubilityModel : integer + { + DuanSun, + SpycherPruess, + Tables + }; +public: CO2Solubility( string const & name, string_array const & inputParams, string_array const & phaseNames, @@ -130,7 +137,6 @@ class CO2Solubility : public FlashModelBase KernelWrapper createKernelWrapper() const; private: - /// Table to compute solubility as a function of pressure and temperature TableFunction const * m_CO2SolubilityTable; @@ -332,6 +338,11 @@ CO2SolubilityUpdate::compute( real64 const & pressure, } } +ENUM_STRINGS( CO2Solubility::SolubilityModel, + "DuanSun", + "SpycherPruess", + "Tables" ); + } // end namespace PVTProps } // end namespace constitutive diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.cpp new file mode 100644 index 00000000000..cc7b5f3bab2 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.cpp @@ -0,0 +1,222 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CO2SolubilityDuanSun.cpp + */ + +#include "constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.hpp" +#include "constitutive/fluid/multifluid/CO2Brine/functions/CO2EOSSolver.hpp" + +#include "common/Units.hpp" + +namespace geos +{ +namespace constitutive +{ +namespace PVTProps +{ + +namespace +{ + +constexpr real64 P_Pa_f = 1e+5; +constexpr real64 P_c = 73.773 * P_Pa_f; +constexpr real64 T_c = 304.1282; +constexpr real64 Rgas = constants::gasConstant; +constexpr real64 V_c = Rgas*T_c/P_c; + +// these coefficients are in Table (A1) of Duan and Sun (2003) +constexpr real64 acoef[] = +{ 8.99288497e-2, -4.94783127e-1, 4.77922245e-2, 1.03808883e-2, -2.82516861e-2, 9.49887563e-2, 5.20600880e-4, + -2.93540971e-4, -1.77265112e-3, -2.51101973e-5, 8.93353441e-5, 7.88998563e-5, -1.66727022e-2, 1.398, 2.96e-2 }; + +real64 co2EOS( real64 const & T, real64 const & P, real64 const & V_r ) +{ + // reduced pressure + real64 const P_r = P*P_Pa_f/P_c; + // reduced temperature + real64 const T_r = units::convertCToK( T )/T_c; + + // CO2 equation of state + // see equation (A1) in Duan and Sun (2003) + real64 const f_Z = 1.0 + + ( acoef[0] + acoef[1]/(T_r * T_r) + acoef[2]/(T_r * T_r * T_r) )/V_r + + ( acoef[3] + acoef[4]/(T_r * T_r) + acoef[5]/(T_r * T_r * T_r) )/(V_r*V_r) + + ( acoef[6] + acoef[7]/(T_r * T_r) + acoef[8]/(T_r * T_r * T_r) )/(V_r*V_r*V_r*V_r) + + ( acoef[9] + acoef[10]/(T_r * T_r) + acoef[11]/(T_r * T_r * T_r) )/(V_r*V_r*V_r*V_r*V_r) + + acoef[12]/(T_r * T_r * T_r)/(V_r * V_r) * (acoef[13] + acoef[14]/(V_r * V_r)) * exp( -acoef[14]/(V_r * V_r)) - P_r * V_r / T_r; + + return f_Z; +} + +real64 PWater( real64 const & T ) +{ + // these coefficients are defined in Table (B1) of Duan and Sun (2003) + constexpr real64 ccoef[] = { -38.640844, 5.8948420, 59.876516, 26.654627, 10.637097 }; + + // H2O critical pressure (bars) + real64 const P_c_w = 220.85; + // H2O critical temperature (K) + real64 const T_c_w = 647.29; + real64 const tt = ( units::convertCToK( T )-T_c_w )/T_c_w; + // Empirical model for water pressure of equation (B1) of Duan and Sun (2003) + real64 const x = ( P_c_w*units::convertCToK( T )/T_c_w ) + * (1 + + ccoef[0]*pow( -tt, 1.9 ) + + ccoef[1]*tt + + ccoef[2]*tt*tt + + ccoef[3]*tt*tt*tt + + ccoef[4]*tt*tt*tt*tt); + + return x; +} + +real64 logF( real64 const & T, real64 const & P, real64 const & V_r ) +{ + // reduced pressure + real64 const P_r = P*P_Pa_f/P_c; + // reduced temperature + real64 const T_r = units::convertCToK( T ) / T_c; + real64 const Z = P_r * V_r/T_r; + + // fugacity coefficient of CO2, equation (A6) of Duan and Sun (2003) + real64 const log_f = Z - 1 - log( Z ) + + ( acoef[0] + acoef[1]/T_r/T_r + acoef[2]/T_r/T_r/T_r )/V_r + + ( acoef[3] + acoef[4]/T_r/T_r + acoef[5]/T_r/T_r/T_r )/2.0/V_r/V_r + + ( acoef[6] + acoef[7]/T_r/T_r + acoef[8]/T_r/T_r/T_r )/4.0/V_r/V_r/V_r/V_r + + ( acoef[9] + acoef[10]/T_r/T_r + acoef[11]/T_r/T_r/T_r )/5.0/V_r/V_r/V_r/V_r/V_r + + acoef[12]/2.0/T_r/T_r/T_r/acoef[14] * ( acoef[13] + 1.0 - (acoef[13] + 1.0 + acoef[14]/V_r/V_r) * exp( -acoef[14]/V_r/V_r ) ); + + return log_f; +} + +real64 Par( real64 const & T, real64 const & P, real64 const * cc ) +{ + // "equation for the parameters", see equation (7) of Duan and Sun (2003) + real64 x = cc[0] + + cc[1]*T + + cc[2]/T + + cc[3]*T*T + + cc[4]/(630.0-T) + + cc[5]*P + + cc[6]*P *log( T ) + + cc[7]*P/T + + cc[8]*P/(630.0-T) + + cc[9]*P*P/(630.0-T)/(630.0-T) + + cc[10]*T *log( P ); + + return x; +} + +real64 CO2SolubilityFunction( string const & name, + real64 const & tolerance, + real64 const & T, + real64 const & P, + real64 (* f)( real64 const & x1, real64 const & x2, real64 const & x3 ) ) +{ + // compute the initial guess for Newton's method + real64 const initialReducedVolume = 0.75*Rgas*units::convertCToK( T )/(P*P_Pa_f)*(1/V_c); + + // define the local solver parameters + // for now, this is hard-coded, but we may want to let the user access the parameters at some point + integer const maxNumNewtonIter = 500; + integer const maxNumBacktrackIter = 8; + real64 const maxAbsUpdate = 1e12; + real64 const minAbsDeriv = 0; + real64 const allowedMinValue = 0.05; // value chosen to match previous implementation + real64 const presMultiplierForReporting = 1e5; // this is because P is in hectopascal in this function + + // solve the CO2 equation of state for this pair of (pres, temp) + // return the reduced volume + return CO2EOSSolver::solve( name, + maxNumNewtonIter, + maxNumBacktrackIter, + tolerance, + minAbsDeriv, + maxAbsUpdate, + allowedMinValue, + initialReducedVolume, + T, + P, + presMultiplierForReporting, + f ); +} + +void calculateCO2Solubility( string const & functionName, + real64 const & tolerance, + PTTableCoordinates const & tableCoords, + real64 const & salinity, + array1d< real64 > const & values ) +{ + // Interaction parameters, see Table 2 of Duan and Sun (2003) + constexpr real64 mu[] = + { 28.9447706, -0.0354581768, -4770.67077, 1.02782768e-5, 33.8126098, 9.04037140e-3, + -1.14934031e-3, -0.307405726, -0.0907301486, 9.32713393e-4, 0 }; + constexpr real64 lambda[] = { -0.411370585, 6.07632013e-4, 97.5347708, 0, 0, 0, 0, -0.0237622469, 0.0170656236, 0, 1.41335834e-5 }; + constexpr real64 zeta[] = { 3.36389723e-4, -1.98298980e-5, 0, 0, 0, 0, 0, 2.12220830e-3, -5.24873303e-3, 0, 0 }; + + localIndex const nPressures = tableCoords.nPressures(); + localIndex const nTemperatures = tableCoords.nTemperatures(); + + for( localIndex i = 0; i < nPressures; ++i ) + { + real64 const P = tableCoords.getPressure( i ) / P_Pa_f; + + for( localIndex j = 0; j < nTemperatures; ++j ) + { + real64 const T = tableCoords.getTemperature( j ); + + // compute reduced volume by solving the CO2 equation of state + real64 const V_r = CO2SolubilityFunction( functionName, tolerance, T, P, &co2EOS ); + + // compute equation (6) of Duan and Sun (2003) + real64 const logK = Par( units::convertCToK( T ), P, mu ) + - logF( T, P, V_r ) + + 2*Par( units::convertCToK( T ), P, lambda ) * salinity + + Par( units::convertCToK( T ), P, zeta ) * salinity * salinity; + real64 const expLogK = exp( logK ); + + // mole fraction of CO2 in vapor phase, equation (4) of Duan and Sun (2003) + real64 const Pw = PWater( T ); + real64 const y_CO2 = (P - Pw)/P; + values[j*nPressures+i] = y_CO2 * P / expLogK; + + GEOS_WARNING_IF( expLogK <= 1e-10, + GEOS_FMT( "CO2Solubility: exp(logK) = {} is too small (logK = {}, P = {}, T = {}, V_r = {}), resulting solubility value is {}", + expLogK, logK, P, T, V_r, values[j*nPressures+i] )); + } + } +} + +} // end namespace + +void CO2SolubilityDuanSun::populateSolubilityTables( string const & functionName, + PTTableCoordinates const & tableCoords, + real64 const & salinity, + real64 const & tolerance, + array1d< real64 > const & co2SolubilityValues, + array1d< real64 > const & h2oSolubilityValues ) +{ + h2oSolubilityValues.zero(); + calculateCO2Solubility( functionName, + tolerance, + tableCoords, + salinity, + co2SolubilityValues ); +} + +} // end namespace PVTProps +} // namespace constitutive +} // end namespace geos diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.hpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.hpp new file mode 100644 index 00000000000..f0db58369b3 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilityDuanSun.hpp @@ -0,0 +1,57 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CO2SolubilityDuanSun.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_CO2BRINE_FUNCTIONS_CO2SOLUBILITYDUANSUN_HPP_ +#define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_CO2BRINE_FUNCTIONS_CO2SOLUBILITYDUANSUN_HPP_ + +#include "constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp" + +namespace geos +{ +namespace constitutive +{ +namespace PVTProps +{ + +struct CO2SolubilityDuanSun +{ +/** + * @brief Create CO2 and H2O solubility table based on Duan and Sun (2003) + * @details Each generated table is a 2D table with lookup properties pressure (in Pa) and + * temperature (in degC). The returned CO2 solubility is in mole of CO2 per kg of + * H2O and the returned water vapourisation is in moles of H2O per kg of CO2. + * @param[in] functionName The name of the model + * @param[in] tableCoords The values of pressure and temperature + * @param[in] salinity The salinity of the brine + * @param[in] tolerance Tolerance to be used in solving for the solubility + * @param[out] co2SolubilityValues The CO2 solubility values (mol/kg) at the given pressures and temperatures + * @param[out] h2oSolubilityValues The H2O solubility values (mol/kg) at the given pressures and temperatures + */ + static void populateSolubilityTables( string const & functionName, + PTTableCoordinates const & tableCoords, + real64 const & salinity, + real64 const & tolerance, + array1d< real64 > const & co2SolubilityValues, + array1d< real64 > const & h2oSolubilityValues ); +}; + +} // end namespace PVTProps +} // end namespace constitutive +} // end namespace geos + +#endif //GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_CO2BRINE_FUNCTIONS_CO2SOLUBILITYDUANSUN_HPP_ diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilitySpycherPruess.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilitySpycherPruess.cpp new file mode 100644 index 00000000000..70c4425a9b3 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilitySpycherPruess.cpp @@ -0,0 +1,202 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CO2SolubilitySpycherPruess.cpp + */ + +#include "CO2SolubilitySpycherPruess.hpp" +#include "constitutive/fluid/multifluid/CO2Brine/functions/SpanWagnerCO2Density.hpp" +#include "common/Units.hpp" + +namespace geos +{ + +namespace +{ +/** + * Physical constants + * These correlations are developed with units: pressure [bar], volume [cm3], temperature [K]. + * The equilibrium constant correlations (equilibriumConstantCO2 and equilibriumConstantH2O) have a temperature in C + */ +static constexpr real64 P_ref = 1.0; // Reference pressure is 1 bar +static constexpr real64 Pa_2_bar = 1.0e-5; // Conversion factor from Pa to Bar +static constexpr real64 R = 10.0*constants::gasConstant; // Universal gas constant in [bar.cm3/mol.K] +static constexpr real64 molarMassCO2 = 44.01e-3; // Molar mass of CO2 [kg/mol] +static constexpr real64 molarMassH2O = 18.01e-3; // Molar mass of H2O [kg/mol] +static constexpr real64 v_av_H2O = 18.1; // Average partial molar volume of H2O [cm3/mol] +static constexpr real64 v_av_CO2 = 32.6; // Average partial molar volume of CO2 [cm3/mol] + +/** + * @brief Calculate the CO2 equilibrium constant from Spycher et al. (2003) + * @details The correlation parameters are given in Table 2 of Spycher et al. (2003) + * @param[in] T Temperature [C] + */ +real64 equilibriumConstantCO2( real64 const T ) +{ + constexpr real64 c[] = {1.189, 1.304e-2, -5.446e-5}; + real64 const logk0_CO2 = c[0] + T * (c[1] + T * c[2]); + return pow( 10.0, logk0_CO2 ); +} + +/** + * @brief Calculate the H2O equilibrium constant from Spycher et al. (2003) + * @details The correlation parameters are given in Table 2 of Spycher et al. (2003) + * @param[in] T Temperature [C] + */ +real64 equilibriumConstantH2O( real64 const T ) +{ + constexpr real64 c[] = {-2.209, 3.097e-2, -1.098e-4, 2.048e-7}; + real64 const logk0_H2O = c[0] + T * (c[1] + T * (c[2] + T* c[3])); + return pow( 10.0, logk0_H2O ); +} + +/** + * @brief Calculate the fugacity coefficient of CO2 from the RK-EOS + * @param[in] P Pressure [bar] + * @param[in] T Temperature [K] + * @param[in] rhoCO2 Density of CO2 [kg/m3] + * @param[in] salinity Salinity of water + */ +real64 fugacityCoefficientCO2( real64 const P, real64 const T, real64 const rhoCO2, real64 const salinity ) +{ + GEOS_UNUSED_VAR( salinity ); + real64 const V = 1.0e6*molarMassCO2/rhoCO2; // Molar volume [cm3/mol] + + // Mixture parameters from the tuned Redlich-Kwong EOS + // These values are given in Table 1 of Spycher et al. (2003) + real64 const a_CO2 = (7.54e7 - 4.13e4*T); // [bar.cm3.K^(1/2)/mol^2] + real64 constexpr b_CO2 = 27.8; // [cm3/mol] + + real64 const lnPhiCO2 = log( V/(V - b_CO2)) + b_CO2/(V - b_CO2) + - 2*a_CO2/(R*pow( T, 1.5 )*b_CO2)*log((V + b_CO2)/V ) + + a_CO2*b_CO2/(R*pow( T, 1.5 )*b_CO2*b_CO2) * (log((V + b_CO2)/V ) - b_CO2/(V + b_CO2)) + - log( P*V/(R*T)); + return exp( lnPhiCO2 ); +} + +/** + * @brief Calculate the fugacity coefficient of CO2 from the RK-EOS + * @param[in] P Pressure [bar] + * @param[in] T Temperature [K] + * @param[in] rhoCO2 Density of CO2 [kg/m3] + * @param[in] salinity Salinity of water + */ +real64 fugacityCoefficientH2O( real64 const P, real64 const T, real64 const rhoCO2, real64 const salinity ) +{ + GEOS_UNUSED_VAR( salinity ); + real64 const V = 1.0e6*molarMassCO2/rhoCO2; // Molar volume [cm3/mol] + + // Mixture parameters from the tuned Redlich-Kwong EOS + // These values are given in Table 1 of Spycher et al. (2003) + real64 const a_CO2 = (7.54e7 - 4.13e4*T); // [bar.cm3.K^(1/2)/mol^2] + real64 constexpr a_CO2_H2O = 7.89e7; // [bar.cm3.K^(1/2)/mol^2] + real64 constexpr b_CO2 = 27.8; // [cm3/mol] + real64 constexpr b_H2O = 18.18; // [cm3/mol] + + real64 const lnPhiH2O = log( V/(V - b_CO2) ) + b_H2O/(V - b_CO2) + - 2.0 * a_CO2_H2O * log( (V + b_CO2)/V ) / ( R*pow( T, 1.5 )*b_CO2 ) + + a_CO2 * b_H2O / ( R*pow( T, 1.5 )*b_CO2*b_CO2 )*( log( (V + b_CO2)/V ) - b_CO2/(V + b_CO2) ) + - log( P*V/(R*T) ); + return exp( lnPhiH2O ); +} + +/** + * @brief Calculate the parameter A from Eq (11) of Spycher et al. (2003) + * @param[in] P Pressure [Pa] + * @param[in] T Temperature [C] + * @param[in] rhoCO2 CO2 density [kg/m3] + * @param[in] salinity Salinity of the water + */ +real64 computeA( real64 const P, real64 const T, real64 const rhoCO2, real64 const salinity ) +{ + real64 const P_in_bar = P * Pa_2_bar; + real64 const deltaP = P_in_bar - P_ref; + real64 const TinK = units::convertCToK( T ); + real64 const k0_H2O = equilibriumConstantH2O( T ); // K-value for H2O at 1 bar + real64 const phi_H2O = fugacityCoefficientH2O( P_in_bar, TinK, rhoCO2, salinity ); // Fugacity coefficient of H2O for the water-CO2 system + real64 const A = k0_H2O/(phi_H2O*P_in_bar) * exp( deltaP*v_av_H2O/(R*TinK)); + return A; +} + +/** + * @brief Calculate the parameter B from Eq (12) of Spycher et al. (2003) + * @param[in] P Pressure [Pa] + * @param[in] T Temperature [C] + * @param[in] rhoCO2 CO2 density [kg/m3] + * @param[in] salinity Salinity of the water + */ +real64 computeB( real64 const P, real64 const T, real64 const rhoCO2, real64 const salinity ) +{ + real64 const P_in_bar = P * Pa_2_bar; + real64 const deltaP = P_in_bar - P_ref; + real64 const TinK = units::convertCToK( T ); + real64 const k0_CO2 = equilibriumConstantCO2( T ); // K-value for CO2 at 1 bar + real64 const phi_CO2 = fugacityCoefficientCO2( P_in_bar, TinK, rhoCO2, salinity ); // Fugacity coefficient of CO2 for the water-CO2 system + real64 const B = phi_CO2*P_in_bar/(55.508*k0_CO2) * exp( -(deltaP*v_av_CO2)/(R*TinK) ); + return B; +} +} // end namespace + +namespace constitutive +{ +namespace PVTProps +{ +void CO2SolubilitySpycherPruess::populateSolubilityTables( string const & functionName, + PTTableCoordinates const & tableCoords, + real64 const & salinity, + real64 const & tolerance, + array1d< real64 > const & co2SolubilityValues, + array1d< real64 > const & h2oSolubilityValues ) +{ + localIndex const nPressures = tableCoords.nPressures(); + localIndex const nTemperatures = tableCoords.nTemperatures(); + + // Calculate the CO2 density + array1d< real64 > densities( nPressures*nTemperatures ); + SpanWagnerCO2Density::calculateCO2Density( GEOS_FMT( "{}_co2_density", functionName ), + tolerance, + tableCoords, + densities ); + + for( localIndex i = 0; i < nPressures; ++i ) + { + real64 const P = tableCoords.getPressure( i ); + + for( localIndex j = 0; j < nTemperatures; ++j ) + { + real64 const T = tableCoords.getTemperature( j ); + + // Get the CO2 density + real64 const rhoCO2 = densities[j*nPressures+i]; + + // Calculate A and B + real64 const A = computeA( P, T, rhoCO2, salinity ); + real64 const B = computeB( P, T, rhoCO2, salinity ); + + // Calculate the mole fractions + // Eqns (13) and (14) from Spycher et al. (2003) + real64 const y_H2O = (1.0 - B)/(1.0/A - B); + real64 const x_CO2 = B*(1.0 - y_H2O); + + // Calculate the solubility + co2SolubilityValues[j*nPressures+i] = x_CO2/((1.0 - x_CO2)*molarMassH2O); + h2oSolubilityValues[j*nPressures+i] = y_H2O/((1.0 - y_H2O)*molarMassCO2); + } + } +} + +} // end namespace PVTProps +} // namespace constitutive +} // end namespace geos diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilitySpycherPruess.hpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilitySpycherPruess.hpp new file mode 100644 index 00000000000..01c71019b86 --- /dev/null +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/CO2SolubilitySpycherPruess.hpp @@ -0,0 +1,59 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CO2SolubilitySpycherPruess.hpp + */ + +#ifndef GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_CO2BRINE_FUNCTIONS_CO2SOLUBILITYSPYCHERPRUESS_HPP_ +#define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_CO2BRINE_FUNCTIONS_CO2SOLUBILITYSPYCHERPRUESS_HPP_ + +#include "constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp" + +namespace geos +{ +namespace constitutive +{ +namespace PVTProps +{ + +struct CO2SolubilitySpycherPruess +{ + +/** + * @brief Create CO2 and H2O solubility table based on Spycher, Pruess, Ennis-King (2003) + * @details Each generated table is a 2D table with lookup properties pressure (in Pa) and + * temperature (in degC). The returned CO2 solubility is in mole of CO2 per kg of + * H2O and the returned water vapourisation is in moles of H2O per kg of CO2. + * @param[in] functionName The name of the model + * @param[in] tableCoords The values of pressure and temperature + * @param[in] salinity The salinity of the brine + * @param[in] tolerance Tolerance to be used in solving for the solubility + * @param[out] co2SolubilityValues The CO2 solubility values (mol/kg) at the given pressures and temperatures + * @param[out] h2oSolubilityValues The H2O solubility values (mol/kg) at the given pressures and temperatures + */ + static void populateSolubilityTables( string const & functionName, + PTTableCoordinates const & tableCoords, + real64 const & salinity, + real64 const & tolerance, + array1d< real64 > const & co2SolubilityValues, + array1d< real64 > const & h2oSolubilityValues ); + +}; + +} // end namespace PVTProps +} // end namespace constitutive +} // end namespace geos + +#endif //GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_CO2BRINE_FUNCTIONS_CO2SOLUBILITYSPYCHERPRUESS_HPP_ diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.cpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.cpp index 951fd1c22f6..f71c6ed46dd 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.cpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.cpp @@ -16,7 +16,6 @@ * @file PVTFunctionHelpers.cpp */ -#include "codingUtilities/StringUtilities.hpp" #include "constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp" #include "LvArray/src/sortedArrayManipulation.hpp" diff --git a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp index fd0cc5d2386..5f0e3f39276 100644 --- a/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp +++ b/src/coreComponents/constitutive/fluid/multifluid/CO2Brine/functions/PVTFunctionHelpers.hpp @@ -19,6 +19,8 @@ #include "common/DataTypes.hpp" #include "common/Units.hpp" +#include "codingUtilities/StringUtilities.hpp" + #ifndef GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_CO2BRINE_FUNCTIONS_PVTFUNCTIONHELPERS_HPP_ #define GEOS_CONSTITUTIVE_FLUID_MULTIFLUID_CO2BRINE_FUNCTIONS_PVTFUNCTIONHELPERS_HPP_ @@ -122,8 +124,8 @@ class PTTableCoordinates localIndex nPressures() const { return coords[coordType::PRES].size(); } localIndex nTemperatures() const { return coords[coordType::TEMP].size(); } - void appendPressure( const real64 & pres ) { coords[coordType::PRES].emplace_back( pres ); } - void appendTemperature( const real64 & temp ) { coords[coordType::TEMP].emplace_back( temp ); } + PTTableCoordinates & appendPressure( const real64 & pres ) { coords[coordType::PRES].emplace_back( pres ); return *this; } + PTTableCoordinates & appendTemperature( const real64 & temp ) { coords[coordType::TEMP].emplace_back( temp ); return *this; } real64 const & getPressure( localIndex i ) const { return coords[coordType::PRES][i]; } real64 const & getTemperature( localIndex i ) const { return coords[coordType::TEMP][i]; } diff --git a/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt b/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt index e8b27218b8f..5f03cf2eaaa 100644 --- a/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt +++ b/src/coreComponents/unitTests/constitutiveTests/CMakeLists.txt @@ -1,9 +1,11 @@ # Specify list of tests set( gtest_geosx_tests + testCapillaryPressure.cpp + testCO2BrinePVTModels.cpp + testCO2SpycherPruessModels.cpp testDamage.cpp testRelPerm.cpp - testRelPermHysteresis.cpp - testCapillaryPressure.cpp ) + testRelPermHysteresis.cpp ) set( gtest_triaxial_xmls testTriaxial_druckerPragerExtended.xml @@ -37,8 +39,7 @@ endif() if( ENABLE_PVTPackage ) list( APPEND gtest_geosx_tests - testMultiFluid.cpp - testCO2BrinePVTModels.cpp ) + testMultiFluid.cpp ) list( APPEND dependencyList PVTPackage ) endif() diff --git a/src/coreComponents/unitTests/constitutiveTests/testCO2SpycherPruessModels.cpp b/src/coreComponents/unitTests/constitutiveTests/testCO2SpycherPruessModels.cpp new file mode 100644 index 00000000000..7a248a4007b --- /dev/null +++ b/src/coreComponents/unitTests/constitutiveTests/testCO2SpycherPruessModels.cpp @@ -0,0 +1,374 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +// Source includes +#include "codingUtilities/UnitTestUtilities.hpp" +#include "constitutive/fluid/multifluid/CO2Brine/functions/CO2Solubility.hpp" +#include "constitutive/fluid/multifluid/MultiFluidUtils.hpp" +#include "mainInterface/GeosxState.hpp" +#include "mainInterface/initialization.hpp" +#include "functions/FunctionManager.hpp" + +// TPL includes +#include + +using namespace geos; +using namespace geos::testing; +using namespace geos::stringutilities; +using namespace geos::constitutive; +using namespace geos::constitutive::multifluid; +using namespace geos::constitutive::PVTProps; + +// Test parameter +using TestParam = std::tuple< + real64 const, // Pressure [Pa] + real64 const, // Temperature [C] + real64 const, // Total co2 mole fraction (z_co2 = 1-z_wat) + real64 const, // Expected vapour fraction (V) + real64 const, // Expected dissolved CO2 fraction (x_co2) + real64 const // Expected vapourised water fraction (y_wat) + >; + +class CO2SolubilitySpycherPruessTestFixture : public ::testing::TestWithParam< TestParam > +{ +protected: + static integer constexpr numPhase = 2; + static integer constexpr numComp = 2; + static integer constexpr numDof = numComp + 2; + static real64 constexpr relTol = 1.0e-5; + static real64 constexpr absTol = 1.0e-7; + static real64 constexpr pertubation = 1.0e-6; + static constexpr char const * flashContent = "FlashModel CO2Solubility 1.0e5 1.0e7 9.9e5 283.15 383.15 10.0 0.15 1.0e-8 SpycherPruess"; + +public: + CO2SolubilitySpycherPruessTestFixture() = default; + ~CO2SolubilitySpycherPruessTestFixture() override = default; + +protected: + static std::unique_ptr< CO2Solubility > makeFlashModel( string const & fileContent ); +}; + +std::unique_ptr< CO2Solubility > +CO2SolubilitySpycherPruessTestFixture::makeFlashModel( string const & fileContent ) +{ + // Define phase names + string_array phaseNames; + phaseNames.resize( 2 ); + phaseNames[0] = "gas"; + phaseNames[1] = "liquid"; + + // Define component names and molar weight + string_array componentNames; + componentNames.resize( 2 ); + componentNames[0] = "co2"; + componentNames[1] = "water"; + + array1d< real64 > componentMolarWeight; + componentMolarWeight.resize( 2 ); + componentMolarWeight[0] = 44.0e-3; + componentMolarWeight[1] = 18.0e-3; + + // Read file parameters + array1d< string > const strs = stringutilities::tokenizeBySpaces< array1d >( fileContent ); + + return std::make_unique< CO2Solubility >( strs[1], + strs, + phaseNames, + componentNames, + componentMolarWeight, + false ); +} + +TEST_P( CO2SolubilitySpycherPruessTestFixture, testExpectedValues ) +{ + auto flashModel = makeFlashModel( CO2SolubilitySpycherPruessTestFixture::flashContent ); + + auto [pressure, temperature, z_co2, expected_V, expected_x_co2, expected_y_wat] = GetParam(); + + StackArray< real64, 1, numComp > composition( 2 ); + composition[0] = z_co2; + composition[1] = 1.0 - z_co2; + + StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseFrac( 1, 1, numPhase ); + StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseFrac( 1, 1, numPhase, numDof ); + MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > + phaseFracAndDeriv { phaseFrac[0][0], dPhaseFrac[0][0] }; + + StackArray< real64, 4, numComp *numPhase, LAYOUT_PHASE_COMP > phaseCompFrac( 1, 1, numPhase, numComp ); + StackArray< real64, 5, numDof *numComp *numPhase, LAYOUT_PHASE_COMP_DC > dPhaseCompFrac( 1, 1, numPhase, numComp, numDof ); + MultiFluidVarSlice< real64, 2, USD_PHASE_COMP - 2, USD_PHASE_COMP_DC - 2 > + phaseCompFracAndDeriv { phaseCompFrac[0][0], dPhaseCompFrac[0][0] }; + + auto flashModelWrapper = flashModel->createKernelWrapper(); + + flashModelWrapper.compute( pressure, + temperature, + composition.toSliceConst(), + phaseFracAndDeriv, + phaseCompFracAndDeriv ); + + real64 const V = phaseFracAndDeriv.value[0]; + real64 const x_co2 = phaseCompFracAndDeriv.value[1][0]; + real64 const y_wat = phaseCompFracAndDeriv.value[0][1]; + + checkRelativeError( V, expected_V, relTol, absTol ); + checkRelativeError( x_co2, expected_x_co2, relTol, absTol ); + checkRelativeError( y_wat, expected_y_wat, relTol, absTol ); +} + +TEST_P( CO2SolubilitySpycherPruessTestFixture, testNumericalDerivatives ) +{ + using Deriv = multifluid::DerivativeOffset; + + auto flashModel = makeFlashModel( CO2SolubilitySpycherPruessTestFixture::flashContent ); + + auto [pressure, temperature, z_co2, expected_V, expected_x_co2, expected_y_wat] = GetParam(); + GEOS_UNUSED_VAR( expected_V, expected_x_co2, expected_y_wat ); + + StackArray< real64, 1, numComp > composition( numComp ); + composition[0] = z_co2; + composition[1] = 1.0 - z_co2; + StackArray< real64, 1, numComp > perturbedComposition( numComp ); + + StackArray< real64, 3, numPhase, LAYOUT_PHASE > phaseFrac( 1, 1, numPhase ); + StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPhaseFrac( 1, 1, numPhase, numDof ); + MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > + phaseFracAndDeriv { phaseFrac[0][0], dPhaseFrac[0][0] }; + + StackArray< real64, 4, numComp *numPhase, LAYOUT_PHASE_COMP > phaseCompFrac( 1, 1, numPhase, numComp ); + StackArray< real64, 5, numDof *numComp *numPhase, LAYOUT_PHASE_COMP_DC > dPhaseCompFrac( 1, 1, numPhase, numComp, numDof ); + MultiFluidVarSlice< real64, 2, USD_PHASE_COMP - 2, USD_PHASE_COMP_DC - 2 > + phaseCompFracAndDeriv { phaseCompFrac[0][0], dPhaseCompFrac[0][0] }; + + // 1) First compute the unperturbed parameters + + auto flashModelWrapper = flashModel->createKernelWrapper(); + + flashModelWrapper.compute( pressure, + temperature, + composition.toSliceConst(), + phaseFracAndDeriv, + phaseCompFracAndDeriv ); + + StackArray< real64, 3, numPhase, LAYOUT_PHASE > perturbedPhaseFrac( 1, 1, numPhase ); + StackArray< real64, 4, numDof *numPhase, LAYOUT_PHASE_DC > dPerturbedPhaseFrac( 1, 1, numPhase, numDof ); + MultiFluidVarSlice< real64, 1, USD_PHASE - 2, USD_PHASE_DC - 2 > + perturbedPhaseFracAndDeriv { perturbedPhaseFrac[0][0], dPerturbedPhaseFrac[0][0] }; + + StackArray< real64, 4, numComp *numPhase, LAYOUT_PHASE_COMP > perturbedPhaseCompFrac( 1, 1, numPhase, numComp ); + StackArray< real64, 5, numDof *numComp *numPhase, LAYOUT_PHASE_COMP_DC > dPerturbedPhaseCompFrac( 1, 1, numPhase, numComp, numDof ); + MultiFluidVarSlice< real64, 2, USD_PHASE_COMP - 2, USD_PHASE_COMP_DC - 2 > + perturbedPhaseCompFracAndDeriv { perturbedPhaseCompFrac[0][0], dPerturbedPhaseCompFrac[0][0] }; + + real64 numericalDerivative = 0.0; + real64 analyticalDerivative = 0.0; + + // 2) Check derivative with respect to pressure + real64 const dP = pertubation * pressure; + flashModelWrapper.compute( pressure + dP, + temperature, + composition.toSliceConst(), + perturbedPhaseFracAndDeriv, + perturbedPhaseCompFracAndDeriv ); + + for( integer i = 0; i < numPhase; ++i ) + { + numericalDerivative = (perturbedPhaseFracAndDeriv.value[i]-phaseFracAndDeriv.value[i])/dP; + analyticalDerivative = perturbedPhaseFracAndDeriv.derivs[i][Deriv::dP]; + checkRelativeError( numericalDerivative, analyticalDerivative, relTol, absTol ); + for( integer j = 0; j < numComp; ++j ) + { + numericalDerivative = (perturbedPhaseCompFracAndDeriv.value[i][j]-phaseCompFracAndDeriv.value[i][j])/dP; + analyticalDerivative = perturbedPhaseCompFracAndDeriv.derivs[i][j][Deriv::dP]; + checkRelativeError( numericalDerivative, analyticalDerivative, relTol, absTol ); + } + } + + // 3) Check derivative with respect to temperature + real64 const dT = pertubation * temperature; + flashModelWrapper.compute( pressure, + temperature + dT, + composition.toSliceConst(), + perturbedPhaseFracAndDeriv, + perturbedPhaseCompFracAndDeriv ); + + for( integer i = 0; i < numPhase; ++i ) + { + numericalDerivative = (perturbedPhaseFracAndDeriv.value[i]-phaseFracAndDeriv.value[i])/dT; + analyticalDerivative = perturbedPhaseFracAndDeriv.derivs[i][Deriv::dT]; + checkRelativeError( numericalDerivative, analyticalDerivative, relTol, absTol ); + + for( integer j = 0; j < numComp; ++j ) + { + numericalDerivative = (perturbedPhaseCompFracAndDeriv.value[i][j]-phaseCompFracAndDeriv.value[i][j])/dT; + analyticalDerivative = perturbedPhaseCompFracAndDeriv.derivs[i][j][Deriv::dT]; + checkRelativeError( numericalDerivative, analyticalDerivative, relTol, absTol ); + } + } + + // 4) Check derivative with respect to composition + for( integer ic=0; ic < numComp; ++ic ) + { + real64 const dC = pertubation; + for( integer k=0; k < numComp; ++k ) + perturbedComposition[k] = composition[k]; + perturbedComposition[ic] += dC; + + flashModelWrapper.compute( pressure, + temperature, + perturbedComposition.toSliceConst(), + perturbedPhaseFracAndDeriv, + perturbedPhaseCompFracAndDeriv ); + + for( integer i = 0; i < numPhase; ++i ) + { + numericalDerivative = (perturbedPhaseFracAndDeriv.value[i]-phaseFracAndDeriv.value[i])/dC; + analyticalDerivative = perturbedPhaseFracAndDeriv.derivs[i][Deriv::dC+ic]; + checkRelativeError( numericalDerivative, analyticalDerivative, relTol, absTol ); + + for( integer j = 0; j < numComp; ++j ) + { + numericalDerivative = (perturbedPhaseCompFracAndDeriv.value[i][j]-phaseCompFracAndDeriv.value[i][j])/dC; + analyticalDerivative = perturbedPhaseCompFracAndDeriv.derivs[i][j][Deriv::dC+ic]; + checkRelativeError( numericalDerivative, analyticalDerivative, relTol, absTol ); + } + } + } +} + +// Test data +std::vector< TestParam > generateTestData() +{ + return { + {1.00000000e+05, 1.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+05, 5.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+05, 8.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+05, 1.10000000e+02, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 1.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 5.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 8.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 1.10000000e+02, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {5.00000000e+06, 1.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {5.00000000e+06, 5.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {5.00000000e+06, 8.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {5.00000000e+06, 1.10000000e+02, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+07, 1.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+07, 5.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+07, 8.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+07, 1.10000000e+02, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+08, 1.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+08, 5.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+08, 8.00000000e+01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+08, 1.10000000e+02, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+05, 1.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+05, 5.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+05, 8.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+05, 1.10000000e+02, 1.00000000e-05, 1.00000000e-05, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 1.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+06, 5.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+06, 8.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+06, 1.10000000e+02, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {5.00000000e+06, 1.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {5.00000000e+06, 5.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {5.00000000e+06, 8.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {5.00000000e+06, 1.10000000e+02, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+07, 1.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+07, 5.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+07, 8.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+07, 1.10000000e+02, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+08, 1.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+08, 5.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+08, 8.00000000e+01, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+08, 1.10000000e+02, 1.00000000e-05, 0.00000000e+00, 1.00000000e-05, 0.00000000e+00}, + {1.00000000e+05, 1.00000000e+01, 3.00000000e-01, 3.03161030e-01, 8.58355906e-04, 1.23998968e-02}, + {1.00000000e+05, 5.00000000e+01, 3.00000000e-01, 3.42216891e-01, 3.10066097e-04, 1.23958952e-01}, + {1.00000000e+05, 8.00000000e+01, 3.00000000e-01, 5.66944658e-01, 1.24160030e-04, 4.70942662e-01}, + {1.00000000e+05, 1.10000000e+02, 3.00000000e-01, 3.00000000e-01, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 1.00000000e+01, 3.00000000e-01, 2.94948990e-01, 8.11504677e-03, 2.27331474e-03}, + {1.00000000e+06, 5.00000000e+01, 3.00000000e-01, 3.04861076e-01, 3.32772200e-03, 2.35330303e-02}, + {1.00000000e+06, 8.00000000e+01, 3.00000000e-01, 3.35758987e-01, 2.14620607e-03, 1.10747849e-01}, + {1.00000000e+06, 1.10000000e+02, 3.00000000e-01, 3.40957999e-01, 1.62203194e-03, 1.23261476e-01}, + {5.00000000e+06, 1.00000000e+01, 3.00000000e-01, 2.79742151e-01, 2.89079228e-02, 2.01367187e-03}, + {5.00000000e+06, 5.00000000e+01, 3.00000000e-01, 2.91378214e-01, 1.36880088e-02, 3.69909366e-03}, + {5.00000000e+06, 8.00000000e+01, 3.00000000e-01, 2.97134400e-01, 9.51668603e-03, 1.28674138e-02}, + {5.00000000e+06, 1.10000000e+02, 3.00000000e-01, 3.05337535e-01, 8.04059852e-03, 3.57736478e-02}, + {1.00000000e+07, 1.00000000e+01, 3.00000000e-01, 2.78859111e-01, 3.01769401e-02, 2.22670319e-03}, + {1.00000000e+07, 5.00000000e+01, 3.00000000e-01, 2.86981820e-01, 1.98935606e-02, 4.06398711e-03}, + {1.00000000e+07, 8.00000000e+01, 3.00000000e-01, 2.91881746e-01, 1.54455782e-02, 9.65816584e-03}, + {1.00000000e+07, 1.10000000e+02, 3.00000000e-01, 2.97372350e-01, 1.38649677e-02, 2.39237441e-02}, + {1.00000000e+08, 1.00000000e+01, 3.00000000e-01, 2.78859111e-01, 3.01769401e-02, 2.22670319e-03}, + {1.00000000e+08, 5.00000000e+01, 3.00000000e-01, 2.86981820e-01, 1.98935606e-02, 4.06398711e-03}, + {1.00000000e+08, 8.00000000e+01, 3.00000000e-01, 2.91881746e-01, 1.54455782e-02, 9.65816584e-03}, + {1.00000000e+08, 1.10000000e+02, 3.00000000e-01, 2.97372350e-01, 1.38649677e-02, 2.39237441e-02}, + {1.00000000e+05, 1.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+05, 5.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+05, 8.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+05, 1.10000000e+02, 9.99990000e-01, 9.99990000e-01, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 1.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+06, 5.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+06, 8.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+06, 1.10000000e+02, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {5.00000000e+06, 1.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {5.00000000e+06, 5.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {5.00000000e+06, 8.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {5.00000000e+06, 1.10000000e+02, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+07, 1.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+07, 5.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+07, 8.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+07, 1.10000000e+02, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+08, 1.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+08, 5.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+08, 8.00000000e+01, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+08, 1.10000000e+02, 9.99990000e-01, 1.00000000e+00, 0.00000000e+00, 1.00000000e-05}, + {1.00000000e+05, 1.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+05, 5.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+05, 8.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+05, 1.10000000e+02, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 1.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 5.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 8.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+06, 1.10000000e+02, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {5.00000000e+06, 1.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {5.00000000e+06, 5.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {5.00000000e+06, 8.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {5.00000000e+06, 1.10000000e+02, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+07, 1.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+07, 5.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+07, 8.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+07, 1.10000000e+02, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+08, 1.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+08, 5.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+08, 8.00000000e+01, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00}, + {1.00000000e+08, 1.10000000e+02, 1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00} + }; +} + +INSTANTIATE_TEST_SUITE_P( + CO2SolubilitySpycherPruessTest, + CO2SolubilitySpycherPruessTestFixture, + ::testing::ValuesIn( generateTestData() ) + ); + +int main( int argc, char * * argv ) +{ + ::testing::InitGoogleTest( &argc, argv ); + + geos::GeosxState state( geos::basicSetup( argc, argv ) ); + + int const result = RUN_ALL_TESTS(); + + geos::basicCleanup(); + + return result; +} diff --git a/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine.xml b/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine.xml index 7363117393c..12ec4e62d8d 100644 --- a/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine.xml +++ b/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine.xml @@ -46,6 +46,28 @@ outputPhaseComposition="1" baseline="testPVT_CO2Brine_testCo2BrineEzrokhiMixtureB.txt" logLevel="1" /> + + + + @@ -79,6 +107,13 @@ componentMolarWeight="{ 44e-3, 18e-3 }" phasePVTParaFiles="{ testPVT_data/carbonDioxidePVT.txt, testPVT_data/brinePVTEzrokhi.txt }" flashModelParaFile="testPVT_data/carbonDioxideFlash.txt" /> + diff --git a/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine_testCo2SpycherPruessBrinePhillipsMixtureA.txt b/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine_testCo2SpycherPruessBrinePhillipsMixtureA.txt new file mode 100644 index 00000000000..0a33d6dc22e --- /dev/null +++ b/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine_testCo2SpycherPruessBrinePhillipsMixtureA.txt @@ -0,0 +1,31 @@ +# column 1 = time +# column 2 = pressure +# column 3 = temperature +# column 4 = density +# column 5 = total compressibility +# columns 6-7 = phase fractions +# columns 8-9 = phase densities +# columns 10-11 = phase viscosities +# columns 12-13 = gas phase fractions [co2, water] +# columns 14-15 = water phase fractions [co2, water] +0.0000e+00 1.0000e+06 3.5000e+02 6.7780e+01 2.1337e-08 2.1773e-01 7.8227e-01 1.5581e+01 1.0033e+03 1.7476e-05 4.1330e-04 9.5589e-01 4.4112e-02 2.2320e-03 9.9777e-01 +5.0000e-02 6.4444e+06 3.5000e+02 4.2328e+02 1.3988e-07 1.9253e-01 8.0747e-01 1.2309e+02 1.0114e+03 1.8998e-05 4.1330e-04 9.8984e-01 1.0158e-02 1.1731e-02 9.8827e-01 +1.0000e-01 1.1889e+07 3.5000e+02 7.1407e+02 6.1644e-08 1.8180e-01 8.1820e-01 3.0497e+02 1.0173e+03 2.5353e-05 4.1330e-04 9.9079e-01 9.2115e-03 1.7035e-02 9.8297e-01 +1.5000e-01 1.7333e+07 3.5000e+02 8.8262e+02 1.8484e-08 1.7661e-01 8.2339e-01 5.3970e+02 1.0219e+03 4.0932e-05 4.1330e-04 9.8874e-01 1.1255e-02 1.9789e-02 9.8021e-01 +2.0000e-01 2.2778e+07 3.5000e+02 9.3903e+02 7.0636e-09 1.7314e-01 8.2686e-01 6.6850e+02 1.0260e+03 5.3806e-05 4.1330e-04 9.8725e-01 1.2749e-02 2.1623e-02 9.7838e-01 +2.5000e-01 2.8222e+07 3.5000e+02 9.6553e+02 3.7969e-09 1.7065e-01 8.2935e-01 7.4098e+02 1.0297e+03 6.3051e-05 4.1330e-04 9.8650e-01 1.3500e-02 2.2909e-02 9.7709e-01 +3.0000e-01 3.3667e+07 3.5000e+02 9.8236e+02 2.6667e-09 1.6883e-01 8.3117e-01 7.9054e+02 1.0333e+03 7.0528e-05 4.1330e-04 9.8604e-01 1.3961e-02 2.3841e-02 9.7616e-01 +3.5000e-01 3.9111e+07 3.5000e+02 9.9478e+02 2.0090e-09 1.6746e-01 8.3254e-01 8.2831e+02 1.0367e+03 7.7028e-05 4.1330e-04 9.8571e-01 1.4290e-02 2.4538e-02 9.7546e-01 +4.0000e-01 4.4556e+07 3.5000e+02 1.0047e+03 1.6650e-09 1.6641e-01 8.3359e-01 8.5890e+02 1.0400e+03 8.2914e-05 4.1330e-04 9.8545e-01 1.4547e-02 2.5069e-02 9.7493e-01 +4.5000e-01 5.0000e+07 3.5000e+02 1.0131e+03 5.1360e-11 1.6561e-01 8.3439e-01 8.8476e+02 1.0432e+03 8.8389e-05 4.1330e-04 9.8524e-01 1.4761e-02 2.5477e-02 9.7452e-01 +5.0000e-01 3.0000e+07 3.0400e+02 9.8816e+02 2.6607e-09 1.6844e-01 8.3156e-01 8.0121e+02 1.0372e+03 7.1971e-05 8.6631e-04 9.8918e-01 1.0819e-02 2.3786e-02 9.7621e-01 +5.5000e-01 1.0000e+07 2.5800e+02 6.6430e+02 8.2527e-08 1.8156e-01 8.1844e-01 2.5862e+02 1.0189e+03 2.2796e-05 1.9511e-03 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +6.0000e-01 1.0000e+07 2.7156e+02 6.6481e+02 8.2643e-08 1.8156e-01 8.1844e-01 2.5862e+02 1.0203e+03 2.2796e-05 1.9588e-03 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +6.5000e-01 1.0000e+07 2.8511e+02 6.6522e+02 8.2733e-08 1.8156e-01 8.1844e-01 2.5862e+02 1.0215e+03 2.2796e-05 1.3681e-03 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +7.0000e-01 1.0000e+07 2.9867e+02 6.6551e+02 8.2800e-08 1.8156e-01 8.1844e-01 2.5862e+02 1.0223e+03 2.2796e-05 9.7022e-04 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +7.5000e-01 1.0000e+07 3.1222e+02 6.6570e+02 8.2842e-08 1.8156e-01 8.1844e-01 2.5862e+02 1.0229e+03 2.2796e-05 7.3691e-04 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +8.0000e-01 1.0000e+07 3.2578e+02 6.6579e+02 8.2862e-08 1.8156e-01 8.1844e-01 2.5862e+02 1.0231e+03 2.2796e-05 5.8345e-04 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +8.5000e-01 1.0000e+07 3.3933e+02 6.6563e+02 8.2888e-08 1.8167e-01 8.1833e-01 2.5862e+02 1.0231e+03 2.2796e-05 4.7645e-04 9.9353e-01 6.4726e-03 1.6875e-02 9.8313e-01 +9.0000e-01 1.0000e+07 3.5289e+02 6.1186e+02 7.9623e-08 1.8498e-01 8.1502e-01 2.2287e+02 1.0133e+03 2.2020e-05 3.9764e-04 9.9066e-01 9.3377e-03 1.5456e-02 9.8454e-01 +9.5000e-01 1.0000e+07 3.6644e+02 5.8928e+02 7.7989e-08 1.8498e-01 8.1502e-01 2.0825e+02 1.0078e+03 2.1800e-05 3.4093e-04 9.9066e-01 9.3377e-03 1.5456e-02 9.8454e-01 +1.0000e+00 1.0000e+07 3.8000e+02 5.8907e+02 7.7935e-08 1.8498e-01 8.1502e-01 2.0825e+02 1.0070e+03 2.1860e-05 2.9672e-04 9.9066e-01 9.3377e-03 1.5456e-02 9.8454e-01 diff --git a/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine_testCo2SpycherPruessBrinePhillipsMixtureB.txt b/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine_testCo2SpycherPruessBrinePhillipsMixtureB.txt new file mode 100644 index 00000000000..762ffae5a75 --- /dev/null +++ b/src/coreComponents/unitTests/constitutiveTests/testPVT_CO2Brine_testCo2SpycherPruessBrinePhillipsMixtureB.txt @@ -0,0 +1,31 @@ +# column 1 = time +# column 2 = pressure +# column 3 = temperature +# column 4 = density +# column 5 = total compressibility +# columns 6-7 = phase fractions +# columns 8-9 = phase densities +# columns 10-11 = phase viscosities +# columns 12-13 = gas phase fractions [co2, water] +# columns 14-15 = water phase fractions [co2, water] +0.0000e+00 1.0000e+06 3.5000e+02 2.9256e+01 8.7493e-09 5.2521e-01 4.7479e-01 1.5581e+01 1.0033e+03 1.7476e-05 4.1330e-04 9.5589e-01 4.4112e-02 2.2320e-03 9.9777e-01 +5.0000e-02 6.4444e+06 3.5000e+02 2.1910e+02 1.8135e-07 5.0111e-01 4.9889e-01 1.2309e+02 1.0114e+03 1.8998e-05 4.1330e-04 9.8984e-01 1.0158e-02 1.1731e-02 9.8827e-01 +1.0000e-01 1.1889e+07 3.5000e+02 4.7218e+02 1.0761e-07 4.9425e-01 5.0575e-01 3.0497e+02 1.0173e+03 2.5353e-05 4.1330e-04 9.9079e-01 9.2115e-03 1.7035e-02 9.8297e-01 +1.5000e-01 1.7333e+07 3.5000e+02 7.1004e+02 3.9438e-08 4.9158e-01 5.0842e-01 5.3970e+02 1.0219e+03 4.0932e-05 4.1330e-04 9.8874e-01 1.1255e-02 1.9789e-02 9.8021e-01 +2.0000e-01 2.2778e+07 3.5000e+02 8.1301e+02 1.5784e-08 4.8984e-01 5.1016e-01 6.6850e+02 1.0260e+03 5.3806e-05 4.1330e-04 9.8725e-01 1.2749e-02 2.1623e-02 9.7838e-01 +2.5000e-01 2.8222e+07 3.5000e+02 8.6506e+02 8.4258e-09 4.8850e-01 5.1150e-01 7.4098e+02 1.0297e+03 6.3051e-05 4.1330e-04 9.8650e-01 1.3500e-02 2.2909e-02 9.7709e-01 +3.0000e-01 3.3667e+07 3.5000e+02 8.9875e+02 5.8443e-09 4.8750e-01 5.1250e-01 7.9054e+02 1.0333e+03 7.0528e-05 4.1330e-04 9.8604e-01 1.3961e-02 2.3841e-02 9.7616e-01 +3.5000e-01 3.9111e+07 3.5000e+02 9.2359e+02 4.2911e-09 4.8675e-01 5.1325e-01 8.2831e+02 1.0367e+03 7.7028e-05 4.1330e-04 9.8571e-01 1.4290e-02 2.4538e-02 9.7546e-01 +4.0000e-01 4.4556e+07 3.5000e+02 9.4329e+02 3.4823e-09 4.8618e-01 5.1382e-01 8.5890e+02 1.0400e+03 8.2914e-05 4.1330e-04 9.8545e-01 1.4547e-02 2.5069e-02 9.7493e-01 +4.5000e-01 5.0000e+07 3.5000e+02 9.5970e+02 2.8240e-11 4.8574e-01 5.1426e-01 8.8476e+02 1.0432e+03 8.8389e-05 4.1330e-04 9.8524e-01 1.4761e-02 2.5477e-02 9.7452e-01 +5.0000e-01 3.0000e+07 3.0400e+02 9.0721e+02 5.8529e-09 4.8643e-01 5.1357e-01 8.0121e+02 1.0372e+03 7.1971e-05 8.6631e-04 9.8918e-01 1.0819e-02 2.3786e-02 9.7621e-01 +5.5000e-01 1.0000e+07 2.5800e+02 4.1580e+02 1.3435e-07 4.9338e-01 5.0662e-01 2.5862e+02 1.0189e+03 2.2796e-05 1.9511e-03 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +6.0000e-01 1.0000e+07 2.7156e+02 4.1592e+02 1.3441e-07 4.9338e-01 5.0662e-01 2.5862e+02 1.0203e+03 2.2796e-05 1.9588e-03 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +6.5000e-01 1.0000e+07 2.8511e+02 4.1602e+02 1.3446e-07 4.9338e-01 5.0662e-01 2.5862e+02 1.0215e+03 2.2796e-05 1.3681e-03 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +7.0000e-01 1.0000e+07 2.9867e+02 4.1609e+02 1.3449e-07 4.9338e-01 5.0662e-01 2.5862e+02 1.0223e+03 2.2796e-05 9.7022e-04 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +7.5000e-01 1.0000e+07 3.1222e+02 4.1614e+02 1.3452e-07 4.9338e-01 5.0662e-01 2.5862e+02 1.0229e+03 2.2796e-05 7.3691e-04 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +8.0000e-01 1.0000e+07 3.2578e+02 4.1616e+02 1.3453e-07 4.9338e-01 5.0662e-01 2.5862e+02 1.0231e+03 2.2796e-05 5.8345e-04 9.9359e-01 6.4066e-03 1.6921e-02 9.8308e-01 +8.5000e-01 1.0000e+07 3.3933e+02 4.1611e+02 1.3454e-07 4.9346e-01 5.0654e-01 2.5862e+02 1.0231e+03 2.2796e-05 4.7645e-04 9.9353e-01 6.4726e-03 1.6875e-02 9.8313e-01 +9.0000e-01 1.0000e+07 3.5289e+02 3.6713e+02 1.2176e-07 4.9624e-01 5.0376e-01 2.2287e+02 1.0133e+03 2.2020e-05 3.9764e-04 9.9066e-01 9.3377e-03 1.5456e-02 9.8454e-01 +9.5000e-01 1.0000e+07 3.6644e+02 3.4688e+02 1.1663e-07 4.9624e-01 5.0376e-01 2.0825e+02 1.0078e+03 2.1800e-05 3.4093e-04 9.9066e-01 9.3377e-03 1.5456e-02 9.8454e-01 +1.0000e+00 1.0000e+07 3.8000e+02 3.4684e+02 1.1661e-07 4.9624e-01 5.0376e-01 2.0825e+02 1.0070e+03 2.1860e-05 2.9672e-04 9.9066e-01 9.3377e-03 1.5456e-02 9.8454e-01 diff --git a/src/coreComponents/unitTests/constitutiveTests/testPVT_data/carbonDioxideSpycherPruessFlash.txt b/src/coreComponents/unitTests/constitutiveTests/testPVT_data/carbonDioxideSpycherPruessFlash.txt new file mode 100644 index 00000000000..98f2cb18927 --- /dev/null +++ b/src/coreComponents/unitTests/constitutiveTests/testPVT_data/carbonDioxideSpycherPruessFlash.txt @@ -0,0 +1 @@ +FlashModel CO2Solubility 0.99e6 5.01e7 4.911e6 339 351 1 0 1.0e-10 SpycherPruess