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 730 : BernoulliPressureVariable::validParams()
18 : {
19 730 : auto params = INSFVPressureVariable::validParams();
20 730 : params += ADFunctorInterface::validParams();
21 1460 : params.addRequiredParam<MooseFunctorName>("u", "The x-component of velocity");
22 1460 : params.addParam<MooseFunctorName>("v", 0, "The y-component of velocity");
23 1460 : params.addParam<MooseFunctorName>("w", 0, "The z-component of velocity");
24 730 : params.addRequiredParam<MooseFunctorName>(NS::porosity, "The porosity");
25 730 : params.addRequiredParam<MooseFunctorName>(NS::density, "The density");
26 1460 : params.addParam<std::vector<BoundaryName>>(
27 : "pressure_drop_sidesets", {}, "Sidesets over which form loss coefficients are to be applied");
28 1460 : 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 2190 : params.addParam<bool>(
33 : "allow_two_term_expansion_on_bernoulli_faces",
34 1460 : 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 1460 : params.addRelationshipManager(
44 : "ElementSideNeighborLayers",
45 : Moose::RelationshipManagerType::GEOMETRIC | Moose::RelationshipManagerType::ALGEBRAIC |
46 : Moose::RelationshipManagerType::COUPLING,
47 735 : [](const InputParameters & /*obj_params*/, InputParameters & rm_params)
48 735 : { rm_params.set<unsigned short>("layers") = 3; });
49 730 : return params;
50 0 : }
51 :
52 390 : BernoulliPressureVariable::BernoulliPressureVariable(const InputParameters & params)
53 : : INSFVPressureVariable(params),
54 : ADFunctorInterface(this),
55 390 : _u(nullptr),
56 390 : _v(nullptr),
57 390 : _w(nullptr),
58 390 : _eps(nullptr),
59 390 : _rho(nullptr),
60 390 : _pressure_drop_sidesets(getParam<std::vector<BoundaryName>>("pressure_drop_sidesets")),
61 390 : _pressure_drop_sideset_ids(this->_mesh.getBoundaryIDs(_pressure_drop_sidesets)),
62 780 : _pressure_drop_form_factors(getParam<std::vector<Real>>("pressure_drop_form_factors")),
63 390 : _allow_two_term_expansion_on_bernoulli_faces(
64 1170 : getParam<bool>("allow_two_term_expansion_on_bernoulli_faces"))
65 : {
66 390 : if (_allow_two_term_expansion_on_bernoulli_faces &&
67 22 : _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 390 : 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 390 : }
76 :
77 : void
78 390 : BernoulliPressureVariable::initialSetup()
79 : {
80 780 : _u = &getFunctor<ADReal>("u", _subproblem);
81 780 : _v = &getFunctor<ADReal>("v", _subproblem);
82 780 : _w = &getFunctor<ADReal>("w", _subproblem);
83 390 : _eps = &getFunctor<ADReal>(NS::porosity, _subproblem);
84 390 : _rho = &getFunctor<ADReal>(NS::density, _subproblem);
85 :
86 390 : INSFVPressureVariable::initialSetup();
87 390 : }
88 :
89 : std::pair<bool, ADRealVectorValue>
90 42570 : BernoulliPressureVariable::elemIsUpwind(const Elem & elem,
91 : const FaceInfo & fi,
92 : const Moose::StateArg & time) const
93 : {
94 42570 : const Moose::FaceArg face{
95 42570 : &fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr};
96 :
97 42570 : const VectorValue<ADReal> vel_face{(*_u)(face, time), (*_v)(face, time), (*_w)(face, time)};
98 42570 : const bool fi_elem_is_upwind = vel_face * fi.normal() > 0;
99 42570 : const bool elem_is_upwind = &elem == &fi.elem() ? fi_elem_is_upwind : !fi_elem_is_upwind;
100 :
101 42570 : return {elem_is_upwind, vel_face};
102 : }
103 :
104 : bool
105 993968 : BernoulliPressureVariable::isExtrapolatedBoundaryFace(const FaceInfo & fi,
106 : const Elem * const elem,
107 : const Moose::StateArg & time) const
108 : {
109 993968 : if (isDirichletBoundaryFace(fi, elem, time))
110 : return false;
111 979055 : if (!isInternalFace(fi))
112 : // We are neither a Dirichlet nor an internal face
113 : return true;
114 :
115 946527 : if (isSeparatorBoundary(fi))
116 : return true;
117 :
118 921232 : if (std::get<0>(NS::isPorosityJumpFace(*_eps, fi, time)))
119 : // We choose to extrapolate for the element that is downwind
120 8490 : return !std::get<0>(elemIsUpwind(*elem, fi, time));
121 :
122 : return false;
123 : }
124 :
125 : bool
126 1465459 : BernoulliPressureVariable::isDirichletBoundaryFace(const FaceInfo & fi,
127 : const Elem * const elem,
128 : const Moose::StateArg & time) const
129 : {
130 1465459 : if (INSFVPressureVariable::isDirichletBoundaryFace(fi, elem, time))
131 : return true;
132 :
133 1452675 : if (isInternalFace(fi))
134 : {
135 1419967 : if (isSeparatorBoundary(fi))
136 : return false;
137 :
138 1394672 : if (std::get<0>(NS::isPorosityJumpFace(*_eps, fi, time)))
139 : // We choose to apply the Dirichlet condition for the upwind element
140 25559 : return std::get<0>(elemIsUpwind(*elem, fi, time));
141 : }
142 :
143 : return false;
144 : }
145 :
146 : ADReal
147 14913 : 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 14913 : if (INSFVPressureVariable::isDirichletBoundaryFace(fi, elem, time))
154 6392 : return INSFVPressureVariable::getDirichletBoundaryFaceValue(fi, elem, time);
155 :
156 8521 : 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 8521 : const Moose::FaceArg face_elem{
165 8521 : &fi, Moose::FV::LimiterType::CentralDifference, true, false, &fi.elem(), nullptr};
166 8521 : const Moose::FaceArg face_neighbor{
167 8521 : &fi, Moose::FV::LimiterType::CentralDifference, true, false, fi.neighborPtr(), nullptr};
168 :
169 8521 : const auto [elem_is_upwind, vel_face] = elemIsUpwind(*elem, fi, time);
170 :
171 8521 : const bool fi_elem_is_upwind = elem == &fi.elem() ? elem_is_upwind : !elem_is_upwind;
172 8521 : const auto & downwind_face = fi_elem_is_upwind ? face_neighbor : face_elem;
173 :
174 8671 : const auto & upwind_elem = makeElemArg(fi_elem_is_upwind ? fi.elemPtr() : fi.neighborPtr());
175 17042 : 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 8521 : 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 8521 : const auto & rho_upwind = (*_rho)(upwind_elem, time);
191 8521 : const auto & rho_downwind = (*_rho)(downwind_elem, time);
192 :
193 17042 : const VectorValue<ADReal> interstitial_vel_upwind = vel_upwind * (1 / eps_upwind);
194 17042 : const VectorValue<ADReal> interstitial_vel_downwind = vel_downwind * (1 / eps_downwind);
195 :
196 8521 : const auto v_dot_n_upwind = interstitial_vel_upwind * fi.normal();
197 8521 : 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 8521 : ADReal factor = 0.0;
202 16606 : for (const auto & bd_id : fi.boundaryIDs())
203 46938 : for (const auto i : index_range(_pressure_drop_sideset_ids))
204 38853 : if (_pressure_drop_sideset_ids[i] == bd_id)
205 : factor += _pressure_drop_form_factors[i];
206 :
207 8521 : auto upwind_bernoulli_vel_chunk = 0.5 * rho_upwind * v_dot_n_upwind * v_dot_n_upwind;
208 17042 : 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 8521 : if (eps_downwind < eps_upwind)
213 : downwind_bernoulli_vel_chunk +=
214 8974 : 0.5 * factor * rho_downwind * v_dot_n_downwind * v_dot_n_downwind;
215 : else
216 8068 : upwind_bernoulli_vel_chunk += -0.5 * factor * rho_upwind * v_dot_n_upwind * v_dot_n_upwind;
217 :
218 : ADReal p_downwind;
219 8521 : if (_allow_two_term_expansion_on_bernoulli_faces)
220 27 : p_downwind = (*this)(downwind_face, time);
221 : else
222 8494 : p_downwind = (*this)(downwind_elem, time);
223 :
224 8521 : return p_downwind + downwind_bernoulli_vel_chunk - upwind_bernoulli_vel_chunk;
225 : }
|
