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 "PorousFlowDispersiveFlux.h"
      11             : 
      12             : #include "MooseVariable.h"
      13             : 
      14             : registerMooseObject("PorousFlowApp", PorousFlowDispersiveFlux);
      15             : 
      16             : InputParameters
      17        1232 : PorousFlowDispersiveFlux::validParams()
      18             : {
      19        1232 :   InputParameters params = Kernel::validParams();
      20        2464 :   params.addParam<unsigned int>(
      21        2464 :       "fluid_component", 0, "The index corresponding to the fluid component for this kernel");
      22        2464 :   params.addRequiredParam<UserObjectName>(
      23             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      24        2464 :   params.addRequiredParam<std::vector<Real>>(
      25             :       "disp_long", "Vector of longitudinal dispersion coefficients for each phase");
      26        2464 :   params.addRequiredParam<std::vector<Real>>(
      27             :       "disp_trans", "Vector of transverse dispersion coefficients for each phase");
      28        2464 :   params.addRequiredParam<RealVectorValue>("gravity",
      29             :                                            "Gravitational acceleration vector downwards (m/s^2)");
      30        1232 :   params.addClassDescription(
      31             :       "Dispersive and diffusive flux of the component given by fluid_component in all phases");
      32        1232 :   return params;
      33           0 : }
      34             : 
      35         676 : PorousFlowDispersiveFlux::PorousFlowDispersiveFlux(const InputParameters & parameters)
      36             :   : Kernel(parameters),
      37             : 
      38         676 :     _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      39        1352 :     _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      40             :         "dPorousFlow_fluid_phase_density_qp_dvar")),
      41        1352 :     _grad_mass_frac(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
      42             :         "PorousFlow_grad_mass_frac_qp")),
      43        1352 :     _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      44             :         "dPorousFlow_mass_frac_qp_dvar")),
      45        1352 :     _porosity_qp(getMaterialProperty<Real>("PorousFlow_porosity_qp")),
      46        1352 :     _dporosity_qp_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_qp_dvar")),
      47        1352 :     _tortuosity(getMaterialProperty<std::vector<Real>>("PorousFlow_tortuosity_qp")),
      48         676 :     _dtortuosity_dvar(
      49         676 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_tortuosity_qp_dvar")),
      50         676 :     _diffusion_coeff(
      51         676 :         getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_diffusion_coeff_qp")),
      52        1352 :     _ddiffusion_coeff_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      53             :         "dPorousFlow_diffusion_coeff_qp_dvar")),
      54         676 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      55        1352 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      56         676 :     _num_phases(_dictator.numPhases()),
      57         676 :     _identity_tensor(RankTwoTensor::initIdentity),
      58         676 :     _relative_permeability(
      59         676 :         getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      60        1352 :     _drelative_permeability_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      61             :         "dPorousFlow_relative_permeability_qp_dvar")),
      62        1352 :     _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      63         676 :     _dfluid_viscosity_dvar(
      64         676 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_qp_dvar")),
      65        1352 :     _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      66         676 :     _dpermeability_dvar(
      67         676 :         getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
      68        1352 :     _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
      69             :         "dPorousFlow_permeability_qp_dgradvar")),
      70        1352 :     _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
      71        1352 :     _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
      72             :         "dPorousFlow_grad_porepressure_qp_dgradvar")),
      73        1352 :     _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
      74             :         "dPorousFlow_grad_porepressure_qp_dvar")),
      75        1352 :     _gravity(getParam<RealVectorValue>("gravity")),
      76        1352 :     _disp_long(getParam<std::vector<Real>>("disp_long")),
      77        1352 :     _disp_trans(getParam<std::vector<Real>>("disp_trans")),
      78        1352 :     _perm_derivs(_dictator.usePermDerivs())
      79             : {
      80         676 :   if (_fluid_component >= _dictator.numComponents())
      81           0 :     paramError(
      82             :         "fluid_component",
      83             :         "The Dictator proclaims that the maximum fluid component index in this simulation is ",
      84           0 :         _dictator.numComponents() - 1,
      85             :         " whereas you have used ",
      86           0 :         _fluid_component,
      87             :         ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
      88             : 
      89             :   // Check that sufficient values of the dispersion coefficients have been entered
      90         676 :   if (_disp_long.size() != _num_phases)
      91           0 :     paramError(
      92             :         "disp_long",
      93             :         "The number of longitudinal dispersion coefficients is not equal to the number of phases");
      94             : 
      95         676 :   if (_disp_trans.size() != _num_phases)
      96           0 :     paramError("disp_trans",
      97             :                "The number of transverse dispersion coefficients disp_trans is not equal to the "
      98             :                "number of phases");
      99         676 : }
     100             : 
     101             : Real
     102     5553920 : PorousFlowDispersiveFlux::computeQpResidual()
     103             : {
     104             :   RealVectorValue flux = 0.0;
     105             :   RealVectorValue velocity;
     106             :   Real velocity_abs;
     107     5553920 :   RankTwoTensor v2;
     108     5553920 :   RankTwoTensor dispersion;
     109             :   dispersion.zero();
     110             :   Real diffusion;
     111             : 
     112    11362880 :   for (unsigned int ph = 0; ph < _num_phases; ++ph)
     113             :   {
     114             :     // Diffusive component
     115     5808960 :     diffusion =
     116     5808960 :         _porosity_qp[_qp] * _tortuosity[_qp][ph] * _diffusion_coeff[_qp][ph][_fluid_component];
     117             : 
     118             :     // Calculate Darcy velocity
     119     5808960 :     velocity = (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity) *
     120     5808960 :                 _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph]);
     121     5808960 :     velocity_abs = std::sqrt(velocity * velocity);
     122             : 
     123     5808960 :     if (velocity_abs > 0.0)
     124             :     {
     125     5111776 :       v2 = RankTwoTensor::selfOuterProduct(velocity);
     126             : 
     127             :       // Add longitudinal dispersion to diffusive component
     128     5111776 :       diffusion += _disp_trans[ph] * velocity_abs;
     129     5111776 :       dispersion = (_disp_long[ph] - _disp_trans[ph]) * v2 / velocity_abs;
     130             :     }
     131             : 
     132     5808960 :     flux += _fluid_density_qp[_qp][ph] * (diffusion * _identity_tensor + dispersion) *
     133     5808960 :             _grad_mass_frac[_qp][ph][_fluid_component];
     134             :   }
     135     5553920 :   return _grad_test[_i][_qp] * flux;
     136             : }
     137             : 
     138             : Real
     139    14157760 : PorousFlowDispersiveFlux::computeQpJacobian()
     140             : {
     141    14157760 :   return computeQpJac(_var.number());
     142             : }
     143             : 
     144             : Real
     145    14157760 : PorousFlowDispersiveFlux::computeQpOffDiagJacobian(unsigned int jvar)
     146             : {
     147    14157760 :   return computeQpJac(jvar);
     148             : }
     149             : 
     150             : Real
     151    28315520 : PorousFlowDispersiveFlux::computeQpJac(unsigned int jvar) const
     152             : {
     153             :   // If the variable is not a valid PorousFlow variable, set the Jacobian to 0
     154    28315520 :   if (_dictator.notPorousFlowVariable(jvar))
     155             :     return 0.0;
     156             : 
     157    28315520 :   const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
     158             : 
     159             :   RealVectorValue velocity;
     160             :   Real velocity_abs;
     161    28315520 :   RankTwoTensor v2;
     162    28315520 :   RankTwoTensor dispersion;
     163             :   dispersion.zero();
     164             :   Real diffusion;
     165             :   RealVectorValue dflux = 0.0;
     166             : 
     167    58292480 :   for (unsigned int ph = 0; ph < _num_phases; ++ph)
     168             :   {
     169             :     // Diffusive component
     170    29976960 :     diffusion =
     171    29976960 :         _porosity_qp[_qp] * _tortuosity[_qp][ph] * _diffusion_coeff[_qp][ph][_fluid_component];
     172             : 
     173             :     // Calculate Darcy velocity
     174    29976960 :     velocity = (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity) *
     175    29976960 :                 _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph]);
     176    29976960 :     velocity_abs = std::sqrt(velocity * velocity);
     177             : 
     178    29976960 :     if (velocity_abs > 0.0)
     179             :     {
     180    26824448 :       v2 = RankTwoTensor::selfOuterProduct(velocity);
     181             : 
     182             :       // Add longitudinal dispersion to diffusive component
     183    26824448 :       diffusion += _disp_trans[ph] * velocity_abs;
     184    26824448 :       dispersion = (_disp_long[ph] - _disp_trans[ph]) * v2 / velocity_abs;
     185             :     }
     186             : 
     187             :     // Derivative of Darcy velocity
     188             :     RealVectorValue dvelocity =
     189    29976960 :         _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
     190    29976960 :                               _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
     191    29976960 :     dvelocity += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
     192             : 
     193    29976960 :     if (_perm_derivs)
     194             :     {
     195           0 :       dvelocity += _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
     196           0 :                    (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
     197             : 
     198           0 :       for (const auto i : make_range(Moose::dim))
     199           0 :         dvelocity += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
     200           0 :                      (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
     201             :     }
     202             : 
     203    29976960 :     dvelocity = dvelocity * _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph] +
     204    29976960 :                 (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity)) *
     205    29976960 :                     (_drelative_permeability_dvar[_qp][ph][pvar] / _fluid_viscosity[_qp][ph] -
     206    29976960 :                      _relative_permeability[_qp][ph] * _dfluid_viscosity_dvar[_qp][ph][pvar] /
     207             :                          std::pow(_fluid_viscosity[_qp][ph], 2)) *
     208             :                     _phi[_j][_qp];
     209             : 
     210    29976960 :     Real dvelocity_abs = 0.0;
     211    29976960 :     if (velocity_abs > 0.0)
     212    26824448 :       dvelocity_abs = velocity * dvelocity / velocity_abs;
     213             : 
     214             :     // Derivative of diffusion term (note: dispersivity is assumed constant)
     215    29976960 :     Real ddiffusion = _phi[_j][_qp] * _dporosity_qp_dvar[_qp][pvar] * _tortuosity[_qp][ph] *
     216    29976960 :                       _diffusion_coeff[_qp][ph][_fluid_component];
     217    29976960 :     ddiffusion += _phi[_j][_qp] * _porosity_qp[_qp] * _dtortuosity_dvar[_qp][ph][pvar] *
     218             :                   _diffusion_coeff[_qp][ph][_fluid_component];
     219    29976960 :     ddiffusion += _phi[_j][_qp] * _porosity_qp[_qp] * _tortuosity[_qp][ph] *
     220    29976960 :                   _ddiffusion_coeff_dvar[_qp][ph][_fluid_component][pvar];
     221    29976960 :     ddiffusion += _disp_trans[ph] * dvelocity_abs;
     222             : 
     223             :     // Derivative of dispersion term (note: dispersivity is assumed constant)
     224    29976960 :     RankTwoTensor ddispersion;
     225             :     ddispersion.zero();
     226    29976960 :     if (velocity_abs > 0.0)
     227             :     {
     228    26824448 :       RankTwoTensor dv2a, dv2b;
     229    26824448 :       dv2a = RankTwoTensor::outerProduct(velocity, dvelocity);
     230    26824448 :       dv2b = RankTwoTensor::outerProduct(dvelocity, velocity);
     231    26824448 :       ddispersion = (_disp_long[ph] - _disp_trans[ph]) * (dv2a + dv2b) / velocity_abs;
     232             :       ddispersion -=
     233    26824448 :           (_disp_long[ph] - _disp_trans[ph]) * v2 * dvelocity_abs / velocity_abs / velocity_abs;
     234             :     }
     235             : 
     236    29976960 :     dflux += _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] *
     237    29976960 :              (diffusion * _identity_tensor + dispersion) *
     238    29976960 :              _grad_mass_frac[_qp][ph][_fluid_component];
     239    29976960 :     dflux += _fluid_density_qp[_qp][ph] * (ddiffusion * _identity_tensor + ddispersion) *
     240    29976960 :              _grad_mass_frac[_qp][ph][_fluid_component];
     241             : 
     242             :     // NOTE: Here we assume that d(grad_mass_frac)/d(var) = d(mass_frac)/d(var) * grad_phi
     243             :     //       This is true for most PorousFlow scenarios, but not for chemical reactions
     244             :     //       where mass_frac is a nonlinear function of the primary MOOSE Variables
     245    29976960 :     dflux += _fluid_density_qp[_qp][ph] * (diffusion * _identity_tensor + dispersion) *
     246    29976960 :              _dmass_frac_dvar[_qp][ph][_fluid_component][pvar] * _grad_phi[_j][_qp];
     247             :   }
     248             : 
     249    28315520 :   return _grad_test[_i][_qp] * dflux;
     250             : }