Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "PorousFlowFluidStateSingleComponent.h"
      11             : #include "PorousFlowCapillaryPressure.h"
      12             : 
      13             : registerMooseObject("PorousFlowApp", PorousFlowFluidStateSingleComponent);
      14             : registerMooseObject("PorousFlowApp", ADPorousFlowFluidStateSingleComponent);
      15             : 
      16             : template <bool is_ad>
      17             : InputParameters
      18        1190 : PorousFlowFluidStateSingleComponentTempl<is_ad>::validParams()
      19             : {
      20        1190 :   InputParameters params = PorousFlowFluidStateBaseMaterialTempl<is_ad>::validParams();
      21        2380 :   params.addRequiredCoupledVar("porepressure",
      22             :                                "Variable that is the porepressure of the liquid phase");
      23        2380 :   params.addRequiredCoupledVar("enthalpy", "Enthalpy of the fluid");
      24        2380 :   params.addRequiredParam<UserObjectName>("fluid_state", "Name of the FluidState UserObject");
      25        1190 :   params.addClassDescription(
      26             :       "Class for single component multiphase fluid state calculations using pressure and enthalpy");
      27        1190 :   return params;
      28           0 : }
      29             : 
      30             : template <bool is_ad>
      31         924 : PorousFlowFluidStateSingleComponentTempl<is_ad>::PorousFlowFluidStateSingleComponentTempl(
      32             :     const InputParameters & parameters)
      33             :   : PorousFlowFluidStateBaseMaterialTempl<is_ad>(parameters),
      34        1848 :     _liquid_porepressure(_nodal_material
      35         462 :                              ? this->template coupledGenericDofValue<is_ad>("porepressure")
      36        1386 :                              : this->template coupledGenericValue<is_ad>("porepressure")),
      37         924 :     _liquid_gradp_qp(this->template coupledGenericGradient<is_ad>("porepressure")),
      38         924 :     _liquid_porepressure_varnum(coupled("porepressure")),
      39         924 :     _pvar(_dictator.isPorousFlowVariable(_liquid_porepressure_varnum)
      40         924 :               ? _dictator.porousFlowVariableNum(_liquid_porepressure_varnum)
      41             :               : 0),
      42         924 :     _enthalpy(_nodal_material ? this->template coupledGenericDofValue<is_ad>("enthalpy")
      43        1386 :                               : this->template coupledGenericValue<is_ad>("enthalpy")),
      44         924 :     _gradh_qp(this->template coupledGenericGradient<is_ad>("enthalpy")),
      45         924 :     _enthalpy_varnum(coupled("enthalpy")),
      46         924 :     _hvar(_dictator.isPorousFlowVariable(_enthalpy_varnum)
      47         924 :               ? _dictator.porousFlowVariableNum(_enthalpy_varnum)
      48             :               : 0),
      49         924 :     _fs(this->template getUserObject<PorousFlowFluidStateSingleComponentBase>("fluid_state")),
      50         924 :     _aqueous_phase_number(_fs.aqueousPhaseIndex()),
      51         924 :     _gas_phase_number(_fs.gasPhaseIndex()),
      52         924 :     _temperature(
      53         924 :         this->template declareGenericProperty<Real, is_ad>("PorousFlow_temperature" + _sfx)),
      54        1848 :     _grad_temperature_qp(_nodal_material
      55         924 :                              ? nullptr
      56        1386 :                              : &this->template declareGenericProperty<RealGradient, is_ad>(
      57             :                                    "PorousFlow_grad_temperature_qp")),
      58         924 :     _dtemperature_dvar(is_ad ? nullptr
      59        1848 :                              : &this->template declareProperty<std::vector<Real>>(
      60             :                                    "dPorousFlow_temperature" + _sfx + "_dvar")),
      61        1848 :     _dgrad_temperature_dgradv(is_ad || _nodal_material
      62         924 :                                   ? nullptr
      63         924 :                                   : &this->template declareProperty<std::vector<Real>>(
      64             :                                         "dPorousFlow_grad_temperature_qp_dgradvar")),
      65        1848 :     _dgrad_temperature_dv(is_ad ? nullptr
      66         924 :                           : _nodal_material
      67             :                               ? nullptr
      68         924 :                               : &this->template declareProperty<std::vector<RealGradient>>(
      69             :                                     "dPorousFlow_grad_temperature_qp_dvar")),
      70         924 :     _pidx(_fs.getPressureIndex()),
      71        1848 :     _hidx(_fs.getEnthalpyIndex())
      72             : {
      73             :   // Check that the number of phases in the fluidstate class is also provided in the Dictator
      74         924 :   if (_fs.numPhases() != _num_phases)
      75           0 :     mooseError(name(),
      76             :                ": only ",
      77           0 :                _fs.numPhases(),
      78             :                " phases are allowed. Please check the number of phases entered in the dictator is "
      79             :                "correct");
      80         924 : }
      81             : 
      82             : template <bool is_ad>
      83             : void
      84      134300 : PorousFlowFluidStateSingleComponentTempl<is_ad>::thermophysicalProperties()
      85             : {
      86      134300 :   _fs.clearFluidStateProperties(_fsp);
      87      134300 :   _fs.thermophysicalProperties(_liquid_porepressure[_qp], _enthalpy[_qp], _qp, _fsp);
      88      134300 : }
      89             : 
      90             : template <bool is_ad>
      91             : void
      92        7140 : PorousFlowFluidStateSingleComponentTempl<is_ad>::initQpStatefulProperties()
      93             : {
      94        7140 :   _is_initqp = true;
      95             :   // Set the size of pressure and saturation vectors
      96        7140 :   PorousFlowFluidStateBaseMaterialTempl<is_ad>::initQpStatefulProperties();
      97             : 
      98             :   // Set the initial values of the temperature at the nodes.
      99             :   // Note: not required for qp materials as no old values at the qps are requested
     100        7140 :   if (_nodal_material)
     101             :   {
     102             :     // Temperature doesn't depend on fluid phase
     103        3570 :     _temperature[_qp] = genericValue(_fsp[_aqueous_phase_number].temperature) - _T_c2k;
     104             :   }
     105        7140 : }
     106             : 
     107             : template <bool is_ad>
     108             : void
     109      127160 : PorousFlowFluidStateSingleComponentTempl<is_ad>::computeQpProperties()
     110             : {
     111             : 
     112      127160 :   PorousFlowFluidStateBaseMaterialTempl<is_ad>::computeQpProperties();
     113             : 
     114             :   // Temperature doesn't depend on fluid phase
     115      127160 :   _temperature[_qp] = genericValue(_fsp[_aqueous_phase_number].temperature) - _T_c2k;
     116             : 
     117             :   if (!is_ad)
     118             :   {
     119      127160 :     (*_dtemperature_dvar)[_qp][_pvar] =
     120      127160 :         _fsp[_aqueous_phase_number].temperature.derivatives()[_pidx];
     121      127160 :     (*_dtemperature_dvar)[_qp][_hvar] =
     122      127160 :         _fsp[_aqueous_phase_number].temperature.derivatives()[_hidx];
     123             :   }
     124             : 
     125             :   // Derivative of pressure, saturation and fluid properties wrt variables
     126             :   if (!is_ad)
     127      381480 :     for (unsigned int ph = 0; ph < _num_phases; ++ph)
     128             :     {
     129      385423 :       (*_dporepressure_dvar)[_qp][ph][_pvar] = _fsp[ph].pressure.derivatives()[_pidx];
     130      254320 :       (*_dporepressure_dvar)[_qp][ph][_hvar] = _fsp[ph].pressure.derivatives()[_hidx];
     131             : 
     132      262224 :       (*_dsaturation_dvar)[_qp][ph][_pvar] = _fsp[ph].saturation.derivatives()[_pidx];
     133      254320 :       (*_dsaturation_dvar)[_qp][ph][_hvar] = _fsp[ph].saturation.derivatives()[_hidx];
     134             : 
     135      385423 :       (*_dfluid_density_dvar)[_qp][ph][_pvar] = _fsp[ph].density.derivatives()[_pidx];
     136      254320 :       (*_dfluid_density_dvar)[_qp][ph][_hvar] = _fsp[ph].density.derivatives()[_hidx];
     137             : 
     138      385423 :       (*_dfluid_viscosity_dvar)[_qp][ph][_pvar] = _fsp[ph].viscosity.derivatives()[_pidx];
     139      254320 :       (*_dfluid_viscosity_dvar)[_qp][ph][_hvar] = _fsp[ph].viscosity.derivatives()[_hidx];
     140             : 
     141      262224 :       (*_dfluid_enthalpy_dvar)[_qp][ph][_pvar] = _fsp[ph].enthalpy.derivatives()[_pidx];
     142      254320 :       (*_dfluid_enthalpy_dvar)[_qp][ph][_hvar] = _fsp[ph].enthalpy.derivatives()[_hidx];
     143             : 
     144      254320 :       (*_dfluid_internal_energy_dvar)[_qp][ph][_pvar] =
     145      254320 :           _fsp[ph].internal_energy.derivatives()[_pidx];
     146      254320 :       (*_dfluid_internal_energy_dvar)[_qp][ph][_hvar] =
     147             :           _fsp[ph].internal_energy.derivatives()[_hidx];
     148             :     }
     149             : 
     150             :   // If the material properties are being evaluated at the qps, calculate the
     151             :   // gradients as well. Note: only nodal properties are evaluated in
     152             :   // initQpStatefulProperties(), so no need to check _is_initqp flag for qp
     153             :   // properties
     154      127160 :   if (!_nodal_material)
     155             :     if constexpr (!is_ad)
     156             :     {
     157             :       // Need to compute second derivatives of properties wrt variables for some of
     158             :       // the gradient derivatives. Use finite differences for now
     159       63454 :       const Real dp = 1.0e-5 * _liquid_porepressure[_qp];
     160       63454 :       const Real dh = 1.0e-5 * _enthalpy[_qp];
     161             : 
     162       63454 :       std::vector<FluidStateProperties> fsp_dp(_num_phases, FluidStateProperties(_num_components));
     163       63454 :       _fs.thermophysicalProperties(_liquid_porepressure[_qp] + dp, _enthalpy[_qp], _qp, fsp_dp);
     164             : 
     165       63454 :       std::vector<FluidStateProperties> fsp_dh(_num_phases, FluidStateProperties(_num_components));
     166       63454 :       _fs.thermophysicalProperties(_liquid_porepressure[_qp], _enthalpy[_qp] + dh, _qp, fsp_dh);
     167             : 
     168             :       // Gradient of temperature (non-zero in all phases)
     169       63454 :       (*_grad_temperature_qp)[_qp] = (*_dtemperature_dvar)[_qp][_pvar] * _liquid_gradp_qp[_qp] +
     170       63454 :                                      (*_dtemperature_dvar)[_qp][_hvar] * _gradh_qp[_qp];
     171       63454 :       (*_dgrad_temperature_dgradv)[_qp][_pvar] = (*_dtemperature_dvar)[_qp][_pvar];
     172       63454 :       (*_dgrad_temperature_dgradv)[_qp][_hvar] = (*_dtemperature_dvar)[_qp][_hvar];
     173             : 
     174      100274 :       const auto d2T_dp2 = (fsp_dp[_aqueous_phase_number].temperature.derivatives()[_pidx] -
     175      126905 :                             _fsp[_aqueous_phase_number].temperature.derivatives()[_pidx]) /
     176             :                            dp;
     177             : 
     178       99010 :       const auto d2T_dh2 = (fsp_dh[_aqueous_phase_number].temperature.derivatives()[_hidx] -
     179       63454 :                             _fsp[_aqueous_phase_number].temperature.derivatives()[_hidx]) /
     180             :                            dh;
     181             : 
     182       63454 :       const auto d2T_dph = (fsp_dp[_aqueous_phase_number].temperature.derivatives()[_hidx] -
     183       63454 :                             _fsp[_aqueous_phase_number].temperature.derivatives()[_hidx]) /
     184       63454 :                                dp +
     185       63454 :                            (fsp_dh[_aqueous_phase_number].temperature.derivatives()[_pidx] -
     186       63454 :                             _fsp[_aqueous_phase_number].temperature.derivatives()[_pidx]) /
     187             :                                dh;
     188             : 
     189       63454 :       (*_dgrad_temperature_dv)[_qp][_pvar] =
     190             :           d2T_dp2 * _liquid_gradp_qp[_qp] + d2T_dph * _gradh_qp[_qp];
     191       63454 :       (*_dgrad_temperature_dv)[_qp][_hvar] =
     192             :           d2T_dph * _liquid_gradp_qp[_qp] + d2T_dh2 * _gradh_qp[_qp];
     193             : 
     194             :       // Gradient of saturation and derivatives
     195       63454 :       (*_grads_qp)[_qp][_gas_phase_number] =
     196             :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_pvar] * _liquid_gradp_qp[_qp] +
     197       63454 :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_hvar] * _gradh_qp[_qp];
     198       63454 :       (*_grads_qp)[_qp][_aqueous_phase_number] = -(*_grads_qp)[_qp][_gas_phase_number];
     199             : 
     200       63454 :       (*_dgrads_qp_dgradv)[_qp][_gas_phase_number][_pvar] =
     201             :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_pvar];
     202       63454 :       (*_dgrads_qp_dgradv)[_qp][_aqueous_phase_number][_pvar] =
     203             :           -(*_dgrads_qp_dgradv)[_qp][_gas_phase_number][_pvar];
     204             : 
     205       63454 :       (*_dgrads_qp_dgradv)[_qp][_gas_phase_number][_hvar] =
     206             :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_hvar];
     207       63454 :       (*_dgrads_qp_dgradv)[_qp][_aqueous_phase_number][_hvar] =
     208             :           -(*_dgrads_qp_dgradv)[_qp][_gas_phase_number][_hvar];
     209             : 
     210       64718 :       const Real d2s_dp2 = (fsp_dp[_gas_phase_number].saturation.derivatives()[_pidx] -
     211       65430 :                             _fsp[_gas_phase_number].saturation.derivatives()[_pidx]) /
     212             :                            dp;
     213             : 
     214       64718 :       const Real d2s_dh2 = (fsp_dh[_gas_phase_number].saturation.derivatives()[_hidx] -
     215       63454 :                             _fsp[_gas_phase_number].saturation.derivatives()[_hidx]) /
     216             :                            dh;
     217             : 
     218       63454 :       const Real d2s_dph = (fsp_dp[_gas_phase_number].saturation.derivatives()[_hidx] -
     219       63454 :                             _fsp[_gas_phase_number].saturation.derivatives()[_hidx]) /
     220       63454 :                                dp +
     221       63454 :                            (fsp_dh[_gas_phase_number].saturation.derivatives()[_pidx] -
     222       63454 :                             _fsp[_gas_phase_number].saturation.derivatives()[_pidx]) /
     223             :                                dh;
     224             : 
     225       63454 :       (*_dgrads_qp_dv)[_qp][_gas_phase_number][_pvar] =
     226             :           d2s_dp2 * _liquid_gradp_qp[_qp] + d2s_dph * _gradh_qp[_qp];
     227       63454 :       (*_dgrads_qp_dv)[_qp][_aqueous_phase_number][_pvar] =
     228             :           -(*_dgrads_qp_dv)[_qp][_gas_phase_number][_pvar];
     229             : 
     230       63454 :       (*_dgrads_qp_dv)[_qp][_gas_phase_number][_hvar] =
     231             :           d2s_dh2 * _gradh_qp[_qp] + d2s_dph * _liquid_gradp_qp[_qp];
     232       63454 :       (*_dgrads_qp_dv)[_qp][_aqueous_phase_number][_hvar] =
     233             :           -(*_dgrads_qp_dv)[_qp][_gas_phase_number][_hvar];
     234             : 
     235             :       // Gradient of porepressure and derivatives
     236             :       // Note: need first and second derivativea of capillary pressure
     237       63454 :       const Real dpc = _pc.dCapillaryPressure(_fsp[_aqueous_phase_number].saturation.value());
     238       63454 :       const Real d2pc = _pc.d2CapillaryPressure(_fsp[_aqueous_phase_number].saturation.value());
     239             : 
     240       63454 :       (*_gradp_qp)[_qp][_aqueous_phase_number] = _liquid_gradp_qp[_qp];
     241       63454 :       (*_gradp_qp)[_qp][_gas_phase_number] =
     242       63454 :           _liquid_gradp_qp[_qp] + dpc * (*_grads_qp)[_qp][_aqueous_phase_number];
     243             : 
     244      190362 :       for (unsigned int ph = 0; ph < _num_phases; ++ph)
     245      126908 :         (*_dgradp_qp_dgradv)[_qp][ph][_pvar] = 1.0;
     246             : 
     247       63454 :       (*_dgradp_qp_dgradv)[_qp][_gas_phase_number][_pvar] +=
     248       63454 :           dpc * (*_dgrads_qp_dgradv)[_qp][_aqueous_phase_number][_pvar];
     249       63454 :       (*_dgradp_qp_dgradv)[_qp][_gas_phase_number][_hvar] =
     250       63454 :           dpc * (*_dgrads_qp_dgradv)[_qp][_aqueous_phase_number][_hvar];
     251             : 
     252       63454 :       (*_dgradp_qp_dv)[_qp][_gas_phase_number][_pvar] =
     253             :           d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     254       63454 :               (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_pvar] +
     255       63454 :           dpc * (*_dgrads_qp_dv)[_qp][_aqueous_phase_number][_pvar];
     256             : 
     257       63454 :       (*_dgradp_qp_dv)[_qp][_gas_phase_number][_hvar] =
     258             :           d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     259             :               (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_hvar] +
     260             :           dpc * (*_dgrads_qp_dv)[_qp][_aqueous_phase_number][_hvar];
     261       63454 :     }
     262      127160 : }
     263             : 
     264             : template <bool is_ad>
     265             : void
     266      134300 : PorousFlowFluidStateSingleComponentTempl<is_ad>::setMaterialVectorSize() const
     267             : {
     268      134300 :   PorousFlowFluidStateBaseMaterialTempl<is_ad>::setMaterialVectorSize();
     269             : 
     270             :   // Derivatives and gradients are not required in initQpStatefulProperties
     271      134300 :   if (!_is_initqp)
     272             :   {
     273             :     // Temperature doesn't depend of fluid phase
     274      127160 :     (*_dtemperature_dvar)[_qp].assign(_num_pf_vars, 0.0);
     275             : 
     276             :     // The gradient of the temperature is needed for qp materials and AD materials
     277      127160 :     if (!_nodal_material || is_ad)
     278       63454 :       (*_grad_temperature_qp)[_qp] = RealGradient();
     279             : 
     280             :     // No derivatives are required for AD materials
     281             :     if (!is_ad)
     282      127160 :       if (!_nodal_material)
     283             :       {
     284       63454 :         (*_dgrad_temperature_dgradv)[_qp].assign(_num_pf_vars, 0.0);
     285       63454 :         (*_dgrad_temperature_dv)[_qp].assign(_num_pf_vars, RealGradient());
     286             :       }
     287             :   }
     288      134300 : }
     289             : 
     290             : template class PorousFlowFluidStateSingleComponentTempl<false>;
     291             : template class PorousFlowFluidStateSingleComponentTempl<true>;