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 "LinearFVTKESourceSink.h"
11 : #include "Assembly.h"
12 : #include "SubProblem.h"
13 : #include "NavierStokesMethods.h"
14 :
15 : registerMooseObject("NavierStokesApp", LinearFVTKESourceSink);
16 :
17 : InputParameters
18 142 : LinearFVTKESourceSink::validParams()
19 : {
20 142 : InputParameters params = LinearFVElementalKernel::validParams();
21 142 : params.addClassDescription("Elemental kernel to compute the production and destruction "
22 : " terms of turbulent kinetic energy (TKE).");
23 284 : params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
24 284 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
25 284 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
26 142 : params.addRequiredParam<MooseFunctorName>(NS::TKED,
27 : "Coupled turbulent kinetic energy dissipation rate.");
28 142 : params.addRequiredParam<MooseFunctorName>(NS::density, "Fluid density");
29 142 : params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
30 142 : params.addRequiredParam<MooseFunctorName>(NS::mu_t, "Turbulent viscosity.");
31 :
32 284 : params.addParam<std::vector<BoundaryName>>(
33 : "walls", {}, "Boundaries that correspond to solid walls.");
34 284 : MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
35 :
36 284 : params.addParam<MooseEnum>("wall_treatment",
37 : wall_treatment,
38 : "The method used for computing the wall functions "
39 : "'eq_newton', 'eq_incremental', 'eq_linearized', 'neq'");
40 284 : params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
41 284 : params.addParam<Real>("C_pl", 10.0, "Production Limiter Constant Multiplier.");
42 :
43 142 : return params;
44 142 : }
45 :
46 71 : LinearFVTKESourceSink::LinearFVTKESourceSink(const InputParameters & params)
47 : : LinearFVElementalKernel(params),
48 71 : _dim(_subproblem.mesh().dimension()),
49 142 : _u_var(getFunctor<Real>("u")),
50 213 : _v_var(params.isParamValid("v") ? &(getFunctor<Real>("v")) : nullptr),
51 71 : _w_var(params.isParamValid("w") ? &(getFunctor<Real>("w")) : nullptr),
52 71 : _epsilon(getFunctor<Real>(NS::TKED)),
53 71 : _rho(getFunctor<Real>(NS::density)),
54 71 : _mu(getFunctor<Real>(NS::mu)),
55 71 : _mu_t(getFunctor<Real>(NS::mu_t)),
56 142 : _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
57 142 : _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
58 142 : _C_mu(getParam<Real>("C_mu")),
59 213 : _C_pl(getParam<Real>("C_pl"))
60 : {
61 71 : if (_dim >= 2 && !_v_var)
62 0 : paramError("v", "In two or more dimensions, the v velocity must be supplied!");
63 :
64 71 : if (_dim >= 3 && !_w_var)
65 0 : paramError("w", "In three or more dimensions, the w velocity must be supplied!");
66 :
67 : // // Strain tensor term requires velocity gradients;
68 71 : if (dynamic_cast<const MooseLinearVariableFV<Real> *>(&_u_var))
69 0 : requestVariableCellGradient(getParam<MooseFunctorName>("u"));
70 71 : if (dynamic_cast<const MooseLinearVariableFV<Real> *>(_v_var))
71 0 : requestVariableCellGradient(getParam<MooseFunctorName>("v"));
72 71 : if (dynamic_cast<const MooseLinearVariableFV<Real> *>(_w_var))
73 0 : requestVariableCellGradient(getParam<MooseFunctorName>("w"));
74 71 : }
75 :
76 : void
77 373 : LinearFVTKESourceSink::initialSetup()
78 : {
79 373 : LinearFVElementalKernel::initialSetup();
80 746 : NS::getWallBoundedElements(
81 373 : _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
82 373 : NS::getWallDistance(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _dist);
83 373 : NS::getElementFaceArgs(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _face_infos);
84 373 : }
85 :
86 : Real
87 633792 : LinearFVTKESourceSink::computeMatrixContribution()
88 : {
89 : /*
90 : Matrix contribution:
91 : - Computes near-wall TKE destruction
92 : - Computes bulk TKE destruction
93 : */
94 :
95 : // Useful variables
96 633792 : const auto state = determineState();
97 633792 : const auto elem_arg = makeElemArg(_current_elem_info->elem());
98 633792 : const Real rho = _rho(elem_arg, state);
99 :
100 633792 : if (_wall_bounded.find(_current_elem_info->elem()) != _wall_bounded.end())
101 : {
102 :
103 : // Useful variables
104 164784 : const Real mu = _mu(elem_arg, state);
105 164784 : const Real TKE = _var.getElemValue(*_current_elem_info, state);
106 :
107 : // Assembly variables
108 : Real production = 0.0;
109 : Real destruction = 0.0;
110 : Real tot_weight = 0.0;
111 : std::vector<Real> y_plus_vec, velocity_grad_norm_vec;
112 :
113 : // Get near-wall element velocity vector
114 164784 : RealVectorValue velocity(_u_var(elem_arg, state));
115 164784 : if (_v_var)
116 164784 : velocity(1) = (*_v_var)(elem_arg, state);
117 164784 : if (_w_var)
118 0 : velocity(2) = (*_w_var)(elem_arg, state);
119 :
120 : // Get faceInfo and distance vector for between all walls and elements
121 164784 : const auto & face_info_vec = libmesh_map_find(_face_infos, _current_elem_info->elem());
122 164784 : const auto & distance_vec = libmesh_map_find(_dist, _current_elem_info->elem());
123 : mooseAssert(distance_vec.size(), "Should have found a distance vector");
124 : mooseAssert(distance_vec.size() == face_info_vec.size(),
125 : "Should be as many distance vectors as face info vectors");
126 :
127 : // Populate y+, near wall velocity gradient, and update weights
128 330208 : for (unsigned int i = 0; i < distance_vec.size(); i++)
129 : {
130 165424 : const auto parallel_speed = NS::computeSpeed<Real>(
131 165424 : velocity - velocity * face_info_vec[i]->normal() * face_info_vec[i]->normal());
132 165424 : const auto distance = distance_vec[i];
133 :
134 : Real y_plus;
135 165424 : if (_wall_treatment == NS::WallTreatmentEnum::NEQ) // Non-equilibrium --> Non-iterative
136 30000 : y_plus = distance * std::sqrt(std::sqrt(_C_mu) * TKE) * rho / mu;
137 : else // Equilibrium --> Iterative
138 135424 : y_plus = NS::findyPlus<Real>(mu, rho, std::max(parallel_speed, 1e-10), distance);
139 :
140 165424 : y_plus_vec.push_back(y_plus);
141 165424 : const Real velocity_grad_norm = parallel_speed / distance;
142 165424 : velocity_grad_norm_vec.push_back(velocity_grad_norm);
143 165424 : tot_weight += 1.0;
144 : }
145 :
146 : // Compute near-wall production and destruction
147 330208 : for (unsigned int i = 0; i < y_plus_vec.size(); i++)
148 : {
149 165424 : const auto y_plus = y_plus_vec[i];
150 :
151 165424 : const auto fi = face_info_vec[i];
152 165424 : const bool defined_on_elem_side = _var.hasFaceSide(*fi, true);
153 165424 : const Elem * const loc_elem = defined_on_elem_side ? &fi->elem() : fi->neighborPtr();
154 165424 : const Moose::FaceArg facearg = {
155 165424 : fi, Moose::FV::LimiterType::CentralDifference, false, false, loc_elem, nullptr};
156 :
157 165424 : const Real wall_mut = _mu_t(facearg, state);
158 165424 : const Real wall_mu = _mu(facearg, state);
159 165424 : const auto tau_w = (wall_mut + wall_mu) * velocity_grad_norm_vec[i];
160 165424 : const auto destruction_visc = 2.0 * wall_mu / Utility::pow<2>(distance_vec[i]) / tot_weight;
161 165424 : const auto destruction_log = std::pow(_C_mu, 0.75) * rho * std::pow(TKE, 0.5) /
162 165424 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
163 :
164 165424 : if (y_plus < 11.25)
165 19477 : destruction += destruction_visc;
166 : else
167 : {
168 145947 : destruction += destruction_log;
169 145947 : production += tau_w * std::pow(_C_mu, 0.25) / std::sqrt(TKE) /
170 145947 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
171 : }
172 : }
173 :
174 : // Assign terms to matrix to solve implicitly (they get multiplied by TKE)
175 164784 : return (destruction - production) * _current_elem_volume;
176 164784 : }
177 : else
178 : {
179 : // Compute destruction
180 : const auto destruction =
181 469008 : rho * _epsilon(elem_arg, state) / _var.getElemValue(*_current_elem_info, state);
182 :
183 : // Assign to matrix to solve implicitly
184 469008 : return destruction * _current_elem_volume;
185 : }
186 : }
187 :
188 : Real
189 633792 : LinearFVTKESourceSink::computeRightHandSideContribution()
190 : {
191 :
192 : // Useful variables
193 633792 : const auto state = determineState();
194 633792 : const auto elem_arg = makeElemArg(_current_elem_info->elem());
195 :
196 633792 : if (_wall_bounded.find(_current_elem_info->elem()) != _wall_bounded.end()) // Wall bounded
197 : return 0.0; // Do nothing
198 : else // Not wall bounded
199 : {
200 : // Compute TKE production
201 469008 : const auto subdomain_id = _current_elem_info->elem()->subdomain_id();
202 469008 : const auto coord_sys = _subproblem.getCoordSystem(subdomain_id);
203 : const unsigned int rz_radial_coord =
204 469008 : coord_sys == Moose::COORD_RZ ? _subproblem.getAxisymmetricRadialCoord() : 0;
205 469008 : const auto symmetric_strain_tensor_sq_norm = NS::computeShearStrainRateNormSquared<Real>(
206 : _u_var, _v_var, _w_var, elem_arg, state, coord_sys, rz_radial_coord);
207 :
208 469008 : auto production = _mu_t(elem_arg, state) * symmetric_strain_tensor_sq_norm;
209 :
210 : // k-Production limiter (needed for flows with stagnation zones)
211 469008 : const Real production_limit = _C_pl * _rho(elem_arg, state) * _epsilon(elem_arg, state);
212 :
213 : // Apply production limiter
214 469008 : production = std::min(production, production_limit);
215 :
216 : // Assign production to RHS
217 469008 : return production * _current_elem_volume;
218 : }
219 : }
|
