LCOV - code coverage report
Current view: top level - src/materials - PorousFlowFluidStateFlashBase.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 182 208 87.5 %
Date: 2017-11-21 14:47:27 Functions: 8 10 80.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 "PorousFlowFluidStateFlashBase.h"
       9             : #include "PorousFlowCapillaryPressure.h"
      10             : 
      11             : template <>
      12             : InputParameters
      13          24 : validParams<PorousFlowFluidStateFlashBase>()
      14             : {
      15          24 :   InputParameters params = validParams<PorousFlowVariableBase>();
      16          72 :   params.addRequiredCoupledVar("gas_porepressure",
      17          24 :                                "Variable that is the porepressure of the gas phase");
      18          72 :   params.addRequiredCoupledVar("z", "Total mass fraction of component i summed over all phases");
      19          72 :   params.addParam<unsigned int>("liquid_phase_number", 0, "The phase number of the liquid phase");
      20          72 :   params.addParam<unsigned int>(
      21          24 :       "liquid_fluid_component", 0, "The fluid component number of the liquid phase");
      22          72 :   MooseEnum unit_choice("Kelvin=0 Celsius=1", "Kelvin");
      23          72 :   params.addParam<MooseEnum>(
      24          24 :       "temperature_unit", unit_choice, "The unit of the temperature variable");
      25          72 :   params.addRequiredParam<UserObjectName>("capillary_pressure",
      26          24 :                                           "Name of the UserObject defining the capillary pressure");
      27          48 :   params.addClassDescription("Base class for fluid state calculations using persistent primary "
      28          24 :                              "variables and a vapor-liquid flash");
      29          24 :   return params;
      30             : }
      31             : 
      32          72 : PorousFlowFluidStateFlashBase::PorousFlowFluidStateFlashBase(const InputParameters & parameters)
      33             :   : PorousFlowVariableBase(parameters),
      34             : 
      35         216 :     _gas_porepressure(_nodal_material ? coupledNodalValue("gas_porepressure")
      36         144 :                                       : coupledValue("gas_porepressure")),
      37         144 :     _gas_gradp_qp(coupledGradient("gas_porepressure")),
      38         144 :     _gas_porepressure_varnum(coupled("gas_porepressure")),
      39          72 :     _pvar(_dictator.isPorousFlowVariable(_gas_porepressure_varnum)
      40          72 :               ? _dictator.porousFlowVariableNum(_gas_porepressure_varnum)
      41             :               : 0),
      42             : 
      43         144 :     _num_z_vars(coupledComponents("z")),
      44             : 
      45         216 :     _aqueous_phase_number(getParam<unsigned int>("liquid_phase_number")),
      46          72 :     _gas_phase_number(1 - _aqueous_phase_number),
      47         216 :     _aqueous_fluid_component(getParam<unsigned int>("liquid_fluid_component")),
      48          72 :     _gas_fluid_component(1 - _aqueous_fluid_component),
      49             : 
      50         216 :     _temperature(_nodal_material ? getMaterialProperty<Real>("PorousFlow_temperature_nodal")
      51         144 :                                  : getMaterialProperty<Real>("PorousFlow_temperature_qp")),
      52         144 :     _gradT_qp(getMaterialProperty<RealGradient>("PorousFlow_grad_temperature_qp")),
      53             :     _dtemperature_dvar(
      54          72 :         _nodal_material
      55         216 :             ? getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")
      56         144 :             : getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
      57             : 
      58          72 :     _mass_frac(_nodal_material
      59         216 :                    ? declareProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")
      60         144 :                    : declareProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
      61         108 :     _grad_mass_frac_qp(_nodal_material ? nullptr
      62             :                                        : &declareProperty<std::vector<std::vector<RealGradient>>>(
      63         144 :                                              "PorousFlow_grad_mass_frac_qp")),
      64         144 :     _dmass_frac_dvar(_nodal_material ? declareProperty<std::vector<std::vector<std::vector<Real>>>>(
      65         144 :                                            "dPorousFlow_mass_frac_nodal_dvar")
      66             :                                      : declareProperty<std::vector<std::vector<std::vector<Real>>>>(
      67         144 :                                            "dPorousFlow_mass_frac_qp_dvar")),
      68          72 :     _saturation_old(_nodal_material
      69         216 :                         ? getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")
      70         144 :                         : getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_qp")),
      71             : 
      72          72 :     _fluid_density(_nodal_material
      73         216 :                        ? declareProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")
      74         144 :                        : declareProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      75         144 :     _dfluid_density_dvar(_nodal_material ? declareProperty<std::vector<std::vector<Real>>>(
      76         144 :                                                "dPorousFlow_fluid_phase_density_nodal_dvar")
      77             :                                          : declareProperty<std::vector<std::vector<Real>>>(
      78         144 :                                                "dPorousFlow_fluid_phase_density_qp_dvar")),
      79          72 :     _fluid_viscosity(_nodal_material
      80         216 :                          ? declareProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
      81         144 :                          : declareProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      82             :     _dfluid_viscosity_dvar(
      83          72 :         _nodal_material
      84         216 :             ? declareProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")
      85         144 :             : declareProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_qp_dvar")),
      86             : 
      87         216 :     _T_c2k(getParam<MooseEnum>("temperature_unit") == 0 ? 0.0 : 273.15),
      88             :     _R(8.3144598),
      89             :     _nr_max_its(42),
      90             :     _nr_tol(1.0e-12),
      91             :     _is_initqp(false),
      92        2016 :     _pc_uo(getUserObject<PorousFlowCapillaryPressure>("capillary_pressure"))
      93             : {
      94             :   // Check that the number of total mass fractions provided as primary variables is correct
      95          72 :   if (_num_z_vars != _num_components - 1)
      96           0 :     mooseError("The number of supplied mass fraction variables should be ",
      97           0 :                _num_components - 1,
      98             :                " in ",
      99           0 :                _name,
     100             :                " but ",
     101             :                _num_z_vars,
     102           0 :                " are supplied");
     103             : 
     104             :   // Store all total mass fractions and associated variable numbers
     105          72 :   _z.resize(_num_z_vars);
     106          72 :   _gradz_qp.resize(_num_z_vars);
     107          72 :   _z_varnum.resize(_num_z_vars);
     108          72 :   _zvar.resize(_num_z_vars);
     109             : 
     110         216 :   for (unsigned int i = 0; i < _num_z_vars; ++i)
     111             :   {
     112         216 :     _z[i] = (_nodal_material ? &coupledNodalValue("z", i) : &coupledValue("z", i));
     113         144 :     _gradz_qp[i] = &coupledGradient("z", i);
     114         144 :     _z_varnum[i] = coupled("z", i);
     115         144 :     _zvar[i] = (_dictator.isPorousFlowVariable(_z_varnum[i])
     116         144 :                     ? _dictator.porousFlowVariableNum(_z_varnum[i])
     117             :                     : 0);
     118             :   }
     119             : 
     120             :   // Set the size of the FluidStateProperties vector
     121          72 :   _fsp.resize(_num_phases);
     122             : 
     123             :   // Set the size of the mass fraction vectors for each phase
     124         360 :   for (unsigned int ph = 0; ph < _num_phases; ++ph)
     125             :   {
     126         288 :     _fsp[ph].mass_fraction.resize(_num_components);
     127         288 :     _fsp[ph].dmass_fraction_dp.resize(_num_components);
     128         288 :     _fsp[ph].dmass_fraction_dT.resize(_num_components);
     129         288 :     _fsp[ph].dmass_fraction_dz.resize(_num_components);
     130             :   }
     131          72 : }
     132             : 
     133             : void
     134        2048 : PorousFlowFluidStateFlashBase::initQpStatefulProperties()
     135             : {
     136        2048 :   _is_initqp = true;
     137             :   // Set the size of pressure and saturation vectors
     138        2048 :   PorousFlowVariableBase::initQpStatefulProperties();
     139             : 
     140             :   // Set the size of all other vectors
     141        2048 :   setMaterialVectorSize();
     142             : 
     143             :   // Set the initial values of the properties at the nodes.
     144             :   // Note: not required for qp materials as no old values at the qps are requested
     145        2048 :   if (_nodal_material)
     146             :   {
     147        1024 :     thermophysicalProperties();
     148             : 
     149        5120 :     for (unsigned int ph = 0; ph < _num_phases; ++ph)
     150             :     {
     151        4096 :       _saturation[_qp][ph] = _fsp[ph].saturation;
     152        4096 :       _porepressure[_qp][ph] = _fsp[ph].pressure;
     153        4096 :       _fluid_density[_qp][ph] = _fsp[ph].fluid_density;
     154        4096 :       _fluid_viscosity[_qp][ph] = _fsp[ph].fluid_viscosity;
     155        4096 :       _mass_frac[_qp][ph] = _fsp[ph].mass_fraction;
     156             :     }
     157             :   }
     158        2048 : }
     159             : 
     160             : void
     161      344390 : PorousFlowFluidStateFlashBase::computeQpProperties()
     162             : {
     163      344390 :   _is_initqp = false;
     164             :   // Prepare the derivative vectors
     165      344390 :   PorousFlowVariableBase::computeQpProperties();
     166             : 
     167             :   // Set the size of all other vectors
     168      344390 :   setMaterialVectorSize();
     169             : 
     170             :   // Calculate all required thermophysical properties
     171      344390 :   thermophysicalProperties();
     172             : 
     173     1721920 :   for (unsigned int ph = 0; ph < _num_phases; ++ph)
     174             :   {
     175     1377536 :     _saturation[_qp][ph] = _fsp[ph].saturation;
     176     1377536 :     _porepressure[_qp][ph] = _fsp[ph].pressure;
     177     1377536 :     _fluid_density[_qp][ph] = _fsp[ph].fluid_density;
     178     1377536 :     _fluid_viscosity[_qp][ph] = _fsp[ph].fluid_viscosity;
     179     1377536 :     _mass_frac[_qp][ph] = _fsp[ph].mass_fraction;
     180             :   }
     181             : 
     182             :   // Derivative of saturation wrt variables
     183     1721920 :   for (unsigned int ph = 0; ph < _num_phases; ++ph)
     184             :   {
     185     1377536 :     _dsaturation_dvar[_qp][ph][_zvar[0]] = _fsp[ph].dsaturation_dz;
     186     1377536 :     _dsaturation_dvar[_qp][ph][_pvar] = _fsp[ph].dsaturation_dp;
     187             :   }
     188             :   // Derivative of capillary pressure
     189      688768 :   Real dpc = _pc_uo.dCapillaryPressure(_fsp[_aqueous_phase_number].saturation);
     190             : 
     191             :   // Derivative of porepressure wrt variables
     192      344384 :   if (_dictator.isPorousFlowVariable(_gas_porepressure_varnum))
     193             :   {
     194     1721920 :     for (unsigned int ph = 0; ph < _num_phases; ++ph)
     195             :     {
     196     1377536 :       _dporepressure_dvar[_qp][ph][_pvar] = 1.0;
     197      688768 :       if (!_nodal_material)
     198      688768 :         (*_dgradp_qp_dgradv)[_qp][ph][_pvar] = 1.0;
     199             :     }
     200             : 
     201      344384 :     if (!_nodal_material)
     202      344384 :       (*_dgradp_qp_dgradv)[_qp][_aqueous_phase_number][_zvar[0]] =
     203      516576 :           -dpc * _fsp[_aqueous_phase_number].dsaturation_dp -
     204      172192 :           dpc * _fsp[_aqueous_phase_number].dsaturation_dz;
     205             : 
     206             :     // The aqueous phase porepressure is also a function of liquid saturation,
     207             :     // which depends on both gas porepressure and z
     208      688768 :     _dporepressure_dvar[_qp][_aqueous_phase_number][_pvar] +=
     209     1033152 :         -dpc * _dsaturation_dvar[_qp][_aqueous_phase_number][_pvar];
     210      688768 :     _dporepressure_dvar[_qp][_aqueous_phase_number][_zvar[0]] =
     211      344384 :         -dpc * _dsaturation_dvar[_qp][_aqueous_phase_number][_zvar[0]];
     212             :   }
     213             : 
     214             :   // Calculate derivatives of material properties wrt primary variables
     215             :   // Derivative of z wrt variables
     216             :   std::vector<Real> dz_dvar;
     217      688768 :   dz_dvar.assign(_num_pf_vars, 0.0);
     218      344384 :   if (_dictator.isPorousFlowVariable(_z_varnum[0]))
     219      688768 :     dz_dvar[_zvar[0]] = 1.0;
     220             : 
     221             :   // Derivatives of properties wrt primary variables
     222     1721920 :   for (unsigned int v = 0; v < _num_pf_vars; ++v)
     223             :   {
     224     3443840 :     for (unsigned int ph = 0; ph < _num_phases; ++ph)
     225             :     {
     226             :       // Derivative of density in each phase
     227     2755072 :       _dfluid_density_dvar[_qp][ph][v] =
     228     2755072 :           _fsp[ph].dfluid_density_dp * _dporepressure_dvar[_qp][ph][v];
     229     2755072 :       _dfluid_density_dvar[_qp][ph][v] += _fsp[ph].dfluid_density_dT * _dtemperature_dvar[_qp][v];
     230     2755072 :       _dfluid_density_dvar[_qp][ph][v] += _fsp[ph].dfluid_density_dz * dz_dvar[v];
     231             : 
     232             :       // Derivative of viscosity in each phase
     233     2755072 :       _dfluid_viscosity_dvar[_qp][ph][v] =
     234     1377536 :           _fsp[ph].dfluid_viscosity_dp * _dporepressure_dvar[_qp][ph][v];
     235     1377536 :       _dfluid_viscosity_dvar[_qp][ph][v] +=
     236     1377536 :           _fsp[ph].dfluid_viscosity_dT * _dtemperature_dvar[_qp][v];
     237     1377536 :       _dfluid_viscosity_dvar[_qp][ph][v] += _fsp[ph].dfluid_viscosity_dz * dz_dvar[v];
     238             : 
     239             :       // The derivative of the mass fractions for each fluid component in each phase
     240     6887680 :       for (unsigned int comp = 0; comp < _num_components; ++comp)
     241             :       {
     242     5510144 :         _dmass_frac_dvar[_qp][ph][comp][v] =
     243     2755072 :             _fsp[ph].dmass_fraction_dp[comp] * _dporepressure_dvar[_qp][ph][v];
     244     2755072 :         _dmass_frac_dvar[_qp][ph][comp][v] +=
     245     2755072 :             _fsp[ph].dmass_fraction_dT[comp] * _dtemperature_dvar[_qp][v];
     246     2755072 :         _dmass_frac_dvar[_qp][ph][comp][v] += _fsp[ph].dmass_fraction_dz[comp] * dz_dvar[v];
     247             :       }
     248             :     }
     249             :   }
     250             : 
     251             :   // If the material properties are being evaluated at the qps, calculate the
     252             :   // gradients as well. Note: only nodal properties are evaluated in
     253             :   // initQpStatefulProperties(), so no need to check _is_initqp flag for qp
     254             :   // properties
     255      344384 :   if (!_nodal_material)
     256             :   {
     257             :     // Second derivative of capillary pressure
     258      344384 :     Real d2pc = _pc_uo.d2CapillaryPressure(_fsp[_aqueous_phase_number].saturation);
     259             : 
     260      344384 :     (*_grads_qp)[_qp][_gas_phase_number] =
     261      516576 :         _dsaturation_dvar[_qp][_gas_phase_number][_pvar] * _gas_gradp_qp[_qp] +
     262      516576 :         _dsaturation_dvar[_qp][_gas_phase_number][_zvar[0]] * (*_gradz_qp[0])[_qp];
     263      344384 :     (*_grads_qp)[_qp][_aqueous_phase_number] = -(*_grads_qp)[_qp][_gas_phase_number];
     264             : 
     265      344384 :     (*_gradp_qp)[_qp][_gas_phase_number] = _gas_gradp_qp[_qp];
     266      172192 :     (*_gradp_qp)[_qp][_aqueous_phase_number] =
     267             :         _gas_gradp_qp[_qp] - dpc * (*_grads_qp)[_qp][_aqueous_phase_number];
     268             : 
     269      344384 :     (*_dgradp_qp_dv)[_qp][_aqueous_phase_number][_zvar[0]] =
     270      172192 :         -d2pc * (*_grads_qp)[_qp][_aqueous_phase_number];
     271             : 
     272      344384 :     (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_aqueous_fluid_component] =
     273      172192 :         _fsp[_aqueous_phase_number].dmass_fraction_dp[_aqueous_fluid_component] *
     274             :             _gas_gradp_qp[_qp] +
     275      172192 :         _fsp[_aqueous_phase_number].dmass_fraction_dz[_aqueous_fluid_component] *
     276             :             (*_gradz_qp[0])[_qp];
     277      344384 :     (*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_gas_fluid_component] =
     278             :         -(*_grad_mass_frac_qp)[_qp][_aqueous_phase_number][_aqueous_fluid_component];
     279      172192 :     (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_aqueous_fluid_component] =
     280      172192 :         _fsp[_gas_phase_number].dmass_fraction_dp[_aqueous_fluid_component] * _gas_gradp_qp[_qp] +
     281      172192 :         _fsp[_gas_phase_number].dmass_fraction_dz[_aqueous_fluid_component] * (*_gradz_qp[0])[_qp];
     282      172192 :     (*_grad_mass_frac_qp)[_qp][_gas_phase_number][_gas_fluid_component] =
     283             :         -(*_grad_mass_frac_qp)[_qp][_gas_phase_number][_aqueous_fluid_component];
     284             :   }
     285      344384 : }
     286             : 
     287             : void
     288      346438 : PorousFlowFluidStateFlashBase::setMaterialVectorSize() const
     289             : {
     290      692876 :   _fluid_density[_qp].assign(_num_phases, 0.0);
     291      692876 :   _fluid_viscosity[_qp].assign(_num_phases, 0.0);
     292      692876 :   _mass_frac[_qp].resize(_num_phases);
     293             : 
     294             :   // Derivatives and gradients are not required in initQpStatefulProperties
     295      346438 :   if (!_is_initqp)
     296             :   {
     297      688780 :     _dfluid_density_dvar[_qp].resize(_num_phases);
     298      688780 :     _dfluid_viscosity_dvar[_qp].resize(_num_phases);
     299      688780 :     _dmass_frac_dvar[_qp].resize(_num_phases);
     300             : 
     301      344390 :     if (!_nodal_material)
     302      344396 :       (*_grad_mass_frac_qp)[_qp].resize(_num_phases);
     303             : 
     304     1721950 :     for (unsigned int ph = 0; ph < _num_phases; ++ph)
     305             :     {
     306     1377560 :       _dfluid_density_dvar[_qp][ph].assign(_num_pf_vars, 0.0);
     307     1377560 :       _dfluid_viscosity_dvar[_qp][ph].assign(_num_pf_vars, 0.0);
     308     1377560 :       _dmass_frac_dvar[_qp][ph].resize(_num_components);
     309             : 
     310     3443900 :       for (unsigned int comp = 0; comp < _num_components; ++comp)
     311     2755120 :         _dmass_frac_dvar[_qp][ph][comp].assign(_num_pf_vars, 0.0);
     312             : 
     313      688780 :       if (!_nodal_material)
     314      688792 :         (*_grad_mass_frac_qp)[_qp][ph].assign(_num_components, RealGradient());
     315             :     }
     316             :   }
     317      346438 : }
     318             : 
     319             : Real
     320           0 : PorousFlowFluidStateFlashBase::rachfordRice(Real x, std::vector<Real> & Ki) const
     321             : {
     322             :   // Check that the size of the equilibrium constant vector is correct
     323           0 :   if (Ki.size() != _num_components)
     324           0 :     mooseError("The number of equilibrium components passed to rachfordRice is not correct");
     325             : 
     326             :   Real f = 0.0;
     327             :   Real z_total = 0.0;
     328             : 
     329           0 :   for (unsigned int i = 0; i < _num_z_vars; ++i)
     330             :   {
     331           0 :     f += (*_z[i])[_qp] * (Ki[i] - 1.0) / (1.0 + x * (Ki[i] - 1.0));
     332           0 :     z_total += (*_z[i])[_qp];
     333             :   }
     334             : 
     335             :   // Add the last component (with total mass fraction = 1 - z_total)
     336           0 :   f += (1.0 - z_total) * (Ki[_num_z_vars] - 1.0) / (1.0 + x * (Ki[_num_z_vars] - 1.0));
     337             : 
     338           0 :   return f;
     339             : }
     340             : 
     341             : Real
     342           0 : PorousFlowFluidStateFlashBase::rachfordRiceDeriv(Real x, std::vector<Real> & Ki) const
     343             : {
     344             :   // Check that the size of the equilibrium constant vector is correct
     345           0 :   if (Ki.size() != _num_components)
     346           0 :     mooseError("The number of equilibrium components passed to rachfordRiceDeriv is not correct");
     347             : 
     348             :   Real df = 0.0;
     349             :   Real z_total = 0.0;
     350             : 
     351           0 :   for (unsigned int i = 0; i < _num_z_vars; ++i)
     352             :   {
     353           0 :     df -= (*_z[i])[_qp] * (Ki[i] - 1.0) * (Ki[i] - 1.0) / (1.0 + x * (Ki[i] - 1.0)) /
     354             :           (1.0 + x * (Ki[i] - 1.0));
     355           0 :     z_total += (*_z[i])[_qp];
     356             :   }
     357             : 
     358             :   // Add the last component (with total mass fraction = 1 - z_total)
     359           0 :   df -= (1.0 - z_total) * (Ki[_num_z_vars] - 1.0) * (Ki[_num_z_vars] - 1.0) /
     360           0 :         (1.0 + x * (Ki[_num_z_vars] - 1.0)) / (1.0 + x * (Ki[_num_z_vars] - 1.0));
     361             : 
     362           0 :   return df;
     363             : }
     364             : 
     365             : Real
     366       25204 : PorousFlowFluidStateFlashBase::vaporMassFraction(std::vector<Real> & Ki) const
     367             : {
     368             :   // Check that the size of the equilibrium constant vector is correct
     369       25204 :   if (Ki.size() != _num_components)
     370           0 :     mooseError("The number of equilibrium components passed to vaporMassFraction is not correct");
     371             : 
     372             :   Real v;
     373             : 
     374             :   // If there are only two components, an analytical solution is possible
     375       25204 :   if (_num_components == 2)
     376       50408 :     v = ((*_z[0])[_qp] * (Ki[1] - Ki[0]) - (Ki[1] - 1.0)) / ((Ki[0] - 1.0) * (Ki[1] - 1.0));
     377             :   else
     378             :   {
     379             :     // More than two components - solve the Rachford-Rice equation using
     380             :     // Newton-Raphson method.
     381             :     // Initial guess for vapor mass fraction
     382             :     Real v0 = 0.5;
     383             :     unsigned int iter = 0;
     384             : 
     385           0 :     while (std::abs(rachfordRice(v0, Ki)) > _nr_tol)
     386             :     {
     387           0 :       v0 = v0 - rachfordRice(v0, Ki) / rachfordRiceDeriv(v0, Ki);
     388           0 :       iter++;
     389             : 
     390           0 :       if (iter > _nr_max_its)
     391             :         break;
     392             :     }
     393             :     v = v0;
     394             :   }
     395       25204 :   return v;
     396        2499 : }

Generated by: LCOV version 1.11