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 "BernoulliPressureVariable.h"
11 : #include "NS.h"
12 : #include "NSFVUtils.h"
13 :
14 : registerMooseObject("NavierStokesApp", BernoulliPressureVariable);
15 :
16 : InputParameters
17 360 : BernoulliPressureVariable::validParams()
18 : {
19 360 : auto params = INSFVPressureVariable::validParams();
20 360 : params += ADFunctorInterface::validParams();
21 720 : params.addRequiredParam<MooseFunctorName>("u", "The x-component of velocity");
22 720 : params.addParam<MooseFunctorName>("v", 0, "The y-component of velocity");
23 720 : params.addParam<MooseFunctorName>("w", 0, "The z-component of velocity");
24 360 : params.addRequiredParam<MooseFunctorName>(NS::porosity, "The porosity");
25 360 : params.addRequiredParam<MooseFunctorName>(NS::density, "The density");
26 720 : params.addParam<std::vector<BoundaryName>>(
27 : "pressure_drop_sidesets", {}, "Sidesets over which form loss coefficients are to be applied");
28 720 : params.addParam<std::vector<Real>>(
29 : "pressure_drop_form_factors",
30 : {},
31 : "User-supplied form loss coefficients to be applied over the sidesets listed above");
32 1080 : params.addParam<bool>(
33 : "allow_two_term_expansion_on_bernoulli_faces",
34 720 : false,
35 : "Switch to enable the two-term extrapolation on porosity jump faces. "
36 : "WARNING: This might lead to crushes in parallel runs if porosity jump faces are connected "
37 : "with one cell (usually corners) due to the insufficient number of ghosted "
38 : "layers.");
39 :
40 : // Assuming that this variable is used for advection problems, due to the
41 : // utilization of the pressure gradient in the advecting velocity
42 : // through the Rhie-Chow interpolation, we have to extend the ghosted layers.
43 720 : params.addRelationshipManager(
44 : "ElementSideNeighborLayers",
45 : Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
46 : Moose::RelationshipManagerType::COUPLING,
47 375 : [](const InputParameters & /*obj_params*/, InputParameters & rm_params)
48 375 : { rm_params.set<unsigned short>("layers") = 3; });
49 360 : return params;
50 0 : }
51 :
52 190 : BernoulliPressureVariable::BernoulliPressureVariable(const InputParameters & params)
53 : : INSFVPressureVariable(params),
54 : ADFunctorInterface(this),
55 190 : _u(nullptr),
56 190 : _v(nullptr),
57 190 : _w(nullptr),
58 190 : _eps(nullptr),
59 190 : _rho(nullptr),
60 190 : _pressure_drop_sidesets(getParam<std::vector<BoundaryName>>("pressure_drop_sidesets")),
61 190 : _pressure_drop_sideset_ids(this->_mesh.getBoundaryIDs(_pressure_drop_sidesets)),
62 380 : _pressure_drop_form_factors(getParam<std::vector<Real>>("pressure_drop_form_factors")),
63 190 : _allow_two_term_expansion_on_bernoulli_faces(
64 570 : getParam<bool>("allow_two_term_expansion_on_bernoulli_faces"))
65 : {
66 190 : if (_allow_two_term_expansion_on_bernoulli_faces &&
67 10 : _face_interp_method == Moose::FV::InterpMethod::SkewCorrectedAverage)
68 0 : paramError(
69 : "allow_two_term_expansion_on_bernoulli_faces",
70 : "Skewness correction with two term extrapolation on Bernoulli faces is not supported!");
71 :
72 190 : if (_pressure_drop_sidesets.size() != _pressure_drop_form_factors.size())
73 0 : paramError("pressure_drop_sidesets",
74 : "The number of sidesets and the number of supplied form losses are not equal!");
75 190 : }
76 :
77 : void
78 190 : BernoulliPressureVariable::initialSetup()
79 : {
80 380 : _u = &getFunctor<ADReal>("u", _subproblem);
81 380 : _v = &getFunctor<ADReal>("v", _subproblem);
82 380 : _w = &getFunctor<ADReal>("w", _subproblem);
83 190 : _eps = &getFunctor<ADReal>(NS::porosity, _subproblem);
84 190 : _rho = &getFunctor<ADReal>(NS::density, _subproblem);
85 :
86 190 : INSFVPressureVariable::initialSetup();
87 190 : }
88 :
89 : std::pair<bool, ADRealVectorValue>
90 28731 : BernoulliPressureVariable::elemIsUpwind(const Elem & elem,
91 : const FaceInfo & fi,
92 : const Moose::StateArg & time) const
93 : {
94 28731 : const Moose::FaceArg face{
95 28731 : &fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr};
96 :
97 28731 : const VectorValue<ADReal> vel_face{(*_u)(face, time), (*_v)(face, time), (*_w)(face, time)};
98 28731 : const bool fi_elem_is_upwind = vel_face * fi.normal() > 0;
99 28731 : const bool elem_is_upwind = &elem == &fi.elem() ? fi_elem_is_upwind : !fi_elem_is_upwind;
100 :
101 28731 : return {elem_is_upwind, vel_face};
102 : }
103 :
104 : bool
105 694903 : BernoulliPressureVariable::isExtrapolatedBoundaryFace(const FaceInfo & fi,
106 : const Elem * const elem,
107 : const Moose::StateArg & time) const
108 : {
109 694903 : if (isDirichletBoundaryFace(fi, elem, time))
110 : return false;
111 684622 : if (!isInternalFace(fi))
112 : // We are neither a Dirichlet nor an internal face
113 : return true;
114 :
115 661670 : if (isSeparatorBoundary(fi))
116 : return true;
117 :
118 643474 : if (std::get<0>(NS::isPorosityJumpFace(*_eps, fi, time)))
119 : // We choose to extrapolate for the element that is downwind
120 5742 : return !std::get<0>(elemIsUpwind(*elem, fi, time));
121 :
122 : return false;
123 : }
124 :
125 : bool
126 1024200 : BernoulliPressureVariable::isDirichletBoundaryFace(const FaceInfo & fi,
127 : const Elem * const elem,
128 : const Moose::StateArg & time) const
129 : {
130 1024200 : if (INSFVPressureVariable::isDirichletBoundaryFace(fi, elem, time))
131 : return true;
132 :
133 1015124 : if (isInternalFace(fi))
134 : {
135 992040 : if (isSeparatorBoundary(fi))
136 : return false;
137 :
138 973844 : if (std::get<0>(NS::isPorosityJumpFace(*_eps, fi, time)))
139 : // We choose to apply the Dirichlet condition for the upwind element
140 17246 : return std::get<0>(elemIsUpwind(*elem, fi, time));
141 : }
142 :
143 : return false;
144 : }
145 :
146 : ADReal
147 10281 : BernoulliPressureVariable::getDirichletBoundaryFaceValue(const FaceInfo & fi,
148 : const Elem * const elem,
149 : const Moose::StateArg & time) const
150 : {
151 : mooseAssert(isDirichletBoundaryFace(fi, elem, time), "This better be a Dirichlet face");
152 :
153 10281 : if (INSFVPressureVariable::isDirichletBoundaryFace(fi, elem, time))
154 4538 : return INSFVPressureVariable::getDirichletBoundaryFaceValue(fi, elem, time);
155 :
156 5743 : const auto [is_jump_face, eps_elem, eps_neighbor] = NS::isPorosityJumpFace(*_eps, fi, time);
157 : #ifndef NDEBUG
158 : mooseAssert(is_jump_face,
159 : "If we are not a traditional Dirichlet face, then we must be a jump face");
160 : #else
161 : libmesh_ignore(is_jump_face);
162 : #endif
163 :
164 5743 : const Moose::FaceArg face_elem{
165 5743 : &fi, Moose::FV::LimiterType::CentralDifference, true, false, &fi.elem(), nullptr};
166 5743 : const Moose::FaceArg face_neighbor{
167 5743 : &fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.neighborPtr(), nullptr};
168 :
169 5743 : const auto [elem_is_upwind, vel_face] = elemIsUpwind(*elem, fi, time);
170 :
171 5743 : const bool fi_elem_is_upwind = elem == &fi.elem() ? elem_is_upwind : !elem_is_upwind;
172 5743 : const auto & downwind_face = fi_elem_is_upwind ? face_neighbor : face_elem;
173 :
174 5793 : const auto & upwind_elem = makeElemArg(fi_elem_is_upwind ? fi.elemPtr() : fi.neighborPtr());
175 11486 : const auto & downwind_elem = makeElemArg(fi_elem_is_upwind ? fi.neighborPtr() : fi.elemPtr());
176 :
177 : const auto & vel_upwind = vel_face;
178 : const auto & vel_downwind = vel_face;
179 :
180 5743 : const auto & eps_upwind = fi_elem_is_upwind ? eps_elem : eps_neighbor;
181 : const auto & eps_downwind = fi_elem_is_upwind ? eps_neighbor : eps_elem;
182 :
183 : /*
184 : If a weakly compressible material is used where the density slightly
185 : depends on pressure. Given that the two-term expansion on Bernoulli faces is
186 : enabled, we take use the downwind face assuming the continuity (minimal jump) of the
187 : density. This protects against an infinite recursion between pressure and
188 : density.
189 : */
190 5743 : const auto & rho_upwind = (*_rho)(upwind_elem, time);
191 5743 : const auto & rho_downwind = (*_rho)(downwind_elem, time);
192 :
193 11486 : const VectorValue<ADReal> interstitial_vel_upwind = vel_upwind * (1 / eps_upwind);
194 11486 : const VectorValue<ADReal> interstitial_vel_downwind = vel_downwind * (1 / eps_downwind);
195 :
196 5743 : const auto v_dot_n_upwind = interstitial_vel_upwind * fi.normal();
197 5743 : const auto v_dot_n_downwind = interstitial_vel_downwind * fi.normal();
198 :
199 : // Iterate through sidesets to see if we have associated irreversible pressure losses
200 : // or not.
201 5743 : ADReal factor = 0.0;
202 11184 : for (const auto & bd_id : fi.boundaryIDs())
203 31486 : for (const auto i : index_range(_pressure_drop_sideset_ids))
204 26045 : if (_pressure_drop_sideset_ids[i] == bd_id)
205 : factor += _pressure_drop_form_factors[i];
206 :
207 5743 : auto upwind_bernoulli_vel_chunk = 0.5 * rho_upwind * v_dot_n_upwind * v_dot_n_upwind;
208 11486 : auto downwind_bernoulli_vel_chunk = 0.5 * rho_downwind * v_dot_n_downwind * v_dot_n_downwind;
209 :
210 : // We add the additional, irreversible pressure loss here.
211 : // If it is a contraction we have to use the downwind values, otherwise the upwind values.
212 5743 : if (eps_downwind < eps_upwind)
213 : downwind_bernoulli_vel_chunk +=
214 6166 : 0.5 * factor * rho_downwind * v_dot_n_downwind * v_dot_n_downwind;
215 : else
216 5320 : upwind_bernoulli_vel_chunk += -0.5 * factor * rho_upwind * v_dot_n_upwind * v_dot_n_upwind;
217 :
218 : ADReal p_downwind;
219 5743 : if (_allow_two_term_expansion_on_bernoulli_faces)
220 18 : p_downwind = (*this)(downwind_face, time);
221 : else
222 5725 : p_downwind = (*this)(downwind_elem, time);
223 :
224 5743 : return p_downwind + downwind_bernoulli_vel_chunk - upwind_bernoulli_vel_chunk;
225 : }
|
