www.mooseframework.org
PorousFlowDispersiveFlux.C
Go to the documentation of this file.
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 
9 
10 // MOOSE includes
11 #include "MooseVariable.h"
12 
13 template <>
14 InputParameters
16 {
17  InputParameters params = validParams<Kernel>();
18  params.addParam<unsigned int>(
19  "fluid_component", 0, "The index corresponding to the fluid component for this kernel");
20  params.addRequiredParam<UserObjectName>(
21  "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
22  params.addRequiredParam<std::vector<Real>>(
23  "disp_long", "Vector of longitudinal dispersion coefficients for each phase");
24  params.addRequiredParam<std::vector<Real>>(
25  "disp_trans", "Vector of transverse dispersion coefficients for each phase");
26  params.addRequiredParam<RealVectorValue>("gravity",
27  "Gravitational acceleration vector downwards (m/s^2)");
28  params.addClassDescription(
29  "Dispersive and diffusive flux of the component given by fluid_component in all phases");
30  return params;
31 }
32 
33 PorousFlowDispersiveFlux::PorousFlowDispersiveFlux(const InputParameters & parameters)
34  : Kernel(parameters),
35 
36  _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
37  _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
38  "dPorousFlow_fluid_phase_density_qp_dvar")),
39  _grad_mass_frac(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
40  "PorousFlow_grad_mass_frac_qp")),
41  _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
42  "dPorousFlow_mass_frac_qp_dvar")),
43  _porosity_qp(getMaterialProperty<Real>("PorousFlow_porosity_qp")),
44  _dporosity_qp_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_qp_dvar")),
45  _tortuosity(getMaterialProperty<std::vector<Real>>("PorousFlow_tortuosity_qp")),
46  _dtortuosity_dvar(
47  getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_tortuosity_qp_dvar")),
48  _diffusion_coeff(
49  getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_diffusion_coeff_qp")),
50  _ddiffusion_coeff_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
51  "dPorousFlow_diffusion_coeff_qp_dvar")),
52  _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
53  _fluid_component(getParam<unsigned int>("fluid_component")),
54  _num_phases(_dictator.numPhases()),
55  _identity_tensor(RankTwoTensor::initIdentity),
56  _relative_permeability(
57  getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
58  _drelative_permeability_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
59  "dPorousFlow_relative_permeability_qp_dvar")),
60  _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
61  _dfluid_viscosity_dvar(
62  getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_qp_dvar")),
63  _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
64  _dpermeability_dvar(
65  getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
66  _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
67  "dPorousFlow_permeability_qp_dgradvar")),
68  _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
69  _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
70  "dPorousFlow_grad_porepressure_qp_dgradvar")),
71  _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
72  "dPorousFlow_grad_porepressure_qp_dvar")),
73  _gravity(getParam<RealVectorValue>("gravity")),
74  _disp_long(getParam<std::vector<Real>>("disp_long")),
75  _disp_trans(getParam<std::vector<Real>>("disp_trans"))
76 {
77  // Check that sufficient values of the dispersion coefficients have been entered
78  if (_disp_long.size() != _num_phases)
79  mooseError("The number of longitudinal dispersion coefficients disp_long in ",
80  _name,
81  " is not equal to the number of phases");
82 
83  if (_disp_trans.size() != _num_phases)
84  mooseError("The number of transverse dispersion coefficients disp_trans in ",
85  _name,
86  " is not equal to the number of phases");
87 }
88 
89 Real
91 {
92  RealVectorValue flux = 0.0;
93  RealVectorValue velocity;
94  Real velocity_abs;
95  RankTwoTensor v2;
96  RankTwoTensor dispersion;
97  dispersion.zero();
98  Real diffusion;
99 
100  for (unsigned int ph = 0; ph < _num_phases; ++ph)
101  {
102  // Diffusive component
103  diffusion =
104  _porosity_qp[_qp] * _tortuosity[_qp][ph] * _diffusion_coeff[_qp][ph][_fluid_component];
105 
106  // Calculate Darcy velocity
107  velocity = (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity) *
108  _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph]);
109  velocity_abs = std::sqrt(velocity * velocity);
110 
111  if (velocity_abs > 0.0)
112  {
113  v2.vectorOuterProduct(velocity, velocity);
114 
115  // Add longitudinal dispersion to diffusive component
116  diffusion += _disp_trans[ph] * velocity_abs;
117  dispersion = (_disp_long[ph] - _disp_trans[ph]) * v2 / velocity_abs;
118  }
119 
120  flux += _fluid_density_qp[_qp][ph] * (diffusion * _identity_tensor + dispersion) *
122  }
123  return _grad_test[_i][_qp] * flux;
124 }
125 
126 Real
128 {
129  return computeQpJac(_var.number());
130 }
131 
132 Real
134 {
135  return computeQpJac(jvar);
136 }
137 
138 Real
140 {
141  // If the variable is not a valid PorousFlow variable, set the Jacobian to 0
143  return 0.0;
144 
145  const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
146 
147  RealVectorValue velocity;
148  Real velocity_abs;
149  RankTwoTensor v2;
150  RankTwoTensor dispersion;
151  dispersion.zero();
152  Real diffusion;
153  RealVectorValue flux = 0.0;
154  RealVectorValue dflux = 0.0;
155 
156  for (unsigned int ph = 0; ph < _num_phases; ++ph)
157  {
158  // Diffusive component
159  diffusion =
160  _porosity_qp[_qp] * _tortuosity[_qp][ph] * _diffusion_coeff[_qp][ph][_fluid_component];
161 
162  // Calculate Darcy velocity
163  velocity = (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity) *
164  _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph]);
165  velocity_abs = std::sqrt(velocity * velocity);
166 
167  if (velocity_abs > 0.0)
168  {
169  v2.vectorOuterProduct(velocity, velocity);
170 
171  // Add longitudinal dispersion to diffusive component
172  diffusion += _disp_trans[ph] * velocity_abs;
173  dispersion = (_disp_long[ph] - _disp_trans[ph]) * v2 / velocity_abs;
174  }
175 
176  // Derivative of Darcy velocity
177  RealVectorValue dvelocity = _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
178  (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
179  for (unsigned i = 0; i < LIBMESH_DIM; ++i)
180  dvelocity += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
181  (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
182  dvelocity +=
183  _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
184  _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
185  dvelocity += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
186 
187  Real dvelocity_abs = 0.0;
188  if (velocity_abs > 0.0)
189  dvelocity_abs = velocity * dvelocity / velocity_abs;
190 
191  // Derivative of diffusion term (note: dispersivity is assumed constant)
192  Real ddiffusion = _phi[_j][_qp] * _dporosity_qp_dvar[_qp][pvar] * _tortuosity[_qp][ph] *
193  _diffusion_coeff[_qp][ph][_fluid_component];
194  ddiffusion += _phi[_j][_qp] * _porosity_qp[_qp] * _dtortuosity_dvar[_qp][ph][pvar] *
195  _diffusion_coeff[_qp][ph][_fluid_component];
196  ddiffusion += _phi[_j][_qp] * _porosity_qp[_qp] * _tortuosity[_qp][ph] *
198  ddiffusion += _disp_trans[ph] * dvelocity_abs;
199 
200  // Derivative of dispersion term (note: dispersivity is assumed constant)
201  RankTwoTensor ddispersion;
202  ddispersion.zero();
203  if (velocity_abs > 0.0)
204  {
205  RankTwoTensor dv2a, dv2b;
206  dv2a.vectorOuterProduct(velocity, dvelocity);
207  dv2b.vectorOuterProduct(dvelocity, velocity);
208  ddispersion = (_disp_long[ph] - _disp_trans[ph]) * (dv2a + dv2b) / velocity_abs;
209  ddispersion -=
210  (_disp_long[ph] - _disp_trans[ph]) * v2 * dvelocity_abs / velocity_abs / velocity_abs;
211  }
212 
213  dflux += _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] *
214  (diffusion * _identity_tensor + dispersion) *
216  dflux += _fluid_density_qp[_qp][ph] * (ddiffusion * _identity_tensor + ddispersion) *
218  dflux += _fluid_density_qp[_qp][ph] * (diffusion * _identity_tensor + dispersion) *
219  _dmass_frac_dvar[_qp][ph][_fluid_component][pvar] * _grad_phi[_j][_qp];
220  }
221 
222  return _grad_test[_i][_qp] * dflux;
223 }
InputParameters validParams< PorousFlowDispersiveFlux >()
const MaterialProperty< std::vector< std::vector< Real > > > & _dtortuosity_dvar
Derivative of tortuosity wrt PorousFlow variables.
const RankTwoTensor _identity_tensor
Identity tensor.
const MaterialProperty< std::vector< RealGradient > > & _grad_p
Gradient of the pore pressure in each phase.
virtual Real computeQpJacobian() override
const MaterialProperty< std::vector< Real > > & _tortuosity
Tortuosity tau_0 * tau_{alpha} for fluid phase alpha.
const MaterialProperty< std::vector< RealTensorValue > > & _dpermeability_dvar
Derivative of permeability wrt PorousFlow variables.
const MaterialProperty< std::vector< std::vector< RealGradient > > > & _dgrad_p_dvar
Derivative of Grad porepressure in each phase wrt PorousFlow variables.
const RealVectorValue _gravity
Gravitational acceleration.
const MaterialProperty< std::vector< Real > > & _fluid_density_qp
Fluid density for each phase (at the qp)
const MaterialProperty< std::vector< std::vector< Real > > > & _diffusion_coeff
Diffusion coefficients of component k in fluid phase alpha.
const MaterialProperty< std::vector< std::vector< RealGradient > > > & _grad_mass_frac
Gradient of mass fraction of each component in each phase.
const MaterialProperty< std::vector< Real > > & _fluid_viscosity
Viscosity of each component in each phase.
const MaterialProperty< std::vector< std::vector< RealTensorValue > > > & _dpermeability_dgradvar
d(permeabiity)/d(grad(porous-flow variable))
const unsigned int _num_phases
The number of fluid phases.
virtual Real computeQpResidual() override
const MaterialProperty< std::vector< Real > > & _dporosity_qp_dvar
Derivative of porosity wrt PorousFlow variables (at the qps)
const PorousFlowDictator & _dictator
PorousFlow Dictator UserObject.
PorousFlowDispersiveFlux(const InputParameters &parameters)
const MaterialProperty< RealTensorValue > & _permeability
Permeability of porous material.
const std::vector< Real > _disp_long
Longitudinal dispersivity for each phase.
Real computeQpJac(unsigned int jvar) const
Derivative of the residual with respect to the PorousFLow Variable with variable number jvar...
const MaterialProperty< std::vector< std::vector< Real > > > & _dfluid_density_qp_dvar
Derivative of the fluid density for each phase wrt PorousFlow variables (at the qp) ...
const MaterialProperty< std::vector< Real > > & _relative_permeability
Relative permeability of each phase.
This holds maps between the nonlinear variables used in a PorousFlow simulation and the variable numb...
const MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _ddiffusion_coeff_dvar
Derivative of the diffusion coefficients wrt PorousFlow variables.
virtual Real computeQpOffDiagJacobian(unsigned int jvar) override
const unsigned int _fluid_component
Index of the fluid component that this kernel acts on.
const MaterialProperty< std::vector< std::vector< std::vector< Real > > > > & _dmass_frac_dvar
Derivative of mass fraction wrt PorousFlow variables.
bool notPorousFlowVariable(unsigned int moose_var_num) const
returns true if moose_var_num is not a porous flow variabe
const MaterialProperty< std::vector< std::vector< Real > > > & _dgrad_p_dgrad_var
Derivative of Grad porepressure in each phase wrt grad(PorousFlow variables)
unsigned int porousFlowVariableNum(unsigned int moose_var_num) const
the PorousFlow variable number
const MaterialProperty< Real > & _porosity_qp
Porosity at the qps.
const std::vector< Real > _disp_trans
Transverse dispersivity for each phase.