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 "WCNSFV2PMomentumAdvectionSlip.h"
11 : #include "NS.h"
12 : #include "FVUtils.h"
13 : #include "INSFVRhieChowInterpolator.h"
14 : #include "SystemBase.h"
15 : #include "NSFVUtils.h"
16 :
17 : registerMooseObject("NavierStokesApp", WCNSFV2PMomentumAdvectionSlip);
18 :
19 : InputParameters
20 164 : WCNSFV2PMomentumAdvectionSlip::validParams()
21 : {
22 164 : InputParameters params = INSFVMomentumAdvection::validParams();
23 164 : params.addClassDescription(
24 : "Computes the slip velocity advection kernel for two-phase mixture model.");
25 328 : params.addRequiredParam<MooseFunctorName>("u_slip", "The slip velocity in the x direction.");
26 328 : params.addParam<MooseFunctorName>("v_slip", "The slip velocity in the y direction.");
27 328 : params.addParam<MooseFunctorName>("w_slip", "The slip velocity in the z direction.");
28 328 : params.addRequiredParam<MooseFunctorName>("rho_d", "Dispersed phase density.");
29 328 : params.addParam<MooseFunctorName>("fd", 0.0, "Fraction dispersed phase.");
30 328 : params.renameParam("fd", "fraction_dispersed", "");
31 164 : params.setDocString(NS::density, "Main phase (not dispersed) density functor");
32 :
33 164 : return params;
34 0 : }
35 :
36 88 : WCNSFV2PMomentumAdvectionSlip::WCNSFV2PMomentumAdvectionSlip(const InputParameters & params)
37 : : INSFVMomentumAdvection(params),
38 88 : _rho_d(getFunctor<ADReal>("rho_d")),
39 88 : _dim(_subproblem.mesh().dimension()),
40 176 : _u_slip(getFunctor<ADReal>("u_slip")),
41 352 : _v_slip(isParamValid("v_slip") ? &getFunctor<ADReal>("v_slip") : nullptr),
42 176 : _w_slip(isParamValid("w_slip") ? &getFunctor<ADReal>("w_slip") : nullptr),
43 264 : _fd(getFunctor<ADReal>("fd"))
44 : {
45 88 : if (_dim >= 2 && !_v_slip)
46 0 : mooseError("In two or more dimensions, the v_slip velocity must be supplied using the 'v_slip' "
47 : "parameter");
48 88 : if (_dim >= 3 && !_w_slip)
49 0 : mooseError(
50 : "In three dimensions, the w_slip velocity must be supplied using the 'w_slip' parameter");
51 88 : }
52 :
53 : void
54 141984 : WCNSFV2PMomentumAdvectionSlip::computeResidualsAndAData(const FaceInfo & fi)
55 : {
56 : mooseAssert(!skipForBoundary(fi),
57 : "We shouldn't get in here if we're supposed to skip for a boundary");
58 :
59 : constexpr Real offset = 1e-15;
60 141984 : _face_info = &fi;
61 141984 : _normal = fi.normal();
62 283968 : _face_type = fi.faceType(std::make_pair(_var.number(), _var.sys().number()));
63 141984 : const auto state = determineState();
64 :
65 : using namespace Moose::FV;
66 :
67 141984 : const bool correct_skewness = _advected_interp_method == InterpMethod::SkewCorrectedAverage;
68 141984 : const bool on_boundary = onBoundary(fi);
69 :
70 141984 : _elem_residual = 0, _neighbor_residual = 0, _ae = 0, _an = 0;
71 :
72 : Moose::FaceArg face_arg;
73 141984 : if (on_boundary)
74 14688 : face_arg = singleSidedFaceArg();
75 : else
76 127296 : face_arg = Moose::FaceArg{
77 : &fi, Moose::FV::LimiterType::CentralDifference, true, false, nullptr, nullptr};
78 :
79 : ADRealVectorValue u_slip_vel_vec;
80 141984 : if (_dim == 1)
81 0 : u_slip_vel_vec(0) = _u_slip(face_arg, state);
82 141984 : if (_dim == 2)
83 283968 : u_slip_vel_vec = ADRealVectorValue(_u_slip(face_arg, state), (*_v_slip)(face_arg, state), 0.0);
84 141984 : if (_dim == 3)
85 0 : u_slip_vel_vec = ADRealVectorValue(
86 0 : _u_slip(face_arg, state), (*_v_slip)(face_arg, state), (*_w_slip)(face_arg, state));
87 :
88 141984 : const auto rho_mix = _rho(face_arg, state) +
89 141984 : (_rho_d(face_arg, state) - _rho(face_arg, state)) * _fd(face_arg, state);
90 :
91 141984 : const auto vdotn = _normal * u_slip_vel_vec;
92 :
93 141984 : if (on_boundary)
94 : {
95 14688 : const auto ssf = singleSidedFaceArg();
96 14688 : const Elem * const sided_elem = ssf.face_side;
97 14688 : const auto dof_number = sided_elem->dof_number(_sys.number(), _var.number(), 0);
98 14688 : const auto rho_face = rho_mix;
99 14688 : const auto eps_face = epsilon()(ssf, state);
100 14688 : const auto u_face = _var(ssf, state);
101 14688 : const Real d_u_face_d_dof = u_face.derivatives()[dof_number];
102 14688 : const auto coeff = vdotn * rho_face / eps_face;
103 :
104 14688 : if (sided_elem == fi.elemPtr())
105 : {
106 14688 : _ae = coeff * d_u_face_d_dof;
107 14688 : _elem_residual = coeff * u_face;
108 14688 : if (_approximate_as)
109 0 : _ae = _cs / fi.elem().n_sides() * rho_face / eps_face;
110 : }
111 : else
112 : {
113 0 : _an = -coeff * d_u_face_d_dof;
114 0 : _neighbor_residual = -coeff * u_face;
115 0 : if (_approximate_as)
116 0 : _an = _cs / fi.neighborPtr()->n_sides() * rho_face / eps_face;
117 : }
118 : }
119 : else // we are an internal fluid flow face
120 : {
121 127296 : const bool elem_is_upwind = MetaPhysicL::raw_value(u_slip_vel_vec) * _normal > 0;
122 127296 : const Moose::FaceArg advected_face_arg{&fi,
123 127296 : limiterType(_advected_interp_method),
124 : elem_is_upwind,
125 : correct_skewness,
126 : nullptr,
127 127296 : nullptr};
128 127296 : if (const auto [is_jump, eps_elem_face, eps_neighbor_face] =
129 127296 : NS::isPorosityJumpFace(epsilon(), fi, state);
130 127296 : is_jump)
131 : {
132 : // For a weakly compressible formulation, the density should not depend on pressure and
133 : // consequently the density should not be impacted by the pressure jump that occurs at a
134 : // porosity jump. Consequently we will allow evaluation of the density using both upstream and
135 : // downstream information
136 : const auto rho_face =
137 0 : _rho_d(advected_face_arg, state) - _rho(advected_face_arg, state) + offset;
138 :
139 : // We set the + and - sides of the superficial velocity equal to the interpolated value
140 0 : const auto & var_elem_face = u_slip_vel_vec(_index);
141 : const auto & var_neighbor_face = u_slip_vel_vec(_index);
142 :
143 0 : const auto elem_dof_number = fi.elem().dof_number(_sys.number(), _var.number(), 0);
144 0 : const auto neighbor_dof_number = fi.neighbor().dof_number(_sys.number(), _var.number(), 0);
145 :
146 0 : const auto d_var_elem_face_d_elem_dof = var_elem_face.derivatives()[elem_dof_number];
147 : const auto d_var_neighbor_face_d_neighbor_dof =
148 0 : var_neighbor_face.derivatives()[neighbor_dof_number];
149 :
150 : // We allow a discontintuity in the advective momentum flux at the jump face. Mainly the
151 : // advective flux is:
152 : // elem side:
153 : // rho * v_superficial / eps_elem * v_superficial = rho * v_interstitial_elem * v_superficial
154 : // neighbor side:
155 : // rho * v_superficial / eps_neigh * v_superficial = rho * v_interstitial_neigh *
156 : // v_superficial
157 0 : const auto elem_coeff = vdotn * rho_face / eps_elem_face;
158 0 : const auto neighbor_coeff = vdotn * rho_face / eps_neighbor_face;
159 0 : _ae = elem_coeff * d_var_elem_face_d_elem_dof;
160 0 : _elem_residual = elem_coeff * var_elem_face;
161 0 : _an = -neighbor_coeff * d_var_neighbor_face_d_neighbor_dof;
162 0 : _neighbor_residual = -neighbor_coeff * var_neighbor_face;
163 0 : if (_approximate_as)
164 : {
165 0 : _ae = _cs / fi.elem().n_sides() * rho_face / eps_elem_face;
166 0 : _an = _cs / fi.neighborPtr()->n_sides() * rho_face / eps_elem_face;
167 : }
168 : }
169 : else
170 : {
171 : const auto [interp_coeffs, advected] =
172 127296 : interpCoeffsAndAdvected(*_rho_u, advected_face_arg, state);
173 :
174 127296 : const auto elem_arg = elemArg();
175 127296 : const auto neighbor_arg = neighborArg();
176 :
177 127296 : const auto rho_elem = _rho_d(elem_arg, state) - _rho(elem_arg, state) + offset;
178 127296 : const auto rho_neighbor = _rho_d(neighbor_arg, state) - _rho(neighbor_arg, state) + offset;
179 127296 : const auto eps_elem = epsilon()(elem_arg, state),
180 127296 : eps_neighbor = epsilon()(neighbor_arg, state);
181 127296 : const auto var_elem = advected.first / rho_elem * eps_elem,
182 127296 : var_neighbor = advected.second / rho_neighbor * eps_neighbor;
183 :
184 127296 : _ae = vdotn * rho_elem / eps_elem * interp_coeffs.first;
185 : // Minus sign because we apply a minus sign to the residual in computeResidual
186 254592 : _an = -vdotn * rho_neighbor / eps_neighbor * interp_coeffs.second;
187 :
188 127296 : _elem_residual = _ae * var_elem - _an * var_neighbor;
189 127296 : _neighbor_residual = -_elem_residual;
190 :
191 127296 : if (_approximate_as)
192 : {
193 0 : _ae = _cs / fi.elem().n_sides() * rho_elem / eps_elem;
194 0 : _an = _cs / fi.neighborPtr()->n_sides() * rho_neighbor / eps_neighbor;
195 : }
196 : }
197 : }
198 :
199 141984 : _ae *= fi.faceArea() * fi.faceCoord();
200 141984 : _an *= fi.faceArea() * fi.faceCoord();
201 141984 : _elem_residual *= fi.faceArea() * fi.faceCoord();
202 141984 : _neighbor_residual *= fi.faceArea() * fi.faceCoord();
203 141984 : }
|
