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         546 : PorousFlowFluidStateSingleComponentTempl<is_ad>::validParams()
      19             : {
      20         546 :   InputParameters params = PorousFlowFluidStateBaseMaterialTempl<is_ad>::validParams();
      21        1092 :   params.addRequiredCoupledVar("porepressure",
      22             :                                "Variable that is the porepressure of the liquid phase");
      23        1092 :   params.addRequiredCoupledVar("enthalpy", "Enthalpy of the fluid");
      24        1092 :   params.addRequiredParam<UserObjectName>("fluid_state", "Name of the FluidState UserObject");
      25         546 :   params.addClassDescription(
      26             :       "Class for single component multiphase fluid state calculations using pressure and enthalpy");
      27         546 :   return params;
      28           0 : }
      29             : 
      30             : template <bool is_ad>
      31         420 : PorousFlowFluidStateSingleComponentTempl<is_ad>::PorousFlowFluidStateSingleComponentTempl(
      32             :     const InputParameters & parameters)
      33             :   : PorousFlowFluidStateBaseMaterialTempl<is_ad>(parameters),
      34         840 :     _liquid_porepressure(_nodal_material
      35         210 :                              ? this->template coupledGenericDofValue<is_ad>("porepressure")
      36         630 :                              : this->template coupledGenericValue<is_ad>("porepressure")),
      37         420 :     _liquid_gradp_qp(this->template coupledGenericGradient<is_ad>("porepressure")),
      38         420 :     _liquid_porepressure_varnum(coupled("porepressure")),
      39         420 :     _pvar(_dictator.isPorousFlowVariable(_liquid_porepressure_varnum)
      40         420 :               ? _dictator.porousFlowVariableNum(_liquid_porepressure_varnum)
      41             :               : 0),
      42         420 :     _enthalpy(_nodal_material ? this->template coupledGenericDofValue<is_ad>("enthalpy")
      43         630 :                               : this->template coupledGenericValue<is_ad>("enthalpy")),
      44         420 :     _gradh_qp(this->template coupledGenericGradient<is_ad>("enthalpy")),
      45         420 :     _enthalpy_varnum(coupled("enthalpy")),
      46         420 :     _hvar(_dictator.isPorousFlowVariable(_enthalpy_varnum)
      47         420 :               ? _dictator.porousFlowVariableNum(_enthalpy_varnum)
      48             :               : 0),
      49         420 :     _fs(this->template getUserObject<PorousFlowFluidStateSingleComponentBase>("fluid_state")),
      50         420 :     _aqueous_phase_number(_fs.aqueousPhaseIndex()),
      51         420 :     _gas_phase_number(_fs.gasPhaseIndex()),
      52         420 :     _temperature(
      53         420 :         this->template declareGenericProperty<Real, is_ad>("PorousFlow_temperature" + _sfx)),
      54         840 :     _grad_temperature_qp(_nodal_material
      55         420 :                              ? nullptr
      56         630 :                              : &this->template declareGenericProperty<RealGradient, is_ad>(
      57             :                                    "PorousFlow_grad_temperature_qp")),
      58         420 :     _dtemperature_dvar(is_ad ? nullptr
      59         840 :                              : &this->template declareProperty<std::vector<Real>>(
      60             :                                    "dPorousFlow_temperature" + _sfx + "_dvar")),
      61         840 :     _dgrad_temperature_dgradv(is_ad || _nodal_material
      62         420 :                                   ? nullptr
      63         420 :                                   : &this->template declareProperty<std::vector<Real>>(
      64             :                                         "dPorousFlow_grad_temperature_qp_dgradvar")),
      65         840 :     _dgrad_temperature_dv(is_ad ? nullptr
      66         420 :                           : _nodal_material
      67             :                               ? nullptr
      68         420 :                               : &this->template declareProperty<std::vector<RealGradient>>(
      69             :                                     "dPorousFlow_grad_temperature_qp_dvar")),
      70         420 :     _pidx(_fs.getPressureIndex()),
      71         840 :     _hidx(_fs.getEnthalpyIndex())
      72             : {
      73             :   // Check that the number of phases in the fluidstate class is also provided in the Dictator
      74         420 :   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         420 : }
      81             : 
      82             : template <bool is_ad>
      83             : void
      84       89485 : PorousFlowFluidStateSingleComponentTempl<is_ad>::thermophysicalProperties()
      85             : {
      86       89485 :   _fs.clearFluidStateProperties(_fsp);
      87       89485 :   _fs.thermophysicalProperties(_liquid_porepressure[_qp], _enthalpy[_qp], _qp, _fsp);
      88       89485 : }
      89             : 
      90             : template <bool is_ad>
      91             : void
      92        4420 : PorousFlowFluidStateSingleComponentTempl<is_ad>::initQpStatefulProperties()
      93             : {
      94        4420 :   _is_initqp = true;
      95             :   // Set the size of pressure and saturation vectors
      96        4420 :   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        4420 :   if (_nodal_material)
     101             :   {
     102             :     // Temperature doesn't depend on fluid phase
     103        2210 :     _temperature[_qp] = genericValue(_fsp[_aqueous_phase_number].temperature) - _T_c2k;
     104             :   }
     105        4420 : }
     106             : 
     107             : template <bool is_ad>
     108             : void
     109       85065 : PorousFlowFluidStateSingleComponentTempl<is_ad>::computeQpProperties()
     110             : {
     111             : 
     112       85065 :   PorousFlowFluidStateBaseMaterialTempl<is_ad>::computeQpProperties();
     113             : 
     114             :   // Temperature doesn't depend on fluid phase
     115       85065 :   _temperature[_qp] = genericValue(_fsp[_aqueous_phase_number].temperature) - _T_c2k;
     116             : 
     117             :   if (!is_ad)
     118             :   {
     119       85065 :     (*_dtemperature_dvar)[_qp][_pvar] =
     120       85065 :         _fsp[_aqueous_phase_number].temperature.derivatives()[_pidx];
     121       85065 :     (*_dtemperature_dvar)[_qp][_hvar] =
     122       85065 :         _fsp[_aqueous_phase_number].temperature.derivatives()[_hidx];
     123             :   }
     124             : 
     125             :   // Derivative of pressure, saturation and fluid properties wrt variables
     126             :   if (!is_ad)
     127      255195 :     for (unsigned int ph = 0; ph < _num_phases; ++ph)
     128             :     {
     129      340260 :       (*_dporepressure_dvar)[_qp][ph][_pvar] = _fsp[ph].pressure.derivatives()[_pidx];
     130      340260 :       (*_dporepressure_dvar)[_qp][ph][_hvar] = _fsp[ph].pressure.derivatives()[_hidx];
     131             : 
     132      340260 :       (*_dsaturation_dvar)[_qp][ph][_pvar] = _fsp[ph].saturation.derivatives()[_pidx];
     133      340260 :       (*_dsaturation_dvar)[_qp][ph][_hvar] = _fsp[ph].saturation.derivatives()[_hidx];
     134             : 
     135      340260 :       (*_dfluid_density_dvar)[_qp][ph][_pvar] = _fsp[ph].density.derivatives()[_pidx];
     136      340260 :       (*_dfluid_density_dvar)[_qp][ph][_hvar] = _fsp[ph].density.derivatives()[_hidx];
     137             : 
     138      340260 :       (*_dfluid_viscosity_dvar)[_qp][ph][_pvar] = _fsp[ph].viscosity.derivatives()[_pidx];
     139      340260 :       (*_dfluid_viscosity_dvar)[_qp][ph][_hvar] = _fsp[ph].viscosity.derivatives()[_hidx];
     140             : 
     141      340260 :       (*_dfluid_enthalpy_dvar)[_qp][ph][_pvar] = _fsp[ph].enthalpy.derivatives()[_pidx];
     142      340260 :       (*_dfluid_enthalpy_dvar)[_qp][ph][_hvar] = _fsp[ph].enthalpy.derivatives()[_hidx];
     143             : 
     144      170130 :       (*_dfluid_internal_energy_dvar)[_qp][ph][_pvar] =
     145      170130 :           _fsp[ph].internal_energy.derivatives()[_pidx];
     146      170130 :       (*_dfluid_internal_energy_dvar)[_qp][ph][_hvar] =
     147      170130 :           _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       85065 :   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       42449 :       const Real dp = 1.0e-5 * _liquid_porepressure[_qp];
     160       42449 :       const Real dh = 1.0e-5 * _enthalpy[_qp];
     161             : 
     162       42449 :       std::vector<FluidStateProperties> fsp_dp(_num_phases, FluidStateProperties(_num_components));
     163       42449 :       _fs.thermophysicalProperties(_liquid_porepressure[_qp] + dp, _enthalpy[_qp], _qp, fsp_dp);
     164             : 
     165       42449 :       std::vector<FluidStateProperties> fsp_dh(_num_phases, FluidStateProperties(_num_components));
     166       42449 :       _fs.thermophysicalProperties(_liquid_porepressure[_qp], _enthalpy[_qp] + dh, _qp, fsp_dh);
     167             : 
     168             :       // Gradient of temperature (non-zero in all phases)
     169       42449 :       (*_grad_temperature_qp)[_qp] = (*_dtemperature_dvar)[_qp][_pvar] * _liquid_gradp_qp[_qp] +
     170       42449 :                                      (*_dtemperature_dvar)[_qp][_hvar] * _gradh_qp[_qp];
     171       42449 :       (*_dgrad_temperature_dgradv)[_qp][_pvar] = (*_dtemperature_dvar)[_qp][_pvar];
     172       42449 :       (*_dgrad_temperature_dgradv)[_qp][_hvar] = (*_dtemperature_dvar)[_qp][_hvar];
     173             : 
     174       42449 :       const auto d2T_dp2 = (fsp_dp[_aqueous_phase_number].temperature.derivatives()[_pidx] -
     175       42449 :                             _fsp[_aqueous_phase_number].temperature.derivatives()[_pidx]) /
     176             :                            dp;
     177             : 
     178       42449 :       const auto d2T_dh2 = (fsp_dh[_aqueous_phase_number].temperature.derivatives()[_hidx] -
     179       42449 :                             _fsp[_aqueous_phase_number].temperature.derivatives()[_hidx]) /
     180             :                            dh;
     181             : 
     182       42449 :       const auto d2T_dph = (fsp_dp[_aqueous_phase_number].temperature.derivatives()[_hidx] -
     183       42449 :                             _fsp[_aqueous_phase_number].temperature.derivatives()[_hidx]) /
     184       42449 :                                dp +
     185       42449 :                            (fsp_dh[_aqueous_phase_number].temperature.derivatives()[_pidx] -
     186       42449 :                             _fsp[_aqueous_phase_number].temperature.derivatives()[_pidx]) /
     187             :                                dh;
     188             : 
     189       42449 :       (*_dgrad_temperature_dv)[_qp][_pvar] =
     190       42449 :           d2T_dp2 * _liquid_gradp_qp[_qp] + d2T_dph * _gradh_qp[_qp];
     191       42449 :       (*_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       42449 :       (*_grads_qp)[_qp][_gas_phase_number] =
     196             :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_pvar] * _liquid_gradp_qp[_qp] +
     197       42449 :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_hvar] * _gradh_qp[_qp];
     198       42449 :       (*_grads_qp)[_qp][_aqueous_phase_number] = -(*_grads_qp)[_qp][_gas_phase_number];
     199             : 
     200       42449 :       (*_dgrads_qp_dgradv)[_qp][_gas_phase_number][_pvar] =
     201             :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_pvar];
     202       42449 :       (*_dgrads_qp_dgradv)[_qp][_aqueous_phase_number][_pvar] =
     203             :           -(*_dgrads_qp_dgradv)[_qp][_gas_phase_number][_pvar];
     204             : 
     205       42449 :       (*_dgrads_qp_dgradv)[_qp][_gas_phase_number][_hvar] =
     206             :           (*_dsaturation_dvar)[_qp][_gas_phase_number][_hvar];
     207       42449 :       (*_dgrads_qp_dgradv)[_qp][_aqueous_phase_number][_hvar] =
     208             :           -(*_dgrads_qp_dgradv)[_qp][_gas_phase_number][_hvar];
     209             : 
     210       42449 :       const Real d2s_dp2 = (fsp_dp[_gas_phase_number].saturation.derivatives()[_pidx] -
     211       42449 :                             _fsp[_gas_phase_number].saturation.derivatives()[_pidx]) /
     212             :                            dp;
     213             : 
     214       42449 :       const Real d2s_dh2 = (fsp_dh[_gas_phase_number].saturation.derivatives()[_hidx] -
     215       42449 :                             _fsp[_gas_phase_number].saturation.derivatives()[_hidx]) /
     216             :                            dh;
     217             : 
     218       42449 :       const Real d2s_dph = (fsp_dp[_gas_phase_number].saturation.derivatives()[_hidx] -
     219       42449 :                             _fsp[_gas_phase_number].saturation.derivatives()[_hidx]) /
     220       42449 :                                dp +
     221       42449 :                            (fsp_dh[_gas_phase_number].saturation.derivatives()[_pidx] -
     222       42449 :                             _fsp[_gas_phase_number].saturation.derivatives()[_pidx]) /
     223             :                                dh;
     224             : 
     225       42449 :       (*_dgrads_qp_dv)[_qp][_gas_phase_number][_pvar] =
     226       42449 :           d2s_dp2 * _liquid_gradp_qp[_qp] + d2s_dph * _gradh_qp[_qp];
     227       42449 :       (*_dgrads_qp_dv)[_qp][_aqueous_phase_number][_pvar] =
     228             :           -(*_dgrads_qp_dv)[_qp][_gas_phase_number][_pvar];
     229             : 
     230       42449 :       (*_dgrads_qp_dv)[_qp][_gas_phase_number][_hvar] =
     231             :           d2s_dh2 * _gradh_qp[_qp] + d2s_dph * _liquid_gradp_qp[_qp];
     232       42449 :       (*_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       42449 :       const Real dpc = _pc.dCapillaryPressure(_fsp[_aqueous_phase_number].saturation.value());
     238       42449 :       const Real d2pc = _pc.d2CapillaryPressure(_fsp[_aqueous_phase_number].saturation.value());
     239             : 
     240       42449 :       (*_gradp_qp)[_qp][_aqueous_phase_number] = _liquid_gradp_qp[_qp];
     241       42449 :       (*_gradp_qp)[_qp][_gas_phase_number] =
     242       42449 :           _liquid_gradp_qp[_qp] + dpc * (*_grads_qp)[_qp][_aqueous_phase_number];
     243             : 
     244      127347 :       for (unsigned int ph = 0; ph < _num_phases; ++ph)
     245       84898 :         (*_dgradp_qp_dgradv)[_qp][ph][_pvar] = 1.0;
     246             : 
     247       42449 :       (*_dgradp_qp_dgradv)[_qp][_gas_phase_number][_pvar] +=
     248       42449 :           dpc * (*_dgrads_qp_dgradv)[_qp][_aqueous_phase_number][_pvar];
     249       42449 :       (*_dgradp_qp_dgradv)[_qp][_gas_phase_number][_hvar] =
     250       42449 :           dpc * (*_dgrads_qp_dgradv)[_qp][_aqueous_phase_number][_hvar];
     251             : 
     252       42449 :       (*_dgradp_qp_dv)[_qp][_gas_phase_number][_pvar] =
     253             :           d2pc * (*_grads_qp)[_qp][_aqueous_phase_number] *
     254       42449 :               (*_dsaturation_dvar)[_qp][_aqueous_phase_number][_pvar] +
     255       42449 :           dpc * (*_dgrads_qp_dv)[_qp][_aqueous_phase_number][_pvar];
     256             : 
     257       42449 :       (*_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       42449 :     }
     262       85065 : }
     263             : 
     264             : template <bool is_ad>
     265             : void
     266       89485 : PorousFlowFluidStateSingleComponentTempl<is_ad>::setMaterialVectorSize() const
     267             : {
     268       89485 :   PorousFlowFluidStateBaseMaterialTempl<is_ad>::setMaterialVectorSize();
     269             : 
     270             :   // Derivatives and gradients are not required in initQpStatefulProperties
     271       89485 :   if (!_is_initqp)
     272             :   {
     273             :     // Temperature doesn't depend of fluid phase
     274       85065 :     (*_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       85065 :     if (!_nodal_material || is_ad)
     278       42449 :       (*_grad_temperature_qp)[_qp] = RealGradient();
     279             : 
     280             :     // No derivatives are required for AD materials
     281             :     if (!is_ad)
     282       85065 :       if (!_nodal_material)
     283             :       {
     284       42449 :         (*_dgrad_temperature_dgradv)[_qp].assign(_num_pf_vars, 0.0);
     285       42449 :         (*_dgrad_temperature_dv)[_qp].assign(_num_pf_vars, RealGradient());
     286             :       }
     287             :   }
     288       89485 : }
     289             : 
     290             : template class PorousFlowFluidStateSingleComponentTempl<false>;
     291             : template class PorousFlowFluidStateSingleComponentTempl<true>;