LCOV - code coverage report
Current view: top level - src/materials - PorousFlowFluidStateWaterNCG.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 168 172 97.7 %
Date: 2017-11-20 10:17:24 Functions: 6 7 85.7 %
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 "PorousFlowFluidStateWaterNCG.h"
       9             : #include "PorousFlowCapillaryPressure.h"
      10             : 
      11             : template <>
      12             : InputParameters
      13          16 : validParams<PorousFlowFluidStateWaterNCG>()
      14             : {
      15          16 :   InputParameters params = validParams<PorousFlowFluidStateFlashBase>();
      16          48 :   params.addRequiredParam<UserObjectName>("water_fp", "The name of the user object for water");
      17          48 :   params.addRequiredParam<UserObjectName>(
      18          16 :       "gas_fp", "The name of the user object for the non-condensable gas");
      19          32 :   params.addClassDescription("Fluid state class for water and non-condensable gas");
      20          16 :   return params;
      21             : }
      22             : 
      23          48 : PorousFlowFluidStateWaterNCG::PorousFlowFluidStateWaterNCG(const InputParameters & parameters)
      24             :   : PorousFlowFluidStateFlashBase(parameters),
      25             : 
      26          96 :     _water_fp(getUserObject<Water97FluidProperties>("water_fp")),
      27          96 :     _ncg_fp(getUserObject<SinglePhaseFluidPropertiesPT>("gas_fp")),
      28          48 :     _Mh2o(_water_fp.molarMass()),
      29          48 :     _Mncg(_ncg_fp.molarMass()),
      30          48 :     _water_triple_temperature(_water_fp.triplePointTemperature()),
      31         288 :     _water_critical_temperature(_water_fp.criticalTemperature())
      32             : {
      33             :   // Check that the correct FluidProperties UserObjects have been provided
      34          96 :   if (_water_fp.fluidName() != "water")
      35           0 :     mooseError("Only a valid water FluidProperties UserObject can be provided in water_fp");
      36          48 : }
      37             : 
      38             : void
      39      341147 : PorousFlowFluidStateWaterNCG::thermophysicalProperties()
      40             : {
      41             :   // The FluidProperty objects use temperature in K
      42      682294 :   Real Tk = _temperature[_qp] + _T_c2k;
      43             : 
      44             :   // Check whether the input temperature is within the region of validity of this equation
      45             :   // of state (T_triple <= T <= T_critical)
      46      341147 :   if (Tk < _water_triple_temperature || Tk > _water_critical_temperature)
      47             :     mooseError("PorousFlowFluidStateWaterNCG: Temperature is outside range 273.16 K <= T "
      48           0 :                "<= 647.096 K");
      49             : 
      50             :   // Equilibrium constants for each component (Henry's law for the NCG
      51             :   // component, and Raoult's law for water).
      52             :   // Note: these are in terms of mole fraction
      53             :   Real psat, dpsat_dT;
      54      341147 :   _water_fp.vaporPressure_dT(Tk, psat, dpsat_dT);
      55      682294 :   Real K0 = _ncg_fp.henryConstant(Tk) / _gas_porepressure[_qp];
      56      341147 :   Real K1 = psat / _gas_porepressure[_qp];
      57             : 
      58             :   // The mole fractions for the NCG component in the two component
      59             :   // case can be expressed in terms of the equilibrium constants only
      60      341147 :   Real x0 = (1.0 - K1) / (K0 - K1);
      61      341147 :   Real y0 = K0 * x0;
      62             : 
      63             :   // Convert mole fractions to mass fractions
      64      341147 :   Real X0 = moleFractionToMassFraction(x0);
      65      341147 :   Real Y0 = moleFractionToMassFraction(y0);
      66             : 
      67             :   // Determine which phases are present based on the value of z
      68             :   bool is_liquid = false;
      69             :   bool is_gas = false;
      70             :   bool is_twophase = false;
      71             : 
      72      682294 :   if ((*_z[0])[_qp] <= X0)
      73             :   {
      74             :     // In this case, there is not enough NCG to form a gas phase,
      75             :     // so only a liquid phase is present
      76             :     is_liquid = true;
      77             :   }
      78       25114 :   else if ((*_z[0])[_qp] > X0 && (*_z[0])[_qp] < Y0)
      79             :   {
      80             :     // Two phases are present
      81             :     is_twophase = true;
      82             :   }
      83             :   else // ((*_z[0])[_qp] >= Y0)
      84             :   {
      85             :     // In this case, there is not enough water to form a liquid
      86             :     // phase, so only a gas phase is present
      87             :     is_gas = true;
      88             :   }
      89             : 
      90             :   // Material must provide the following properties
      91             :   Real vapor_mass_fraction = 0.0;
      92             :   Real gas_saturation, liquid_saturation;
      93             :   Real gas_density, liquid_density;
      94             :   Real gas_viscosity, liquid_viscosity;
      95             : 
      96             :   // And the following derivatives
      97             :   Real dX0_dp = 0.0, dX0_dT = 0.0, dX0_dz = 0.0;
      98             :   Real dY0_dp = 0.0, dY0_dT = 0.0, dY0_dz = 0.0;
      99             :   Real dgas_density_dp = 0.0, dgas_density_dT = 0.0, dgas_density_dz = 0.0;
     100      341147 :   Real dliquid_density_dp = 0.0, dliquid_density_dT = 0.0, dliquid_density_dz = 0.0;
     101             :   Real dgas_viscosity_dp = 0.0, dgas_viscosity_dT = 0.0, dgas_viscosity_dz = 0.0;
     102      341147 :   Real dliquid_viscosity_dp = 0.0, dliquid_viscosity_dT = 0.0, dliquid_viscosity_dz = 0.0;
     103             : 
     104      341147 :   if (is_liquid)
     105             :   {
     106      316033 :     X0 = (*_z[0])[_qp];
     107             :     Y0 = 0.0;
     108             :     gas_saturation = 0.0;
     109             :     dX0_dz = 1.0;
     110             :     gas_density = 0.0;
     111             :     gas_viscosity = 1.0; // To guard against division by 0
     112             :   }
     113       25114 :   else if (is_twophase)
     114             :   {
     115             :     // Two phase are present. Set mass equilibrium constants used in the
     116             :     // calculation of vapor mass fraction
     117       23666 :     std::vector<Real> Ki(_num_components);
     118       23666 :     Ki[0] = Y0 / X0;
     119       23666 :     Ki[1] = (1.0 - Y0) / (1.0 - X0);
     120       23666 :     vapor_mass_fraction = vaporMassFraction(Ki);
     121             : 
     122             :     // Derivatives of mass fractions wrt PorousFlow variables
     123             :     Real Kh, dKh_dT;
     124       23666 :     _ncg_fp.henryConstant_dT(Tk, Kh, dKh_dT);
     125       47332 :     Real dK0_dp = -Kh / _gas_porepressure[_qp] / _gas_porepressure[_qp];
     126       23666 :     Real dK0_dT = dKh_dT / _gas_porepressure[_qp];
     127             : 
     128       23666 :     Real dK1_dp = -psat / _gas_porepressure[_qp] / _gas_porepressure[_qp];
     129       23666 :     Real dK1_dT = dpsat_dT / _gas_porepressure[_qp];
     130             : 
     131       23666 :     Real dx0_dp = ((K1 - 1.0) * dK0_dp + (1 - K0) * dK1_dp) / (K0 - K1) / (K0 - K1);
     132       23666 :     Real dy0_dp = x0 * dK0_dp + K0 * dx0_dp;
     133       23666 :     Real dx0_dT = ((K1 - 1.0) * dK0_dT + (1 - K0) * dK1_dT) / (K0 - K1) / (K0 - K1);
     134       23666 :     Real dy0_dT = x0 * dK0_dT + K0 * dx0_dT;
     135             : 
     136             :     Real dX0_dx0 =
     137       23666 :         _Mncg * _Mh2o / (x0 * _Mncg + (1.0 - x0) * _Mh2o) / (x0 * _Mncg + (1.0 - x0) * _Mh2o);
     138             :     Real dY0_dy0 =
     139       23666 :         _Mncg * _Mh2o / (y0 * _Mncg + (1.0 - y0) * _Mh2o) / (y0 * _Mncg + (1.0 - y0) * _Mh2o);
     140             : 
     141       23666 :     dX0_dp = dX0_dx0 * dx0_dp;
     142       23666 :     dX0_dT = dX0_dx0 * dx0_dT;
     143       23666 :     dY0_dp = dY0_dy0 * dy0_dp;
     144       23666 :     dY0_dT = dY0_dy0 * dy0_dT;
     145             :   }
     146             :   else // if (is_gas)
     147             :   {
     148             :     X0 = 0.0;
     149        1448 :     Y0 = (*_z[0])[_qp];
     150             :     vapor_mass_fraction = 1.0;
     151             :     gas_saturation = 1.0;
     152             :     dY0_dz = 1.0;
     153        1448 :     liquid_density = 0.0;
     154        1448 :     liquid_viscosity = 1.0; // To guard against division by 0
     155             :   }
     156             : 
     157             :   // Calculate the gas density and viscosity in the single phase gas or two
     158             :   // phase region
     159      341147 :   if (is_gas || is_twophase)
     160             :   {
     161             :     Real ncg_density, dncg_density_dp, dncg_density_dT;
     162             :     Real vapor_density, dvapor_density_dp, dvapor_density_dT;
     163             :     // NCG density calculated using partial pressure Y0 * gas_poreressure (Dalton's law)
     164       50228 :     _ncg_fp.rho_dpT(Y0 * _gas_porepressure[_qp], Tk, ncg_density, dncg_density_dp, dncg_density_dT);
     165             :     // Vapor density calculated using partial pressure X1 * psat (Raoult's law)
     166       25114 :     _water_fp.rho_dpT((1.0 - X0) * psat, Tk, vapor_density, dvapor_density_dp, dvapor_density_dT);
     167             : 
     168             :     // The derivatives wrt pressure above must be multiplied by the derivative of the pressure
     169             :     // variable using the chain rule
     170       25114 :     gas_density = ncg_density + vapor_density;
     171       75342 :     dgas_density_dp = (Y0 + dY0_dp * _gas_porepressure[_qp]) * dncg_density_dp -
     172       25114 :                       dX0_dp * psat * dvapor_density_dp;
     173       25114 :     dgas_density_dT = dncg_density_dT + dvapor_density_dT;
     174             : 
     175             :     Real ncg_viscosity, dncg_viscosity_drho, dncg_viscosity_dT;
     176             :     Real vapor_viscosity, dvapor_viscosity_drho, dvapor_viscosity_dT;
     177       25114 :     _ncg_fp.mu_drhoT_from_rho_T(
     178       25114 :         ncg_density, Tk, dncg_density_dT, ncg_viscosity, dncg_viscosity_drho, dncg_viscosity_dT);
     179       25114 :     _water_fp.mu_drhoT_from_rho_T(vapor_density,
     180             :                                   Tk,
     181             :                                   dvapor_density_dT,
     182             :                                   vapor_viscosity,
     183             :                                   dvapor_viscosity_drho,
     184       25114 :                                   dvapor_viscosity_dT);
     185             : 
     186             :     // Assume that the viscosity of the gas phase is a weighted sum of the
     187             :     // individual viscosities
     188       25114 :     gas_viscosity = Y0 * ncg_viscosity + (1.0 - Y0) * vapor_viscosity;
     189       25114 :     dgas_viscosity_dp =
     190       50228 :         dY0_dp * (ncg_viscosity - vapor_viscosity) +
     191       50228 :         Y0 * (Y0 + dY0_dp * _gas_porepressure[_qp]) * dncg_viscosity_drho * dncg_density_dp -
     192       25114 :         dX0_dp * psat * (1.0 - Y0) * dvapor_viscosity_drho * dvapor_density_dp;
     193       50228 :     dgas_viscosity_dT = dY0_dT * (ncg_viscosity - vapor_viscosity) + Y0 * dncg_viscosity_dT +
     194       25114 :                         (1.0 - Y0) * dvapor_density_dT;
     195             : 
     196             :     // Also calculate derivatives wrt z in the gas phase (these are 0 in the two phase region)
     197       25114 :     if (is_gas)
     198             :     {
     199        1448 :       dgas_density_dz =
     200        1448 :           dY0_dz * _gas_porepressure[_qp] * dncg_density_dp - dX0_dz * psat * dvapor_density_dp;
     201             : 
     202        1448 :       dgas_viscosity_dz =
     203        2896 :           dY0_dz * (ncg_viscosity - vapor_viscosity) +
     204        1448 :           Y0 * dncg_viscosity_drho * dncg_density_dp * dY0_dz * _gas_porepressure[_qp] -
     205        1448 :           dX0_dz * psat * (1.0 - Y0) * dvapor_viscosity_drho * dvapor_density_dp;
     206             :     }
     207             :   }
     208             : 
     209             :   // Calculate the saturation in the two phase case using the vapor mass fraction
     210      341147 :   if (is_twophase)
     211             :   {
     212             :     // Liquid density is approximated by the water density.
     213             :     // Use old value of gas saturation to estimate liquid saturation
     214             :     Real liqsat = 1.0;
     215       23666 :     if (!_is_initqp)
     216       47252 :       liqsat -= _saturation_old[_qp][_gas_phase_number];
     217       47332 :     liquid_density = _water_fp.rho(_gas_porepressure[_qp] + _pc_uo.capillaryPressure(liqsat), Tk);
     218             : 
     219             :     // The gas saturation in the two phase case
     220       47332 :     gas_saturation = vapor_mass_fraction * liquid_density /
     221       23666 :                      (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
     222             :   }
     223             : 
     224             :   // Calculate the saturations and pressures for each phase
     225      341147 :   liquid_saturation = 1.0 - gas_saturation;
     226      682294 :   Real liquid_porepressure = _gas_porepressure[_qp] - _pc_uo.capillaryPressure(liquid_saturation);
     227             : 
     228             :   // Calculate liquid density and viscosity if in the two phase or single phase
     229             :   // liquid region, assuming they are not affected by the presence of dissolved
     230             :   // NCG. Note: the (small) contribution due to derivative of capillary pressure
     231             :   //  wrt pressure (using the chain rule) is not implemented.
     232      341147 :   if (is_liquid || is_twophase)
     233             :   {
     234             :     Real dliquid_viscosity_drho;
     235      339699 :     _water_fp.rho_dpT(
     236      339699 :         liquid_porepressure, Tk, liquid_density, dliquid_density_dp, dliquid_density_dT);
     237      339696 :     _water_fp.mu_drhoT_from_rho_T(liquid_density,
     238             :                                   Tk,
     239             :                                   dliquid_density_dT,
     240             :                                   liquid_viscosity,
     241             :                                   dliquid_viscosity_drho,
     242      339696 :                                   dliquid_viscosity_dT);
     243             : 
     244             :     // The derivative of viscosity wrt pressure is given by the chain rule
     245      339696 :     dliquid_viscosity_dp = dliquid_viscosity_drho * dliquid_density_dp;
     246             :   }
     247             : 
     248             :   // Save properties in FluidStateProperties vector
     249      341144 :   auto & liquid = _fsp[_aqueous_phase_number];
     250      341144 :   auto & gas = _fsp[_gas_phase_number];
     251             : 
     252      341144 :   liquid.saturation = liquid_saturation;
     253      341144 :   gas.saturation = gas_saturation;
     254      341144 :   liquid.pressure = liquid_porepressure;
     255      682288 :   gas.pressure = _gas_porepressure[_qp];
     256             : 
     257      682288 :   liquid.mass_fraction[_aqueous_fluid_component] = 1.0 - X0;
     258      682288 :   liquid.mass_fraction[_gas_fluid_component] = X0;
     259      341144 :   gas.mass_fraction[_aqueous_fluid_component] = 1.0 - Y0;
     260      341144 :   gas.mass_fraction[_gas_fluid_component] = Y0;
     261             : 
     262      341144 :   liquid.fluid_density = liquid_density;
     263      341144 :   gas.fluid_density = gas_density;
     264      341144 :   liquid.fluid_viscosity = liquid_viscosity;
     265      341144 :   gas.fluid_viscosity = gas_viscosity;
     266             : 
     267             :   // Derivatives wrt PorousFlow variables are required by the kernels
     268             :   // Note: these don't need to be stateful so don't calculate them in
     269             :   // initQpStatefulProperties
     270      341144 :   if (!_is_initqp)
     271             :   {
     272             :     // Derivative of gas saturation wrt variables
     273             :     Real ds_dp = 0.0, ds_dT = 0.0, ds_dz = 0.0;
     274      340224 :     if (is_twophase)
     275             :     {
     276       23626 :       K0 = Y0 / X0;
     277       23626 :       K1 = (1.0 - Y0) / (1.0 - X0);
     278       23626 :       Real dv_dz = (K1 - K0) / ((K0 - 1.0) * (K1 - 1.0));
     279       47252 :       ds_dz = gas_density * liquid_density * dv_dz +
     280       47252 :               vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
     281       23626 :                   (gas_density * dliquid_density_dz - dgas_density_dz * liquid_density);
     282       23626 :       ds_dz /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
     283             :                (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
     284             : 
     285             :       Real Kh, dKh_dT;
     286       23626 :       _ncg_fp.henryConstant_dT(Tk, Kh, dKh_dT);
     287       47252 :       Real dK0_dp = -Kh / _gas_porepressure[_qp] / _gas_porepressure[_qp];
     288       23626 :       Real dK0_dT = dKh_dT / _gas_porepressure[_qp];
     289             : 
     290       23626 :       Real dK1_dp = -psat / _gas_porepressure[_qp] / _gas_porepressure[_qp];
     291       23626 :       Real dK1_dT = dpsat_dT / _gas_porepressure[_qp];
     292             : 
     293       47252 :       Real dv_dp = (*_z[0])[_qp] * dK1_dp / (K1 - 1.0) / (K1 - 1.0) +
     294       23626 :                    (1.0 - (*_z[0])[_qp]) * dK0_dp / (K0 - 1.0) / (K0 - 1.0);
     295       47252 :       ds_dp = gas_density * liquid_density * dv_dp +
     296       23626 :               vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
     297       23626 :                   (gas_density * dliquid_density_dp - dgas_density_dp * liquid_density);
     298             : 
     299       23626 :       ds_dp /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
     300             :                (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
     301             : 
     302       23626 :       Real dv_dT = (*_z[0])[_qp] * dK1_dT / (K1 - 1.0) / (K1 - 1.0) +
     303       23626 :                    (1.0 - (*_z[0])[_qp]) * dK0_dT / (K0 - 1.0) / (K0 - 1.0);
     304       47252 :       ds_dT = gas_density * liquid_density * dv_dT +
     305       23626 :               vapor_mass_fraction * (1.0 - vapor_mass_fraction) *
     306       23626 :                   (gas_density * dliquid_density_dT - dgas_density_dT * liquid_density);
     307             : 
     308       23626 :       ds_dT /= (gas_density + vapor_mass_fraction * (liquid_density - gas_density)) *
     309             :                (gas_density + vapor_mass_fraction * (liquid_density - gas_density));
     310             :     }
     311             : 
     312      340224 :     liquid.dsaturation_dp = -ds_dp;
     313      340224 :     liquid.dsaturation_dT = -ds_dT;
     314      340224 :     liquid.dsaturation_dz = -ds_dz;
     315      340224 :     gas.dsaturation_dp = ds_dp;
     316      340224 :     gas.dsaturation_dT = ds_dT;
     317      340224 :     gas.dsaturation_dz = ds_dz;
     318             : 
     319      680448 :     liquid.dmass_fraction_dp[_aqueous_fluid_component] = -dX0_dp;
     320      680448 :     liquid.dmass_fraction_dp[_gas_fluid_component] = dX0_dp;
     321      340224 :     liquid.dmass_fraction_dT[_aqueous_fluid_component] = -dX0_dT;
     322      340224 :     liquid.dmass_fraction_dT[_gas_fluid_component] = dX0_dT;
     323      340224 :     liquid.dmass_fraction_dz[_aqueous_fluid_component] = -dX0_dz;
     324      340224 :     liquid.dmass_fraction_dz[_gas_fluid_component] = dX0_dz;
     325             : 
     326      340224 :     gas.dmass_fraction_dp[_aqueous_fluid_component] = -dY0_dp;
     327      340224 :     gas.dmass_fraction_dp[_gas_fluid_component] = dY0_dp;
     328      340224 :     gas.dmass_fraction_dT[_aqueous_fluid_component] = -dY0_dT;
     329      340224 :     gas.dmass_fraction_dT[_gas_fluid_component] = dY0_dT;
     330      340224 :     gas.dmass_fraction_dz[_aqueous_fluid_component] = -dY0_dz;
     331      340224 :     gas.dmass_fraction_dz[_gas_fluid_component] = dY0_dz;
     332             : 
     333      340224 :     liquid.dfluid_density_dp = dliquid_density_dp;
     334      340224 :     liquid.dfluid_density_dT = dliquid_density_dT;
     335      340224 :     liquid.dfluid_density_dz = dliquid_density_dz;
     336      340224 :     gas.dfluid_density_dp = dgas_density_dp;
     337      340224 :     gas.dfluid_density_dT = dgas_density_dT;
     338      340224 :     gas.dfluid_density_dz = dgas_density_dz;
     339             : 
     340      340224 :     liquid.dfluid_viscosity_dp = dliquid_viscosity_dp;
     341      340224 :     liquid.dfluid_viscosity_dT = dliquid_viscosity_dT;
     342      340224 :     liquid.dfluid_viscosity_dz = dliquid_viscosity_dz;
     343      340224 :     gas.dfluid_viscosity_dp = dgas_viscosity_dp;
     344      340224 :     gas.dfluid_viscosity_dT = dgas_viscosity_dT;
     345      340224 :     gas.dfluid_viscosity_dz = dgas_viscosity_dz;
     346             :   }
     347      341144 : }
     348             : 
     349             : Real
     350           0 : PorousFlowFluidStateWaterNCG::enthalpyOfDissolution(Real temperature, Real Kh, Real dKh_dT) const
     351             : {
     352           0 :   return -_R * temperature * temperature * _Mncg * dKh_dT / Kh;
     353             : }
     354             : 
     355             : Real
     356      682294 : PorousFlowFluidStateWaterNCG::moleFractionToMassFraction(Real xmol) const
     357             : {
     358      682294 :   return xmol * _Mncg / (xmol * _Mncg + (1.0 - xmol) * _Mh2o);
     359        2499 : }

Generated by: LCOV version 1.11