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 "FVPorousFlowDispersiveFlux.h"
11 : #include "PorousFlowDictator.h"
12 : #include "MooseUtils.h"
13 :
14 : registerADMooseObject("PorousFlowApp", FVPorousFlowDispersiveFlux);
15 :
16 : InputParameters
17 246 : FVPorousFlowDispersiveFlux::validParams()
18 : {
19 246 : InputParameters params = FVFluxKernel::validParams();
20 246 : params += FVDiffusionInterpolationInterface::validParams();
21 : RealVectorValue g(0, 0, -9.81);
22 492 : params.addParam<RealVectorValue>("gravity", g, "Gravity vector. Defaults to (0, 0, -9.81)");
23 492 : params.addRequiredParam<UserObjectName>("PorousFlowDictator",
24 : "The PorousFlowDictator UserObject");
25 492 : params.addParam<unsigned int>("fluid_component", 0, "The fluid component");
26 492 : params.addRequiredParam<std::vector<Real>>(
27 : "disp_long", "Vector of longitudinal dispersion coefficients for each phase");
28 492 : params.addRequiredParam<std::vector<Real>>(
29 : "disp_trans", "Vector of transverse dispersion coefficients for each phase");
30 246 : params.addClassDescription(
31 : "Dispersive and diffusive flux of the component given by fluid_component in all phases");
32 246 : params.set<unsigned short>("ghost_layers") = 2;
33 :
34 : // We add the relationship manager here, this will select the right number of
35 : // ghosting layers depending on the chosen interpolation method
36 492 : params.addRelationshipManager(
37 : "ElementSideNeighborLayers",
38 : Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
39 : Moose::RelationshipManagerType::COUPLING,
40 : [](const InputParameters & obj_params, InputParameters & rm_params)
41 342 : { FVRelationshipManagerInterface::setRMParamsDiffusion(obj_params, rm_params, 3); });
42 :
43 246 : params.addClassDescription("Advective Darcy flux");
44 246 : return params;
45 0 : }
46 :
47 132 : FVPorousFlowDispersiveFlux::FVPorousFlowDispersiveFlux(const InputParameters & params)
48 : : FVFluxKernel(params),
49 : FVDiffusionInterpolationInterface(params),
50 132 : _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
51 132 : _num_phases(_dictator.numPhases()),
52 264 : _fluid_component(getParam<unsigned int>("fluid_component")),
53 264 : _density(getADMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
54 264 : _viscosity(getADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
55 132 : _viscosity_neighbor(
56 132 : getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
57 264 : _relperm(getADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
58 132 : _relperm_neighbor(
59 132 : getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
60 132 : _mass_fractions(
61 132 : getADMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
62 132 : _mass_fractions_neighbor(
63 132 : getNeighborADMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_qp")),
64 264 : _grad_mass_frac(getADMaterialProperty<std::vector<std::vector<RealGradient>>>(
65 : "PorousFlow_grad_mass_frac_qp")),
66 264 : _permeability(getADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
67 132 : _permeability_neighbor(
68 132 : getNeighborADMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
69 264 : _pressure(getADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
70 132 : _pressure_neighbor(
71 132 : getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")),
72 264 : _grad_p(getADMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
73 264 : _porosity(getADMaterialProperty<Real>("PorousFlow_porosity_qp")),
74 264 : _porosity_neighbor(getNeighborADMaterialProperty<Real>("PorousFlow_porosity_qp")),
75 264 : _tortuosity(getADMaterialProperty<std::vector<Real>>("PorousFlow_tortuosity_qp")),
76 132 : _tortuosity_neighbor(
77 132 : getNeighborADMaterialProperty<std::vector<Real>>("PorousFlow_tortuosity_qp")),
78 132 : _diffusion_coeff(
79 132 : getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_diffusion_coeff_qp")),
80 264 : _diffusion_coeff_neighbor(getNeighborMaterialProperty<std::vector<std::vector<Real>>>(
81 : "PorousFlow_diffusion_coeff_qp")),
82 264 : _gravity(getParam<RealVectorValue>("gravity")),
83 132 : _identity_tensor(ADRankTwoTensor::initIdentity),
84 264 : _disp_long(getParam<std::vector<Real>>("disp_long")),
85 396 : _disp_trans(getParam<std::vector<Real>>("disp_trans"))
86 : {
87 132 : if (_fluid_component >= _dictator.numComponents())
88 0 : paramError(
89 : "fluid_component",
90 : "The Dictator proclaims that the maximum fluid component index in this simulation is ",
91 0 : _dictator.numComponents() - 1,
92 : " whereas you have used ",
93 0 : _fluid_component,
94 : ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
95 :
96 : // Check that sufficient values of the dispersion coefficients have been entered
97 132 : if (_disp_long.size() != _num_phases)
98 0 : paramError(
99 : "disp_long",
100 : "The number of longitudinal dispersion coefficients is not equal to the number of phases");
101 :
102 132 : if (_disp_trans.size() != _num_phases)
103 0 : paramError("disp_trans",
104 : "The number of transverse dispersion coefficients disp_trans is not equal to the "
105 : "number of phases");
106 132 : }
107 :
108 : ADReal
109 194708 : FVPorousFlowDispersiveFlux::computeQpResidual()
110 : {
111 194708 : ADReal flux = 0.0;
112 194708 : ADRankTwoTensor dispersion;
113 : ADReal diffusion;
114 :
115 551616 : for (unsigned int p = 0; p < _num_phases; ++p)
116 : {
117 : // Diffusive component
118 : ADRealGradient gradX;
119 :
120 356908 : const auto state = determineState();
121 356908 : const auto dudn = gradUDotNormal(state, _correct_skewness);
122 :
123 : // If we are on a boundary face, use the reconstructed gradient computed in _grad_mass_frac
124 356908 : if (onBoundary(*_face_info))
125 : {
126 3096 : gradX = -1 * _grad_mass_frac[_qp][p][_fluid_component];
127 6192 : diffusion = _porosity[_qp] * _tortuosity[_qp][p] * _density[_qp][p] *
128 6192 : _diffusion_coeff[_qp][p][_fluid_component];
129 : }
130 : else
131 : {
132 : // If we are on an internal face, calculate the gradient explicitly
133 353812 : const auto & X_elem = _mass_fractions[_qp][p][_fluid_component];
134 353812 : const auto & X_neighbor = _mass_fractions_neighbor[_qp][p][_fluid_component];
135 :
136 707624 : gradX = (X_elem - X_neighbor) * _face_info->eCN() / _face_info->dCNMag();
137 :
138 : const auto coeff =
139 353812 : _porosity[_qp] * _tortuosity[_qp][p] * _diffusion_coeff[_qp][p][_fluid_component];
140 :
141 353812 : const auto coeff_neighbor = _porosity_neighbor[_qp] * _tortuosity_neighbor[_qp][p] *
142 353812 : _diffusion_coeff_neighbor[_qp][p][_fluid_component];
143 :
144 353812 : interpolate(
145 353812 : Moose::FV::InterpMethod::Average, diffusion, coeff, coeff_neighbor, *_face_info, true);
146 : }
147 :
148 : // Dispersive component. Calculate Darcy velocity
149 : ADRealGradient gradp;
150 : ADRealTensorValue mobility;
151 :
152 : // If we are on a boundary face, use the gradient computed in _grad_p
153 356908 : if (onBoundary(*_face_info))
154 : {
155 3096 : gradp = -_grad_p[_qp][p] + _density[_qp][p] * _gravity;
156 6192 : mobility = _relperm[_qp][p] * _permeability[_qp] * _density[_qp][p] / _viscosity[_qp][p];
157 : }
158 : else
159 : {
160 : // If we are on an internal face, calculate the gradient explicitly
161 353812 : const auto & p_elem = _pressure[_qp][p];
162 353812 : const auto & p_neighbor = _pressure_neighbor[_qp][p];
163 :
164 707624 : gradp = (p_elem - p_neighbor) * _face_info->eCN() / _face_info->dCNMag() +
165 353812 : _density[_qp][p] * _gravity;
166 :
167 353812 : const auto mobility_element = _relperm[_qp][p] * _permeability[_qp] / _viscosity[_qp][p];
168 :
169 : const auto mobility_neighbor =
170 353812 : _relperm_neighbor[_qp][p] * _permeability_neighbor[_qp] / _viscosity_neighbor[_qp][p];
171 :
172 353812 : interpolate(Moose::FV::InterpMethod::Upwind,
173 : mobility,
174 : mobility_element,
175 : mobility_neighbor,
176 : gradp,
177 353812 : *_face_info,
178 : true);
179 : }
180 :
181 356908 : const auto velocity = mobility * gradp;
182 356908 : const auto velocity_abs = MetaPhysicL::raw_value(velocity).norm();
183 :
184 356908 : if (!MooseUtils::isZero(velocity_abs))
185 : {
186 194284 : const auto v2 = ADRankTwoTensor::selfOuterProduct(velocity);
187 :
188 : // Add longitudinal dispersion to diffusive component
189 194284 : diffusion += _disp_trans[p] * velocity_abs;
190 388568 : dispersion = (_disp_long[p] - _disp_trans[p]) * v2 / velocity_abs;
191 : }
192 :
193 356908 : flux += (diffusion * _identity_tensor + dispersion) * _density[_qp][p] * gradX * _normal;
194 : }
195 :
196 194708 : return flux;
197 : }
|
