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 "LinearFVTKEDSourceSink.h"
11 : #include "Assembly.h"
12 : #include "SubProblem.h"
13 : #include "NavierStokesMethods.h"
14 :
15 : registerMooseObject("NavierStokesApp", LinearFVTKEDSourceSink);
16 :
17 : InputParameters
18 142 : LinearFVTKEDSourceSink::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 dissipation (TKED).");
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::TKE, "Coupled turbulent kinetic energy.");
27 142 : params.addRequiredParam<MooseFunctorName>(NS::density, "fluid density");
28 142 : params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
29 142 : params.addRequiredParam<MooseFunctorName>(NS::mu_t, "Turbulent viscosity.");
30 :
31 284 : params.addParam<std::vector<BoundaryName>>(
32 : "walls", {}, "Boundaries that correspond to solid walls.");
33 284 : MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
34 284 : params.addParam<MooseEnum>("wall_treatment",
35 : wall_treatment,
36 : "The method used for computing the wall functions "
37 : "'eq_newton', 'eq_incremental', 'eq_linearized', 'neq'");
38 :
39 284 : params.addParam<MooseFunctorName>("C1_eps", 1.44, "First epsilon coefficient");
40 284 : params.addParam<MooseFunctorName>("C2_eps", 1.92, "Second epsilon coefficient");
41 284 : params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
42 284 : params.addParam<Real>("C_pl", 10.0, "Production limiter constant multiplier.");
43 :
44 142 : return params;
45 142 : }
46 :
47 71 : LinearFVTKEDSourceSink::LinearFVTKEDSourceSink(const InputParameters & params)
48 : : LinearFVElementalKernel(params),
49 71 : _dim(_subproblem.mesh().dimension()),
50 142 : _u_var(getFunctor<Real>("u")),
51 213 : _v_var(params.isParamValid("v") ? &(getFunctor<Real>("v")) : nullptr),
52 71 : _w_var(params.isParamValid("w") ? &(getFunctor<Real>("w")) : nullptr),
53 71 : _k(getFunctor<Real>(NS::TKE)),
54 71 : _rho(getFunctor<Real>(NS::density)),
55 71 : _mu(getFunctor<Real>(NS::mu)),
56 71 : _mu_t(getFunctor<Real>(NS::mu_t)),
57 142 : _wall_boundary_names(getParam<std::vector<BoundaryName>>("walls")),
58 142 : _wall_treatment(getParam<MooseEnum>("wall_treatment").getEnum<NS::WallTreatmentEnum>()),
59 142 : _C1_eps(getFunctor<Real>("C1_eps")),
60 142 : _C2_eps(getFunctor<Real>("C2_eps")),
61 142 : _C_mu(getParam<Real>("C_mu")),
62 213 : _C_pl(getParam<Real>("C_pl"))
63 : {
64 71 : if (_dim >= 2 && !_v_var)
65 0 : paramError("v", "In two or more dimensions, the v velocity must be supplied!");
66 :
67 71 : if (_dim >= 3 && !_w_var)
68 0 : paramError("w", "In three or more dimensions, the w velocity must be supplied!");
69 :
70 : // Strain tensor term requires velocity gradients;
71 71 : if (dynamic_cast<const MooseLinearVariableFV<Real> *>(&_u_var))
72 0 : requestVariableCellGradient(getParam<MooseFunctorName>("u"));
73 71 : if (dynamic_cast<const MooseLinearVariableFV<Real> *>(_v_var))
74 0 : requestVariableCellGradient(getParam<MooseFunctorName>("v"));
75 71 : if (dynamic_cast<const MooseLinearVariableFV<Real> *>(_w_var))
76 0 : requestVariableCellGradient(getParam<MooseFunctorName>("w"));
77 71 : }
78 :
79 : void
80 373 : LinearFVTKEDSourceSink::initialSetup()
81 : {
82 373 : LinearFVElementalKernel::initialSetup();
83 746 : NS::getWallBoundedElements(
84 373 : _wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _wall_bounded);
85 373 : NS::getWallDistance(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _dist);
86 373 : NS::getElementFaceArgs(_wall_boundary_names, _fe_problem, _subproblem, blockIDs(), _face_infos);
87 373 : }
88 :
89 : Real
90 633792 : LinearFVTKEDSourceSink::computeMatrixContribution()
91 : {
92 633792 : if (_wall_bounded.find(_current_elem_info->elem()) != _wall_bounded.end())
93 : // TKED value for near wall element will be directly assigned for this cell
94 : return 1.0;
95 : else
96 : {
97 : // Convenient definitions
98 469008 : const auto state = determineState();
99 469008 : const auto elem_arg = makeElemArg(_current_elem_info->elem());
100 469008 : const Real rho = _rho(elem_arg, state);
101 469008 : const Real TKE = _k(elem_arg, state);
102 469008 : const auto epsilon = _var.getElemValue(*_current_elem_info, state);
103 :
104 : // Compute destruction
105 469008 : const auto destruction = _C2_eps(elem_arg, state) * rho * epsilon / TKE;
106 :
107 : // Assign to matrix (term gets multiplied by TKED)
108 469008 : return destruction * _current_elem_volume;
109 : }
110 : }
111 :
112 : Real
113 633792 : LinearFVTKEDSourceSink::computeRightHandSideContribution()
114 : {
115 633792 : if (_wall_bounded.find(_current_elem_info->elem()) != _wall_bounded.end())
116 : {
117 : // Convenient definitions
118 164784 : const auto state = determineState();
119 164784 : const auto elem_arg = makeElemArg(_current_elem_info->elem());
120 164784 : const Real rho = _rho(elem_arg, state);
121 164784 : const Real mu = _mu(elem_arg, state);
122 164784 : const Real TKE = _k(elem_arg, state);
123 :
124 : // Convenient variables
125 : Real destruction = 0.0;
126 : std::vector<Real> y_plus_vec, velocity_grad_norm_vec;
127 : Real tot_weight = 0.0;
128 :
129 : // Get velocity vector
130 164784 : RealVectorValue velocity(_u_var(elem_arg, state));
131 164784 : if (_v_var)
132 164784 : velocity(1) = (*_v_var)(elem_arg, state);
133 164784 : if (_w_var)
134 0 : velocity(2) = (*_w_var)(elem_arg, state);
135 :
136 : // Get near wall faceInfo and distances from cell center to every wall
137 164784 : const auto & face_info_vec = libmesh_map_find(_face_infos, _current_elem_info->elem());
138 164784 : const auto & distance_vec = libmesh_map_find(_dist, _current_elem_info->elem());
139 : mooseAssert(distance_vec.size(), "Should have found a distance vector");
140 : mooseAssert(distance_vec.size() == face_info_vec.size(),
141 : "Should be as many distance vectors as face info vectors");
142 :
143 : // Update y+ and wall face cell
144 330208 : for (unsigned int i = 0; i < distance_vec.size(); i++)
145 : {
146 165424 : const auto distance = distance_vec[i];
147 : mooseAssert(distance > 0, "Should be at a non-zero distance");
148 :
149 : Real y_plus;
150 165424 : if (_wall_treatment == NS::WallTreatmentEnum::NEQ) // Non-equilibrium / Non-iterative
151 30000 : y_plus = distance * std::sqrt(std::sqrt(_C_mu) * TKE) * rho / mu;
152 : else // Equilibrium / Iterative
153 : {
154 135424 : const auto parallel_speed = NS::computeSpeed<Real>(
155 135424 : velocity - velocity * face_info_vec[i]->normal() * face_info_vec[i]->normal());
156 135424 : y_plus = NS::findyPlus<Real>(mu, rho, std::max(parallel_speed, 1e-10), distance);
157 : }
158 :
159 165424 : y_plus_vec.push_back(y_plus);
160 165424 : tot_weight += 1.0;
161 : }
162 :
163 : // Compute near wall epsilon value
164 330208 : for (const auto i : index_range(y_plus_vec))
165 : {
166 165424 : const auto y_plus = y_plus_vec[i];
167 :
168 165424 : if (y_plus < 11.25)
169 19429 : destruction += 2.0 * TKE * mu / rho / Utility::pow<2>(distance_vec[i]) / tot_weight;
170 : else
171 145995 : destruction += std::pow(_C_mu, 0.75) * std::pow(TKE, 1.5) /
172 145995 : (NS::von_karman_constant * distance_vec[i]) / tot_weight;
173 : }
174 :
175 : // Assign the computed value of TKED for element near the wall
176 : return destruction;
177 164784 : }
178 : else
179 : {
180 : // Convenient definitions
181 469008 : const auto state = determineState();
182 469008 : const auto elem_arg = makeElemArg(_current_elem_info->elem());
183 469008 : const Real rho = _rho(elem_arg, state);
184 469008 : const Real TKE = _k(elem_arg, state);
185 469008 : const Real TKED = _var.getElemValue(*_current_elem_info, state);
186 :
187 : // Compute production of TKE
188 469008 : const auto subdomain_id = _current_elem_info->elem()->subdomain_id();
189 469008 : const auto coord_sys = _subproblem.getCoordSystem(subdomain_id);
190 : const unsigned int rz_radial_coord =
191 469008 : coord_sys == Moose::COORD_RZ ? _subproblem.getAxisymmetricRadialCoord() : 0;
192 469008 : const auto symmetric_strain_tensor_sq_norm = NS::computeShearStrainRateNormSquared<Real>(
193 : _u_var, _v_var, _w_var, elem_arg, state, coord_sys, rz_radial_coord);
194 469008 : Real production = _mu_t(elem_arg, state) * symmetric_strain_tensor_sq_norm;
195 :
196 : // Limit TKE production (needed for flows with stagnation zones)
197 469008 : const Real production_limit = _C_pl * rho * TKED;
198 469008 : production = std::min(production, production_limit);
199 :
200 : // Compute production - recasted with mu_t definition to avoid division by epsilon
201 469008 : const auto production_epsilon = _C1_eps(elem_arg, state) * production * TKED / TKE;
202 :
203 : // Assign to matrix (term gets multiplied by TKED)
204 469008 : return production_epsilon * _current_elem_volume;
205 : }
206 : }
|
