LCOV - code coverage report
Current view: top level - src/materials - PorousFlowFluidStateBrineCO2.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 315 317 99.4 %
Date: 2017-11-16 18:50:22 Functions: 12 12 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /****************************************************************/
       2             : /* MOOSE - Multiphysics Object Oriented Simulation Environment  */
       3             : /*                                                              */
       4             : /*          All contents are licensed under LGPL V2.1           */
       5             : /*             See LICENSE for full restrictions                */
       6             : /****************************************************************/
       7             : 
       8             : #include "PorousFlowFluidStateBrineCO2.h"
       9             : #include "BrineFluidProperties.h"
      10             : #include "SinglePhaseFluidPropertiesPT.h"
      11             : #include "PorousFlowCapillaryPressure.h"
      12             : 
      13             : template <>
      14             : InputParameters
      15           8 : validParams<PorousFlowFluidStateBrineCO2>()
      16             : {
      17           8 :   InputParameters params = validParams<PorousFlowFluidStateFlashBase>();
      18          24 :   params.addRequiredParam<UserObjectName>("brine_fp", "The name of the user object for brine");
      19          24 :   params.addRequiredParam<UserObjectName>("co2_fp", "The name of the user object for CO2");
      20          24 :   params.addCoupledVar("xnacl", 0, "The salt mass fraction in the brine (kg/kg)");
      21          16 :   params.addClassDescription("Fluid state class for brine and CO2");
      22           8 :   return params;
      23             : }
      24             : 
      25          24 : PorousFlowFluidStateBrineCO2::PorousFlowFluidStateBrineCO2(const InputParameters & parameters)
      26             :   : PorousFlowFluidStateFlashBase(parameters),
      27             : 
      28          72 :     _xnacl(_nodal_material ? coupledNodalValue("xnacl") : coupledValue("xnacl")),
      29          48 :     _brine_fp(getUserObject<BrineFluidProperties>("brine_fp")),
      30          48 :     _co2_fp(getUserObject<SinglePhaseFluidPropertiesPT>("co2_fp")),
      31          24 :     _water_fp(_brine_fp.getComponent(BrineFluidProperties::WATER)),
      32          24 :     _Mh2o(_brine_fp.molarMassH2O()),
      33          24 :     _invMh2o(1.0 / _Mh2o),
      34          24 :     _Mco2(_co2_fp.molarMass()),
      35          24 :     _Mnacl(_brine_fp.molarMassNaCl()),
      36         216 :     _Rbar(_R * 10.0)
      37             : {
      38             :   // Check that the correct FluidProperties UserObjects have been provided
      39          48 :   if (_brine_fp.fluidName() != "brine")
      40           0 :     mooseError("Only a valid Brine FluidProperties UserObject can be provided in brine_fp");
      41             : 
      42          48 :   if (_co2_fp.fluidName() != "co2")
      43           0 :     mooseError("Only a valid CO2 FluidProperties UserObject can be provided in co2_fp");
      44          24 : }
      45             : 
      46             : void
      47        4267 : PorousFlowFluidStateBrineCO2::thermophysicalProperties()
      48             : {
      49             :   // The FluidProperty objects use temperature in K
      50        8534 :   Real Tk = _temperature[_qp] + _T_c2k;
      51        8534 :   Real pressure = _gas_porepressure[_qp];
      52        8534 :   Real xnacl = _xnacl[_qp];
      53             : 
      54             :   // Mass fraction of CO2 in liquid and H2O in gas phases
      55             :   Real X0, dX0_dp, dX0_dT, Y0, Y1, dY1_dp, dY1_dT;
      56        4267 :   massFractions(pressure, Tk, xnacl, X0, dX0_dp, dX0_dT, Y1, dY1_dp, dY1_dT);
      57        4264 :   Y0 = 1.0 - Y1;
      58             : 
      59             :   // Determine which phases are present based on the value of z
      60             :   bool is_liquid = false;
      61             :   bool is_gas = false;
      62             :   bool is_twophase = false;
      63             : 
      64        8528 :   if ((*_z[0])[_qp] <= X0)
      65             :   {
      66             :     // In this case, there is not enough CO2 to form a gas phase,
      67             :     // so only a liquid phase is present
      68             :     is_liquid = true;
      69             :   }
      70        2946 :   else if ((*_z[0])[_qp] > X0 && (*_z[0])[_qp] < Y0)
      71             :   {
      72             :     // Two phases are present
      73             :     is_twophase = true;
      74             :   }
      75             :   else // ((*_z[0])[_qp] >= Y0)
      76             :   {
      77             :     // In this case, there is not enough brine to form a liquid
      78             :     // phase, so only a gas phase is present
      79             :     is_gas = true;
      80             :   }
      81             : 
      82             :   // Material must provide the following properties
      83             :   Real vapor_mass_fraction = 0.0;
      84             :   Real gas_saturation, liquid_saturation;
      85             :   Real gas_density, liquid_density;
      86             :   Real gas_viscosity, liquid_viscosity;
      87             : 
      88             :   // And the following derivatives
      89             :   Real dX0_dz = 0.0, dY0_dp = 0.0, dY0_dT = 0.0, dY0_dz = 0.0;
      90             :   Real dgas_density_dp = 0.0, dgas_density_dT = 0.0, dgas_density_dz = 0.0;
      91             :   Real dliquid_density_dp = 0.0, dliquid_density_dT = 0.0, dliquid_density_dz = 0.0;
      92             :   Real dgas_viscosity_dp = 0.0, dgas_viscosity_dT = 0.0, dgas_viscosity_dz = 0.0;
      93        4264 :   Real dliquid_viscosity_dp = 0.0, dliquid_viscosity_dT = 0.0, dliquid_viscosity_dz = 0.0;
      94             : 
      95        4264 :   if (is_liquid)
      96             :   {
      97        1318 :     X0 = (*_z[0])[_qp];
      98             :     Y0 = 0.0;
      99             :     gas_saturation = 0.0;
     100        1318 :     dX0_dp = 0.0;
     101        1318 :     dX0_dT = 0.0;
     102             :     dX0_dz = 1.0;
     103             :     gas_density = 0.0;
     104             :     gas_viscosity = 1.0; // To guard against division by 0
     105             :   }
     106        2946 :   else if (is_twophase)
     107             :   {
     108             :     // Two phase are present. Set mass equilibrium constants used in the
     109             :     // calculation of vapor mass fraction
     110             :     Y0 = 1.0 - Y1;
     111        1538 :     std::vector<Real> Ki(_num_components);
     112        1538 :     Ki[0] = Y0 / X0;
     113        1538 :     Ki[1] = (1.0 - Y0) / (1.0 - X0);
     114        1538 :     vapor_mass_fraction = vaporMassFraction(Ki);
     115             : 
     116             :     // Derivatives of mass fractions wrt PorousFlow variables
     117        1538 :     dY0_dp = -dY1_dp;
     118        1538 :     dY0_dT = -dY1_dT;
     119             :   }
     120             :   else // if (is_gas)
     121             :   {
     122        1408 :     X0 = 0.0;
     123        1408 :     dX0_dp = 0.0;
     124        1408 :     dX0_dT = 0.0;
     125        1408 :     Y0 = (*_z[0])[_qp];
     126             :     dY0_dz = 1.0;
     127             :     vapor_mass_fraction = 1.0;
     128             :     gas_saturation = 1.0;
     129             :     liquid_density = 0.0;
     130        1408 :     liquid_viscosity = 1.0; // To guard against division by 0
     131             :   }
     132             : 
     133             :   // Update all remaining mass fractions and derivatives
     134             :   // Mass fraction of CO2 in gas and H2O in liquid phases
     135        4264 :   Real X1 = 1.0 - X0;
     136        4264 :   Y1 = 1.0 - Y0;
     137             : 
     138             :   // Calculate the gas density and viscosity in the single phase gas or two
     139             :   // phase region
     140        4264 :   if (is_gas || is_twophase)
     141             :   {
     142             :     Real co2_density, dco2_density_dp, dco2_density_dT;
     143        2946 :     _co2_fp.rho_dpT(pressure, Tk, co2_density, dco2_density_dp, dco2_density_dT);
     144             : 
     145             :     // Gas density is given by the CO2 density - no correction due to the small amount of
     146             :     // brine vapor is made
     147        2946 :     gas_density = co2_density;
     148        2946 :     dgas_density_dp = dco2_density_dp;
     149        2946 :     dgas_density_dT = dco2_density_dT;
     150             :     dgas_density_dz = 0.0;
     151             : 
     152             :     Real co2_viscosity, dco2_viscosity_drho, dco2_viscosity_dT;
     153        2946 :     _co2_fp.mu_drhoT_from_rho_T(
     154        2946 :         co2_density, Tk, dco2_density_dT, co2_viscosity, dco2_viscosity_drho, dco2_viscosity_dT);
     155             : 
     156             :     // Assume that the viscosity of the gas phase is a weighted sum of the
     157             :     // individual viscosities
     158        2946 :     gas_viscosity = co2_viscosity;
     159        2946 :     dgas_viscosity_dp = dco2_viscosity_drho * dco2_density_dp;
     160        2946 :     dgas_viscosity_dT = dco2_viscosity_dT;
     161             :     dgas_viscosity_dz = 0.0;
     162             :   }
     163             : 
     164             :   // Calculate the saturation in the two phase case using the vapor mass
     165             :   // fraction
     166        4264 :   if (is_twophase)
     167             :   {
     168             :     // Liquid density
     169             :     Real co2_partial_density, dco2_partial_density_dT;
     170        1538 :     partialDensityCO2(Tk, co2_partial_density, dco2_partial_density_dT);
     171             :     // Use old value of gas saturation to estimate liquid saturation
     172             :     Real liqsat = 1.0;
     173        1538 :     if (!_is_initqp)
     174        2988 :       liqsat -= _saturation_old[_qp][_gas_phase_number];
     175             : 
     176        1538 :     liquid_density =
     177        3076 :         1.0 / (X0 / co2_partial_density +
     178        1538 :                X1 / _brine_fp.rho(pressure - _pc_uo.capillaryPressure(liqsat), Tk, xnacl));
     179             : 
     180             :     // The gas saturation in the two phase case
     181        3076 :     gas_saturation = vapor_mass_fraction * liquid_density /
     182        1538 :                      (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
     183             :   }
     184             : 
     185             :   // Calculate the saturations and pressures for each phase
     186        4264 :   liquid_saturation = 1.0 - gas_saturation;
     187        4264 :   Real liquid_porepressure = pressure - _pc_uo.capillaryPressure(liquid_saturation);
     188             : 
     189             :   // Calculate liquid density and viscosity if in the two phase or single phase
     190             :   // liquid region, including a density correction due to the presence of dissolved
     191             :   // CO2
     192        4264 :   if (is_liquid || is_twophase)
     193             :   {
     194             :     Real brine_density, dbrine_density_dp, dbrine_density_dT, dbrine_density_dx;
     195             :     Real dliquid_viscosity_drho, dliquid_viscosity_dx;
     196             : 
     197        2856 :     _brine_fp.rho_dpTx(liquid_porepressure,
     198             :                        Tk,
     199             :                        xnacl,
     200             :                        brine_density,
     201             :                        dbrine_density_dp,
     202             :                        dbrine_density_dT,
     203        2856 :                        dbrine_density_dx);
     204             : 
     205             :     // The liquid density
     206             :     Real co2_partial_density, dco2_partial_density_dT;
     207        2856 :     partialDensityCO2(Tk, co2_partial_density, dco2_partial_density_dT);
     208        2856 :     liquid_density = 1.0 / (X0 / co2_partial_density + X1 / brine_density);
     209             : 
     210        2856 :     dliquid_density_dp =
     211        5712 :         (dX0_dp / brine_density + X1 * dbrine_density_dp / brine_density / brine_density -
     212        5712 :          dX0_dp / co2_partial_density) *
     213             :         liquid_density * liquid_density;
     214             : 
     215        2856 :     dliquid_density_dT =
     216        5712 :         (dX0_dT / brine_density + X1 * dbrine_density_dT / brine_density / brine_density -
     217        5712 :          dX0_dT / co2_partial_density +
     218        5712 :          X0 * dco2_partial_density_dT / co2_partial_density / co2_partial_density) *
     219             :         liquid_density * liquid_density;
     220             : 
     221        2856 :     dliquid_density_dz =
     222        2856 :         (dX0_dz / brine_density - dX0_dz / co2_partial_density) * liquid_density * liquid_density;
     223             : 
     224             :     // Liquid viscosity is just the brine viscosity
     225             :     // Note: brine viscosity (and derivatives) requires water density (and derivatives)
     226             :     Real water_density, dwater_density_dp, dwater_density_dT;
     227        2856 :     _water_fp.rho_dpT(liquid_porepressure, Tk, water_density, dwater_density_dp, dwater_density_dT);
     228             : 
     229        2856 :     _brine_fp.mu_drhoTx(water_density,
     230             :                         Tk,
     231             :                         xnacl,
     232             :                         dwater_density_dT,
     233             :                         liquid_viscosity,
     234             :                         dliquid_viscosity_drho,
     235             :                         dliquid_viscosity_dT,
     236        2856 :                         dliquid_viscosity_dx);
     237             : 
     238             :     // The derivative of viscosity wrt pressure is given by the chain rule
     239        2856 :     dliquid_viscosity_dp = dliquid_viscosity_drho * dwater_density_dp;
     240             :   }
     241             : 
     242             :   // Save properties in FluidStateProperties vector
     243        4264 :   auto & liquid = _fsp[_aqueous_phase_number];
     244        4264 :   auto & gas = _fsp[_gas_phase_number];
     245             : 
     246        4264 :   liquid.saturation = liquid_saturation;
     247        4264 :   gas.saturation = gas_saturation;
     248        4264 :   liquid.pressure = liquid_porepressure;
     249        8528 :   gas.pressure = _gas_porepressure[_qp];
     250             : 
     251        8528 :   liquid.mass_fraction[_aqueous_fluid_component] = 1.0 - X0;
     252        8528 :   liquid.mass_fraction[_gas_fluid_component] = X0;
     253        4264 :   gas.mass_fraction[_aqueous_fluid_component] = 1.0 - Y0;
     254        4264 :   gas.mass_fraction[_gas_fluid_component] = Y0;
     255             : 
     256        4264 :   liquid.fluid_density = liquid_density;
     257        4264 :   gas.fluid_density = gas_density;
     258        4264 :   liquid.fluid_viscosity = liquid_viscosity;
     259        4264 :   gas.fluid_viscosity = gas_viscosity;
     260             : 
     261             :   // Derivatives wrt PorousFlow variables are required by the kernels
     262             :   // Note: these don't need to be stateful so don't calculate them in
     263             :   // initQpStatefulProperties
     264        4264 :   if (!_is_initqp)
     265             :   {
     266             :     // Derivative of gas saturation wrt variables
     267             :     Real ds_dp = 0.0, ds_dT = 0.0, ds_dz = 0.0;
     268        4160 :     if (is_twophase)
     269             :     {
     270        1494 :       Real K0 = Y0 / X0;
     271        1494 :       Real K1 = (1.0 - Y0) / (1.0 - X0);
     272        1494 :       Real dv_dz = (K1 - K0) / ((K0 - 1.0) * (K1 - 1.0));
     273             : 
     274        1494 :       ds_dz = gas_density * liquid_density * dv_dz;
     275        1494 :       ds_dz /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
     276             :                (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
     277             : 
     278        1494 :       Real dK0_dp = (-Y0 * dX0_dp - X0 * dY1_dp) / X0 / X0;
     279        1494 :       Real dK0_dT = (-Y0 * dX0_dT - X0 * dY1_dT) / X0 / X0;
     280             : 
     281        1494 :       Real dK1_dp = (Y1 * dX0_dp + (1.0 - X0) * dY1_dp) / (1.0 - X0) / (1.0 - X0);
     282        1494 :       Real dK1_dT = (Y1 * dX0_dT + (1.0 - X0) * dY1_dT) / (1.0 - X0) / (1.0 - X0);
     283             : 
     284        2988 :       Real dv_dp = (*_z[0])[_qp] * dK1_dp / (K1 - 1.0) / (K1 - 1.0) +
     285        1494 :                    (1.0 - (*_z[0])[_qp]) * dK0_dp / (K0 - 1.0) / (K0 - 1.0);
     286             : 
     287        2988 :       ds_dp = gas_density * liquid_density * dv_dp +
     288        2988 :               vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
     289        1494 :                   (gas_density * dliquid_density_dp - dgas_density_dp * liquid_density);
     290        1494 :       ds_dp /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
     291             :                (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
     292             : 
     293        1494 :       Real dv_dT = (*_z[0])[_qp] * dK1_dT / (K1 - 1.0) / (K1 - 1.0) +
     294        1494 :                    (1.0 - (*_z[0])[_qp]) * dK0_dT / (K0 - 1.0) / (K0 - 1.0);
     295             : 
     296        2988 :       ds_dT = gas_density * liquid_density * dv_dT +
     297        1494 :               vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
     298        1494 :                   (gas_density * dliquid_density_dT - dgas_density_dT * liquid_density);
     299        1494 :       ds_dT /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
     300             :                (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
     301             :     }
     302             : 
     303             :     // Save derivatives in FluidStateProperties vector
     304        4160 :     liquid.dsaturation_dp = -ds_dp;
     305        4160 :     liquid.dsaturation_dT = -ds_dT;
     306        4160 :     liquid.dsaturation_dz = -ds_dz;
     307        4160 :     gas.dsaturation_dp = ds_dp;
     308        4160 :     gas.dsaturation_dT = ds_dT;
     309        4160 :     gas.dsaturation_dz = ds_dz;
     310             : 
     311        4160 :     liquid.dmass_fraction_dp[_aqueous_fluid_component] = -dX0_dp;
     312        4160 :     liquid.dmass_fraction_dp[_gas_fluid_component] = dX0_dp;
     313        4160 :     liquid.dmass_fraction_dT[_aqueous_fluid_component] = -dX0_dT;
     314        4160 :     liquid.dmass_fraction_dT[_gas_fluid_component] = dX0_dT;
     315        4160 :     liquid.dmass_fraction_dz[_aqueous_fluid_component] = -dX0_dz;
     316        4160 :     liquid.dmass_fraction_dz[_gas_fluid_component] = dX0_dz;
     317        4160 :     gas.dmass_fraction_dp[_aqueous_fluid_component] = -dY0_dp;
     318        4160 :     gas.dmass_fraction_dp[_gas_fluid_component] = dY0_dp;
     319        4160 :     gas.dmass_fraction_dT[_aqueous_fluid_component] = -dY0_dT;
     320        4160 :     gas.dmass_fraction_dT[_gas_fluid_component] = dY0_dT;
     321        4160 :     gas.dmass_fraction_dz[_aqueous_fluid_component] = -dY0_dz;
     322        4160 :     gas.dmass_fraction_dz[_gas_fluid_component] = dY0_dz;
     323             : 
     324        4160 :     liquid.dfluid_density_dp = dliquid_density_dp;
     325        4160 :     liquid.dfluid_density_dT = dliquid_density_dT;
     326        4160 :     liquid.dfluid_density_dz = dliquid_density_dz;
     327        4160 :     gas.dfluid_density_dp = dgas_density_dp;
     328        4160 :     gas.dfluid_density_dT = dgas_density_dT;
     329        4160 :     gas.dfluid_density_dz = dgas_density_dz;
     330             : 
     331        4160 :     liquid.dfluid_viscosity_dp = dliquid_viscosity_dp;
     332        4160 :     liquid.dfluid_viscosity_dT = dliquid_viscosity_dT;
     333        4160 :     liquid.dfluid_viscosity_dz = dliquid_viscosity_dz;
     334        4160 :     gas.dfluid_viscosity_dp = dgas_viscosity_dp;
     335        4160 :     gas.dfluid_viscosity_dT = dgas_viscosity_dT;
     336        4160 :     gas.dfluid_viscosity_dz = dgas_viscosity_dz;
     337             :   }
     338        4264 : }
     339             : 
     340             : void
     341        4267 : PorousFlowFluidStateBrineCO2::massFractions(Real pressure,
     342             :                                             Real temperature,
     343             :                                             Real xnacl,
     344             :                                             Real & xco2l,
     345             :                                             Real & dxco2l_dp,
     346             :                                             Real & dxco2l_dT,
     347             :                                             Real & xh2og,
     348             :                                             Real & dxh2og_dp,
     349             :                                             Real & dxh2og_dT) const
     350             : {
     351             :   // Pressure in bar
     352        4267 :   Real pbar = pressure * 1.0e-5;
     353             :   // Pressure minus 1 bar
     354        4267 :   Real delta_pbar = pbar - 1.0;
     355             : 
     356             :   // Average partial molar volumes (cm^3/mol) as given by Sypcher, Pruess and Ennis-King (2003)
     357             :   const Real vCO2 = 32.6;
     358             :   const Real vH2O = 18.1;
     359             : 
     360             :   // NaCl molality (mol/kg)
     361        4267 :   Real mnacl = xnacl / (1.0 - xnacl) / _Mnacl;
     362             : 
     363             :   // Equilibrium constants
     364             :   Real K0H2O, dK0H2O_dT, K0CO2, dK0CO2_dT;
     365        4267 :   equilibriumConstantH2O(temperature, K0H2O, dK0H2O_dT);
     366        4267 :   equilibriumConstantCO2(temperature, K0CO2, dK0CO2_dT);
     367             : 
     368             :   // Fugacity coefficients
     369             :   Real phiH2O, dphiH2O_dp, dphiH2O_dT;
     370             :   Real phiCO2, dphiCO2_dp, dphiCO2_dT;
     371        4267 :   fugacityCoefficientCO2(pressure, temperature, phiCO2, dphiCO2_dp, dphiCO2_dT);
     372        4264 :   fugacityCoefficientH2O(pressure, temperature, phiH2O, dphiH2O_dp, dphiH2O_dT);
     373             : 
     374             :   // Activity coefficient
     375             :   Real gamma, dgamma_dp, dgamma_dT;
     376        4264 :   activityCoefficient(pressure, temperature, xnacl, gamma, dgamma_dp, dgamma_dT);
     377             : 
     378        4264 :   Real A = K0H2O / (phiH2O * pbar) * std::exp(delta_pbar * vH2O / (_Rbar * temperature));
     379             :   Real B =
     380        4264 :       phiCO2 * pbar / (_invMh2o * K0CO2) * std::exp(-delta_pbar * vCO2 / (_Rbar * temperature));
     381             : 
     382        8528 :   Real dA_dp = (-1.0e-5 * K0H2O / pbar + 1.0e-5 * vH2O * K0H2O / (_Rbar * temperature) -
     383        8528 :                 K0H2O * dphiH2O_dp / phiH2O) *
     384        8528 :                std::exp(delta_pbar * vH2O / (_Rbar * temperature)) / (pbar * phiH2O);
     385        8528 :   Real dB_dp = (1.0e-5 * phiCO2 + pbar * dphiCO2_dp -
     386        8528 :                 1.0e-5 * vCO2 * pbar * phiCO2 / (_Rbar * temperature)) *
     387        8528 :                std::exp(-delta_pbar * vCO2 / (_Rbar * temperature)) / (_invMh2o * K0CO2);
     388             : 
     389        8528 :   Real dA_dT = (dK0H2O_dT - dphiH2O_dT * K0H2O / phiH2O -
     390        8528 :                 delta_pbar * vH2O * K0H2O / (_Rbar * temperature * temperature)) *
     391        8528 :                std::exp(delta_pbar * vH2O / (_Rbar * temperature)) / (pbar * phiH2O);
     392        8528 :   Real dB_dT = (-pbar * phiCO2 * dK0CO2_dT / K0CO2 + pbar * dphiCO2_dT +
     393        8528 :                 delta_pbar * vCO2 * pbar * phiCO2 / (_Rbar * temperature * temperature)) *
     394        8528 :                std::exp(-delta_pbar * vCO2 / (_Rbar * temperature)) / (_invMh2o * K0CO2);
     395             : 
     396             :   // The mole fraction of H2O in the CO2-rich gas phase is then
     397        4264 :   Real yH2O = (1.0 - B) / (1.0 / A - B);
     398             :   // The mole fraction of CO2 in the H2O-rich liquid phase is then (note: no salinty effect)
     399        4264 :   Real xCO2 = B * (1.0 - yH2O);
     400             :   // The molality of CO2 in the H2O-rich liquid phase is (note: no salinty effect)
     401        4264 :   Real mCO2 = xCO2 * _invMh2o / (1.0 - xCO2);
     402             : 
     403             :   // The molality of CO2 in brine is then given by
     404        4264 :   Real mCO2b = mCO2 / gamma;
     405             :   // The mole fraction of CO2 in brine is then
     406        4264 :   Real xCO2b = mCO2b / (2.0 * mnacl + _invMh2o + mCO2b);
     407             :   // The mole fraction of H2O in the CO2-rich gas phase corrected for NaCl mole fraction is
     408        4264 :   Real yH2Ob = A * (1.0 - xCO2b - 2.0 * mnacl / (2.0 * mnacl + _invMh2o + mCO2b));
     409             : 
     410             :   // Convert the mole fractions to mass fractions and then update referenced values
     411             :   // The mass fraction of H2O in gas (assume no salt in gas phase)
     412        4264 :   xh2og = yH2Ob * _Mh2o / (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2);
     413             : 
     414             :   // The number of moles of CO2 in 1kg of H2O
     415        4264 :   Real nco2 = xCO2b * (2.0 * mnacl + _invMh2o) / (1.0 - xCO2b);
     416             : 
     417             :   // The mass fraction of CO2 in liquid
     418        4264 :   xco2l = nco2 * _Mco2 / (1.0 + mnacl * _Mnacl + nco2 * _Mco2);
     419             : 
     420             :   // The derivatives of the mass fractions wrt pressure
     421        4264 :   Real dyH2O_dp = ((1.0 - B) * dA_dp + (A - 1.0) * A * dB_dp) / (1.0 - A * B) / (1.0 - A * B);
     422        4264 :   Real dxCO2_dp = dB_dp * (1.0 - yH2O) - B * dyH2O_dp;
     423             : 
     424        4264 :   Real dmCO2_dp = _invMh2o * dxCO2_dp / (1.0 - xCO2) / (1.0 - xCO2);
     425        4264 :   Real dmCO2b_dp = dmCO2_dp / gamma - mCO2 * dgamma_dp / gamma / gamma;
     426        4264 :   Real dxCO2b_dp = (2.0 * mnacl + _invMh2o) * dmCO2b_dp / (2.0 * mnacl + _invMh2o + mCO2b) /
     427        4264 :                    (2.0 * mnacl + _invMh2o + mCO2b);
     428             : 
     429        8528 :   Real dyH2Ob_dp = (1.0 - xCO2b - 2.0 * mnacl / (2.0 * mnacl + _invMh2o + mCO2b)) * dA_dp -
     430        4264 :                    A * dxCO2b_dp +
     431        4264 :                    2.0 * A * mnacl * dmCO2b_dp / (2.0 * mnacl + _invMh2o + mCO2b) /
     432        4264 :                        (2.0 * mnacl + _invMh2o + mCO2b);
     433             : 
     434        4264 :   dxh2og_dp = _Mco2 * _Mh2o * dyH2Ob_dp / (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2) /
     435             :               (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2);
     436             : 
     437        4264 :   Real dnco2_dp = dxCO2b_dp * (2.0 * mnacl + _invMh2o) / (1.0 - xCO2b) / (1.0 - xCO2b);
     438             : 
     439        4264 :   dxco2l_dp = (1.0 + mnacl * _Mnacl) * _Mco2 * dnco2_dp / (1.0 + mnacl * _Mnacl + nco2 * _Mco2) /
     440             :               (1.0 + mnacl * _Mnacl + nco2 * _Mco2);
     441             : 
     442             :   // The derivatives of the mass fractions wrt temperature
     443        4264 :   Real dyH2O_dT = ((1.0 - B) * dA_dT + (A - 1.0) * A * dB_dT) / (1.0 - A * B) / (1.0 - A * B);
     444        4264 :   Real dxCO2_dT = dB_dT * (1.0 - yH2O) - B * dyH2O_dT;
     445             : 
     446        4264 :   Real dmCO2_dT = _invMh2o * dxCO2_dT / (1.0 - xCO2) / (1.0 - xCO2);
     447        4264 :   Real dmCO2b_dT = dmCO2_dT / gamma - mCO2 * dgamma_dT / gamma / gamma;
     448        4264 :   Real dxCO2b_dT = (2.0 * mnacl + _invMh2o) * dmCO2b_dT / (2.0 * mnacl + _invMh2o + mCO2b) /
     449        4264 :                    (2.0 * mnacl + _invMh2o + mCO2b);
     450             : 
     451        8528 :   Real dyH2Ob_dT = (1.0 - xCO2b - 2.0 * mnacl / (2.0 * mnacl + _invMh2o + mCO2b)) * dA_dT -
     452        4264 :                    A * dxCO2b_dT +
     453        4264 :                    2.0 * A * mnacl * dmCO2b_dT / (2.0 * mnacl + _invMh2o + mCO2b) /
     454        4264 :                        (2.0 * mnacl + _invMh2o + mCO2b);
     455             : 
     456        4264 :   dxh2og_dT = _Mco2 * _Mh2o * dyH2Ob_dT / (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2) /
     457             :               (yH2Ob * _Mh2o + (1.0 - yH2Ob) * _Mco2);
     458             : 
     459        4264 :   Real dnco2_dT = dxCO2b_dT * (2.0 * mnacl + _invMh2o) / (1.0 - xCO2b) / (1.0 - xCO2b);
     460             : 
     461        4264 :   dxco2l_dT = (1.0 + mnacl * _Mnacl) * _Mco2 * dnco2_dT / (1.0 + mnacl * _Mnacl + nco2 * _Mco2) /
     462             :               (1.0 + mnacl * _Mnacl + nco2 * _Mco2);
     463        4264 : }
     464             : 
     465             : void
     466        4267 : PorousFlowFluidStateBrineCO2::fugacityCoefficientCO2(
     467             :     Real pressure, Real temperature, Real & fco2, Real & dfco2_dp, Real & dfco2_dT) const
     468             : {
     469             :   // Need pressure in bar
     470        4267 :   Real pbar = pressure * 1.0e-5;
     471             :   // CO2 density and derivatives wrt pressure and temperature
     472             :   Real gas_density, dgas_density_dp, dgas_density_dT;
     473        4267 :   _co2_fp.rho_dpT(pressure, temperature, gas_density, dgas_density_dp, dgas_density_dT);
     474             :   // Molar volume in cm^3/mol
     475        4264 :   Real V = _Mco2 / gas_density * 1.0e6;
     476             : 
     477             :   // Redlich-Kwong parameters
     478        4264 :   Real aCO2 = 7.54e7 - 4.13e4 * temperature;
     479             :   Real bCO2 = 27.8;
     480             : 
     481        8528 :   Real term1 = std::log(V / (V - bCO2)) + bCO2 / (V - bCO2) -
     482        4264 :                2.0 * aCO2 / (_Rbar * std::pow(temperature, 1.5) * bCO2) * std::log((V + bCO2) / V) +
     483        8528 :                aCO2 / (_Rbar * std::pow(temperature, 1.5) * bCO2) *
     484        8528 :                    (std::log((V + bCO2) / V) - bCO2 / (V + bCO2));
     485             : 
     486        4264 :   Real lnPhiCO2 = term1 - std::log(pbar * V / (_Rbar * temperature));
     487        4264 :   fco2 = std::exp(lnPhiCO2);
     488             : 
     489             :   // The derivative of the fugacity coefficient wrt pressure
     490        4264 :   Real dV_dp = -_Mco2 / gas_density / gas_density * dgas_density_dp * 1.0e6;
     491             :   Real dterm1_dV =
     492        8528 :       (bCO2 * (bCO2 - 2.0 * V) / (bCO2 - V) / (bCO2 - V) +
     493        4264 :        aCO2 * (bCO2 + 2.0 * V) / (_Rbar * std::pow(temperature, 1.5) * (bCO2 + V) * (bCO2 + V))) /
     494        4264 :       V;
     495        4264 :   dfco2_dp = (dterm1_dV * dV_dp - 1.0e-5 / pbar - 1.0 / V * dV_dp) * fco2;
     496             : 
     497             :   // The derivative of the fugacity coefficient wrt temperature
     498        4264 :   Real dV_dT = -_Mco2 / gas_density / gas_density * dgas_density_dT * 1.0e6;
     499        8528 :   Real dterm1_dT = 3.0 * aCO2 * _Rbar * std::sqrt(temperature) * bCO2 * std::log((V + bCO2) / V) /
     500        8528 :                        (_Rbar * std::pow(temperature, 1.5) * bCO2) /
     501        8528 :                        (_Rbar * std::pow(temperature, 1.5) * bCO2) +
     502        8528 :                    8.26e4 / (_Rbar * std::pow(temperature, 1.5) * bCO2) * std::log((V + bCO2) / V) -
     503        8528 :                    1.5 * aCO2 * _Rbar * std::sqrt(temperature) * bCO2 *
     504        8528 :                        (std::log((V + bCO2) / V) - bCO2 / (V + bCO2)) /
     505        8528 :                        (_Rbar * std::pow(temperature, 1.5) * bCO2) /
     506        4264 :                        (_Rbar * std::pow(temperature, 1.5) * bCO2) -
     507        8528 :                    4.13e4 / (_Rbar * std::pow(temperature, 1.5) * bCO2) *
     508        8528 :                        (std::log((V + bCO2) / V) - bCO2 / (V + bCO2));
     509        4264 :   dfco2_dT = (dterm1_dT + dterm1_dV * dV_dT - dV_dT / V + 1.0 / temperature) * fco2;
     510        4264 : }
     511             : 
     512             : void
     513        4264 : PorousFlowFluidStateBrineCO2::fugacityCoefficientH2O(
     514             :     Real pressure, Real temperature, Real & fh2o, Real & dfh2o_dp, Real & dfh2o_dT) const
     515             : {
     516             :   // Need pressure in bar
     517        4264 :   Real pbar = pressure * 1.0e-5;
     518             :   // CO2 density and derivatives wrt pressure and temperature
     519             :   Real gas_density, dgas_density_dp, dgas_density_dT;
     520        4264 :   _co2_fp.rho_dpT(pressure, temperature, gas_density, dgas_density_dp, dgas_density_dT);
     521             :   // Molar volume in cm^3/mol
     522        4264 :   Real V = _Mco2 / gas_density * 1.0e6;
     523             : 
     524             :   // Redlich-Kwong parameters
     525        4264 :   Real aCO2 = 7.54e7 - 4.13e4 * temperature;
     526             :   Real bCO2 = 27.8;
     527             :   Real aCO2H2O = 7.89e7;
     528             :   Real bH2O = 18.18;
     529             : 
     530             :   Real term1 =
     531        8528 :       std::log(V / (V - bCO2)) + bH2O / (V - bCO2) -
     532        4264 :       2.0 * aCO2H2O / (_Rbar * std::pow(temperature, 1.5) * bCO2) * std::log((V + bCO2) / V) +
     533        8528 :       aCO2 * bH2O / (_Rbar * std::pow(temperature, 1.5) * bCO2 * bCO2) *
     534        8528 :           (std::log((V + bCO2) / V) - bCO2 / (V + bCO2));
     535             : 
     536        4264 :   Real lnPhiH2O = term1 - std::log(pbar * V / (_Rbar * temperature));
     537        4264 :   fh2o = std::exp(lnPhiH2O);
     538             : 
     539             :   // The derivative of the fugacity coefficient wrt pressure
     540        4264 :   Real dV_dp = -_Mco2 / gas_density / gas_density * dgas_density_dp * 1.0e6;
     541        8528 :   Real dterm1_dV = ((bCO2 * bCO2 - (bCO2 + bH2O) * V) / (bCO2 - V) / (bCO2 - V) +
     542        8528 :                     (2.0 * aCO2H2O * (bCO2 + V) - aCO2 * bH2O) /
     543        4264 :                         (_Rbar * std::pow(temperature, 1.5) * (bCO2 + V) * (bCO2 + V))) /
     544        4264 :                    V;
     545        4264 :   dfh2o_dp = (dterm1_dV * dV_dp - 1.0e-5 / pbar - dV_dp / V) * fh2o;
     546             : 
     547             :   // The derivative of the fugacity coefficient wrt temperature
     548        4264 :   Real dV_dT = -_Mco2 / gas_density / gas_density * dgas_density_dT * 1.0e6;
     549        8528 :   Real dterm1_dT = 3.0 * _Rbar * std::sqrt(temperature) * bCO2 * aCO2H2O *
     550       12792 :                        std::log((V + bCO2) / V) / (_Rbar * std::pow(temperature, 1.5) * bCO2) /
     551        8528 :                        (_Rbar * std::pow(temperature, 1.5) * bCO2) -
     552        8528 :                    1.5 * aCO2 * bH2O * _Rbar * std::sqrt(temperature) * bCO2 * bCO2 *
     553        8528 :                        (std::log((V + bCO2) / V) - bCO2 / (V + bCO2)) /
     554        8528 :                        (_Rbar * std::pow(temperature, 1.5) * bCO2 * bCO2) /
     555        4264 :                        (_Rbar * std::pow(temperature, 1.5) * bCO2 * bCO2) -
     556        8528 :                    4.13e4 * bH2O * (std::log((V + bCO2) / V) - bCO2 / (V + bCO2)) /
     557        8528 :                        (_Rbar * std::pow(temperature, 1.5) * bCO2 * bCO2);
     558        4264 :   dfh2o_dT = (dterm1_dT + dterm1_dV * dV_dT - dV_dT / V + 1.0 / temperature) * fh2o;
     559        4264 : }
     560             : 
     561             : void
     562        4264 : PorousFlowFluidStateBrineCO2::activityCoefficient(Real pressure,
     563             :                                                   Real temperature,
     564             :                                                   Real xnacl,
     565             :                                                   Real & gamma,
     566             :                                                   Real & dgamma_dp,
     567             :                                                   Real & dgamma_dT) const
     568             : {
     569             :   // Need pressure in bar
     570        4264 :   Real pbar = pressure * 1.0e-5;
     571             :   // Need NaCl molality (mol/kg)
     572        4264 :   Real mnacl = xnacl / (1.0 - xnacl) / _Mnacl;
     573             : 
     574        8528 :   Real lambda = -0.411370585 + 6.07632013e-4 * temperature + 97.5347708 / temperature -
     575        8528 :                 0.0237622469 * pbar / temperature + 0.0170656236 * pbar / (630.0 - temperature) +
     576        4264 :                 1.41335834e-5 * temperature * std::log(pbar);
     577             : 
     578        4264 :   Real xi = 3.36389723e-4 - 1.9829898e-5 * temperature + 2.12220830e-3 * pbar / temperature -
     579        4264 :             5.24873303e-3 * pbar / (630.0 - temperature);
     580             : 
     581        4264 :   gamma = std::exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
     582             : 
     583             :   // Derivative wrt pressure
     584             :   Real dlambda_dp, dxi_dp;
     585        8528 :   dlambda_dp = -0.0237622469 / temperature + 0.0170656236 / (630.0 - temperature) +
     586        4264 :                1.41335834e-5 * temperature / pbar;
     587        4264 :   dxi_dp = 2.12220830e-3 / temperature - 5.24873303e-3 / (630.0 - temperature);
     588       12792 :   dgamma_dp = (2.0 * mnacl * dlambda_dp + mnacl * mnacl * dxi_dp) *
     589        8528 :               std::exp(2.0 * lambda * mnacl + xi * mnacl * mnacl) * 1.0e-5;
     590             : 
     591             :   // Derivative wrt temperature
     592             :   Real dlambda_dT, dxi_dT;
     593       12792 :   dlambda_dT = 6.07632013e-4 - 97.5347708 / temperature / temperature +
     594        8528 :                0.0237622469 * pbar / temperature / temperature +
     595        4264 :                0.0170656236 * pbar / (630.0 - temperature) / (630.0 - temperature) +
     596        4264 :                1.41335834e-5 * std::log(pbar);
     597        8528 :   dxi_dT = -1.9829898e-5 - 2.12220830e-3 * pbar / temperature / temperature -
     598        4264 :            5.24873303e-3 * pbar / (630.0 - temperature) / (630.0 - temperature);
     599        8528 :   dgamma_dT = (2.0 * mnacl * dlambda_dT + mnacl * mnacl * dxi_dT) *
     600        4264 :               std::exp(2.0 * lambda * mnacl + xi * mnacl * mnacl);
     601        4264 : }
     602             : 
     603             : void
     604        4267 : PorousFlowFluidStateBrineCO2::equilibriumConstantH2O(Real temperature,
     605             :                                                      Real & kh2o,
     606             :                                                      Real & dkh2o_dT) const
     607             : {
     608             :   // Uses temperature in Celcius
     609        4267 :   Real Tc = temperature - _T_c2k;
     610             : 
     611        4267 :   Real logK0H2O = -2.209 + 3.097e-2 * Tc - 1.098e-4 * Tc * Tc + 2.048e-7 * Tc * Tc * Tc;
     612        4267 :   Real dlogK0H2O = 3.097e-2 - 2.196e-4 * Tc + 6.144e-7 * Tc * Tc;
     613             : 
     614        4267 :   kh2o = std::pow(10.0, logK0H2O);
     615        4267 :   dkh2o_dT = std::log(10.0) * dlogK0H2O * kh2o;
     616        4267 : }
     617             : 
     618             : void
     619        4267 : PorousFlowFluidStateBrineCO2::equilibriumConstantCO2(Real temperature,
     620             :                                                      Real & kco2,
     621             :                                                      Real & dkco2_dT) const
     622             : {
     623             :   // Uses temperature in Celcius
     624        4267 :   Real Tc = temperature - _T_c2k;
     625             : 
     626        4267 :   Real logK0CO2 = 1.189 + 1.304e-2 * Tc - 5.446e-5 * Tc * Tc;
     627        4267 :   Real dlogK0CO2 = 1.304e-2 - 1.0892e-4 * Tc;
     628             : 
     629        4267 :   kco2 = std::pow(10.0, logK0CO2);
     630        4267 :   dkco2_dT = std::log(10.0) * dlogK0CO2 * kco2;
     631        4267 : }
     632             : 
     633             : void
     634        4394 : PorousFlowFluidStateBrineCO2::partialDensityCO2(Real temperature,
     635             :                                                 Real & partial_density,
     636             :                                                 Real & dpartial_density_dT) const
     637             : {
     638             :   // This correlation uses temperature in C
     639        4394 :   Real Tc = temperature - _T_c2k;
     640             :   // The parial molar volume
     641        4394 :   Real V = 37.51 - 9.585e-2 * Tc + 8.74e-4 * Tc * Tc - 5.044e-7 * Tc * Tc * Tc;
     642        4394 :   Real dV_dT = -9.585e-2 + 1.748e-3 * Tc - 1.5132e-6 * Tc * Tc;
     643             : 
     644        4394 :   partial_density = 1.0e6 * _Mco2 / V;
     645        4394 :   dpartial_density_dT = -1.0e6 * _Mco2 * dV_dT / V / V;
     646        6893 : }

Generated by: LCOV version 1.11