Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
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 "LinearWCNSFV2PMomentumDriftFlux.h"
11 : #include "NS.h"
12 : #include "RhieChowMassFlux.h"
13 : #include "LinearFVBoundaryCondition.h"
14 : #include "LinearFVAdvectionDiffusionBC.h"
15 :
16 : registerMooseObject("NavierStokesApp", LinearWCNSFV2PMomentumDriftFlux);
17 :
18 : InputParameters
19 64 : LinearWCNSFV2PMomentumDriftFlux ::validParams()
20 : {
21 64 : auto params = LinearFVFluxKernel::validParams();
22 64 : params.addClassDescription("Implements the drift momentum flux source.");
23 128 : params.addRequiredParam<UserObjectName>(
24 : "rhie_chow_user_object",
25 : "The rhie-chow user-object which is used to determine the face velocity.");
26 128 : params.addRequiredParam<MooseFunctorName>("u_slip", "The slip velocity in the x direction.");
27 128 : params.addParam<MooseFunctorName>("v_slip", "The slip velocity in the y direction.");
28 128 : params.addParam<MooseFunctorName>("w_slip", "The slip velocity in the z direction.");
29 128 : params.addRequiredParam<MooseFunctorName>("rho_d", "Dispersed phase density.");
30 128 : params.addParam<MooseFunctorName>("fd", 0.0, "Fraction dispersed phase.");
31 128 : params.renameParam("fd", "fraction_dispersed", "");
32 :
33 128 : params.addParam<bool>(
34 128 : "force_boundary_execution", true, "This kernel should execute on boundaries by default");
35 128 : MooseEnum momentum_component("x=0 y=1 z=2");
36 128 : params.addRequiredParam<MooseEnum>(
37 : "momentum_component",
38 : momentum_component,
39 : "The component of the momentum equation that this kernel applies to.");
40 :
41 128 : MooseEnum coeff_interp_method("average harmonic", "harmonic");
42 128 : params.addParam<MooseEnum>("density_interp_method",
43 : coeff_interp_method,
44 : "Switch that can select face interpolation method for the density.");
45 :
46 64 : return params;
47 64 : }
48 :
49 32 : LinearWCNSFV2PMomentumDriftFlux ::LinearWCNSFV2PMomentumDriftFlux(const InputParameters & params)
50 : : LinearFVFluxKernel(params),
51 32 : _dim(_subproblem.mesh().dimension()),
52 32 : _mass_flux_provider(getUserObject<RhieChowMassFlux>("rhie_chow_user_object")),
53 64 : _rho_d(getFunctor<Real>("rho_d")),
54 64 : _f_d(getFunctor<Real>("fd")),
55 64 : _u_slip(getFunctor<Real>("u_slip")),
56 128 : _v_slip(isParamValid("v_slip") ? &getFunctor<Real>("v_slip") : nullptr),
57 64 : _w_slip(isParamValid("w_slip") ? &getFunctor<Real>("w_slip") : nullptr),
58 64 : _index(getParam<MooseEnum>("momentum_component")),
59 32 : _density_interp_method(
60 128 : Moose::FV::selectInterpolationMethod(getParam<MooseEnum>("density_interp_method")))
61 : {
62 32 : if (_dim >= 2 && !_v_slip)
63 0 : mooseError("In two or more dimensions, the v_slip velocity must be supplied using the 'v_slip' "
64 : "parameter");
65 32 : if (_dim >= 3 && !_w_slip)
66 0 : mooseError(
67 : "In three dimensions, the w_slip velocity must be supplied using the 'w_slip' parameter");
68 :
69 : // Note that since this is used in a segregated solver at this time, we don't need to declare
70 : // that this kernel may depend on a phase fraction variable
71 32 : }
72 :
73 : void
74 191760 : LinearWCNSFV2PMomentumDriftFlux::computeFlux()
75 : {
76 191760 : const auto & normal = _current_face_info->normal();
77 191760 : const auto state = determineState();
78 191760 : const bool on_boundary = Moose::FV::onBoundary(*this, *_current_face_info);
79 :
80 : Moose::FaceArg face_arg;
81 191760 : if (on_boundary)
82 57120 : face_arg = singleSidedFaceArg(_current_face_info);
83 : else
84 134640 : face_arg = makeCDFace(*_current_face_info);
85 :
86 : RealVectorValue u_slip_vel_vec;
87 191760 : if (_dim == 1)
88 0 : u_slip_vel_vec = RealVectorValue(_u_slip(face_arg, state), 0.0, 0.0);
89 191760 : else if (_dim == 2)
90 191760 : u_slip_vel_vec = RealVectorValue(_u_slip(face_arg, state), (*_v_slip)(face_arg, state), 0.0);
91 : else
92 0 : u_slip_vel_vec = RealVectorValue(
93 0 : _u_slip(face_arg, state), (*_v_slip)(face_arg, state), (*_w_slip)(face_arg, state));
94 :
95 : const auto uslipdotn = normal * u_slip_vel_vec;
96 :
97 : Real face_rho_fd;
98 191760 : if (on_boundary)
99 57120 : face_rho_fd = _rho_d(face_arg, state) * _f_d(face_arg, state);
100 : else
101 : {
102 134640 : const auto elem_arg = makeElemArg(_current_face_info->elemPtr());
103 269280 : const auto neigh_arg = makeElemArg(_current_face_info->neighborPtr());
104 :
105 134640 : Moose::FV::interpolate(_density_interp_method,
106 : face_rho_fd,
107 134640 : _rho_d(elem_arg, state) * _f_d(elem_arg, state),
108 269280 : _rho_d(neigh_arg, state) * _f_d(neigh_arg, state),
109 134640 : *_current_face_info,
110 : true);
111 : }
112 :
113 191760 : _face_flux = -face_rho_fd * uslipdotn * u_slip_vel_vec(_index);
114 191760 : }
115 :
116 : Real
117 134640 : LinearWCNSFV2PMomentumDriftFlux::computeElemMatrixContribution()
118 : {
119 134640 : const auto u_old = _var(makeCDFace(*_current_face_info), Moose::previousNonlinearState()).value();
120 134640 : if (std::abs(u_old) > 1e-6)
121 123320 : return _velocity_interp_coeffs.first * _face_flux / u_old * _current_face_area;
122 : else
123 : return 0.;
124 : }
125 :
126 : Real
127 134640 : LinearWCNSFV2PMomentumDriftFlux::computeNeighborMatrixContribution()
128 : {
129 134640 : const auto u_old = _var(makeCDFace(*_current_face_info), Moose::previousNonlinearState()).value();
130 134640 : if (std::abs(u_old) > 1e-6)
131 123320 : return _velocity_interp_coeffs.second * _face_flux / u_old * _current_face_area;
132 : else
133 : // place term on RHS if u_old is too close to 0
134 : return 0.;
135 : }
136 :
137 : Real
138 134640 : LinearWCNSFV2PMomentumDriftFlux::computeElemRightHandSideContribution()
139 : {
140 : // Get old velocity
141 134640 : const auto u_old = _var(makeCDFace(*_current_face_info), Moose::previousNonlinearState()).value();
142 134640 : const auto u = _var(makeCDFace(*_current_face_info), Moose::currentState()).value();
143 134640 : if (std::abs(u_old) > 1e-6)
144 123320 : return _velocity_interp_coeffs.first * _face_flux * (u / u_old - 1) * _current_face_area;
145 : else
146 11320 : return -_face_flux * _current_face_area;
147 : }
148 :
149 : Real
150 134640 : LinearWCNSFV2PMomentumDriftFlux::computeNeighborRightHandSideContribution()
151 : {
152 : // Get old velocity
153 134640 : const auto u_old = _var(makeCDFace(*_current_face_info), Moose::previousNonlinearState()).value();
154 134640 : const auto u = _var(makeCDFace(*_current_face_info), Moose::currentState()).value();
155 134640 : if (std::abs(u_old) > 1e-6)
156 123320 : return _velocity_interp_coeffs.second * _face_flux * (u / u_old - 1) * _current_face_area;
157 : else
158 11320 : return -_face_flux * _current_face_area;
159 : }
160 :
161 : Real
162 57120 : LinearWCNSFV2PMomentumDriftFlux::computeBoundaryRHSContribution(
163 : const LinearFVBoundaryCondition & /*bc*/)
164 : {
165 : // Lagging the whole term for now
166 : // TODO: make sure this only gets called once, and not once per BC
167 57120 : return -_face_flux * _current_face_area;
168 : }
169 :
170 : void
171 191760 : LinearWCNSFV2PMomentumDriftFlux::setupFaceData(const FaceInfo * face_info)
172 : {
173 191760 : LinearFVFluxKernel::setupFaceData(face_info);
174 :
175 : // Caching the mass flux on the face which will be reused in the advection term's matrix and right
176 : // hand side contributions
177 191760 : computeFlux();
178 :
179 : // Caching the interpolation coefficients so they will be reused for the matrix and right hand
180 : // side terms
181 : _velocity_interp_coeffs =
182 383520 : Moose::FV::interpCoeffs(_density_interp_method,
183 191760 : *_current_face_info,
184 : true,
185 191760 : _mass_flux_provider.getMassFlux(*_current_face_info));
186 191760 : }
|
