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 "INSFVTurbulentTemperatureWallFunction.h"
11 : #include "NavierStokesMethods.h"
12 : #include "NS.h"
13 :
14 : registerADMooseObject("NavierStokesApp", INSFVTurbulentTemperatureWallFunction);
15 :
16 : InputParameters
17 366 : INSFVTurbulentTemperatureWallFunction::validParams()
18 : {
19 366 : InputParameters params = FVFluxBC::validParams();
20 :
21 366 : params.addClassDescription("Adds turbulent temperature wall function.");
22 732 : params.addRequiredParam<MooseFunctorName>("T_w", "The wall temperature.");
23 732 : params.addRequiredParam<MooseFunctorName>("u", "The velocity in the x direction.");
24 732 : params.addParam<MooseFunctorName>("v", "The velocity in the y direction.");
25 732 : params.addParam<MooseFunctorName>("w", "The velocity in the z direction.");
26 366 : params.addRequiredParam<MooseFunctorName>(NS::density, "Density");
27 366 : params.addRequiredParam<MooseFunctorName>(NS::mu, "Dynamic viscosity.");
28 366 : params.addRequiredParam<MooseFunctorName>(NS::cp, "The specific heat at constant pressure.");
29 366 : params.addParam<MooseFunctorName>(NS::kappa, "The thermal conductivity.");
30 732 : params.addParam<MooseFunctorName>("Pr_t", 0.58, "The turbulent Prandtl number.");
31 366 : params.addParam<MooseFunctorName>(NS::TKE, "Turbulent kinetic energy functor.");
32 732 : params.addParam<Real>("C_mu", 0.09, "Coupled turbulent kinetic energy closure.");
33 732 : MooseEnum wall_treatment("eq_newton eq_incremental eq_linearized neq", "neq");
34 732 : params.addParam<MooseEnum>(
35 : "wall_treatment", wall_treatment, "The method used for computing the wall functions");
36 :
37 732 : params.addParam<bool>("newton_solve", false, "Whether a Newton nonlinear solve is being used");
38 732 : params.addParamNamesToGroup("newton_solve", "Advanced");
39 366 : return params;
40 366 : }
41 :
42 198 : INSFVTurbulentTemperatureWallFunction::INSFVTurbulentTemperatureWallFunction(
43 198 : const InputParameters & parameters)
44 : : FVFluxBC(parameters),
45 198 : _dim(_subproblem.mesh().dimension()),
46 396 : _T_w(getFunctor<ADReal>("T_w")),
47 396 : _u_var(getFunctor<ADReal>("u")),
48 594 : _v_var(parameters.isParamValid("v") ? &(getFunctor<ADReal>("v")) : nullptr),
49 198 : _w_var(parameters.isParamValid("w") ? &(getFunctor<ADReal>("w")) : nullptr),
50 198 : _rho(getFunctor<ADReal>(NS::density)),
51 198 : _mu(getFunctor<ADReal>(NS::mu)),
52 198 : _cp(getFunctor<ADReal>(NS::cp)),
53 198 : _kappa(getFunctor<ADReal>(NS::kappa)),
54 396 : _Pr_t(getFunctor<ADReal>("Pr_t")),
55 198 : _k(getFunctor<ADReal>(NS::TKE)),
56 396 : _C_mu(getParam<Real>("C_mu")),
57 396 : _wall_treatment(getParam<MooseEnum>("wall_treatment")),
58 594 : _newton_solve(getParam<bool>("newton_solve"))
59 : {
60 198 : }
61 :
62 : ADReal
63 16984 : INSFVTurbulentTemperatureWallFunction::computeQpResidual()
64 : {
65 : using std::abs, std::log, std::max, std::sqrt, std::pow, std::exp;
66 :
67 : // Useful parameters
68 16984 : const FaceInfo & fi = *_face_info;
69 : // if on an internal face (internal to the mesh, but an external boundary of the flow area),
70 : // we have to select the element on which the flow variables / temperature are defined
71 16984 : const bool use_elem = (_face_type == FaceInfo::VarFaceNeighbors::ELEM);
72 : const Real wall_dist = abs(
73 16984 : ((use_elem ? fi.elemCentroid() : fi.neighborCentroid()) - fi.faceCentroid()) * fi.normal());
74 16984 : const Elem * const elem_ptr = use_elem ? fi.elemPtr() : fi.neighborPtr();
75 16984 : const auto elem_arg = makeElemArg(elem_ptr);
76 16984 : const auto state = determineState();
77 : const auto old_state = Moose::StateArg(1, Moose::SolutionIterationType::Nonlinear);
78 16984 : const auto mu = _mu(elem_arg, state);
79 16984 : const auto rho = _rho(elem_arg, state);
80 16984 : const auto cp = _cp(elem_arg, state);
81 16984 : const auto kappa = _kappa(elem_arg, state);
82 :
83 : // Get the velocity vector
84 33968 : ADRealVectorValue velocity(_u_var(elem_arg, state));
85 16984 : if (_v_var)
86 16984 : velocity(1) = (*_v_var)(elem_arg, state);
87 16984 : if (_w_var)
88 0 : velocity(2) = (*_w_var)(elem_arg, state);
89 :
90 : // Compute the velocity and direction of the velocity component that is parallel to the wall
91 : const ADReal parallel_speed =
92 50952 : NS::computeSpeed<ADReal>(velocity - velocity * (fi.normal()) * (fi.normal()));
93 :
94 : // Computing friction velocity and y+
95 : ADReal u_tau, y_plus;
96 :
97 16984 : if (_wall_treatment == "eq_newton")
98 : {
99 : // Full Newton-Raphson solve to find the wall quantities from the law of the wall
100 0 : u_tau = NS::findUStar<ADReal>(mu, rho, parallel_speed, wall_dist);
101 0 : y_plus = wall_dist * u_tau * rho / mu;
102 : }
103 16984 : else if (_wall_treatment == "eq_incremental")
104 : {
105 : // Incremental solve on y_plus to get the near-wall quantities
106 0 : y_plus = NS::findyPlus<ADReal>(mu, rho, max(parallel_speed, 1e-10), wall_dist);
107 0 : u_tau = parallel_speed /
108 0 : (log(max(NS::E_turb_constant * y_plus, 1.0 + 1e-4)) / NS::von_karman_constant);
109 : }
110 16984 : else if (_wall_treatment == "eq_linearized")
111 : {
112 : // Linearized approximation to the wall function to find the near-wall quantities faster
113 16984 : const ADReal a_c = 1 / NS::von_karman_constant;
114 : const ADReal b_c =
115 84920 : 1.0 / NS::von_karman_constant * (log(NS::E_turb_constant * max(wall_dist, 1.0) / mu) + 1.0);
116 16984 : const ADReal c_c = parallel_speed;
117 :
118 84920 : u_tau = (-b_c + sqrt(pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c);
119 16984 : y_plus = wall_dist * u_tau * rho / mu;
120 : }
121 0 : else if (_wall_treatment == "neq")
122 : {
123 : // Assign non-equilibrium wall function value
124 0 : y_plus = wall_dist * sqrt(sqrt(_C_mu) * _k(elem_arg, old_state)) * rho / mu;
125 0 : u_tau = parallel_speed /
126 0 : (log(max(NS::E_turb_constant * y_plus, 1.0 + 1e-4)) / NS::von_karman_constant);
127 : }
128 :
129 : ADReal alpha;
130 16984 : if (y_plus <= 5.0) // sub-laminar layer
131 : {
132 16723 : alpha = kappa / (rho * cp);
133 : }
134 261 : else if (y_plus >= 30.0) // log-layer
135 : {
136 0 : const auto Pr = cp * mu / kappa;
137 0 : const auto Pr_ratio = Pr / _Pr_t(elem_arg, state);
138 : const auto jayatilleke_P =
139 0 : 9.24 * (pow(Pr_ratio, 0.75) - 1.0) * (1.0 + 0.28 * exp(-0.007 * Pr_ratio));
140 : const auto wall_scaling =
141 0 : 1.0 / NS::von_karman_constant * log(NS::E_turb_constant * y_plus) + jayatilleke_P;
142 0 : alpha = u_tau * wall_dist / (_Pr_t(elem_arg, state) * wall_scaling);
143 : }
144 : else // buffer layer
145 : {
146 261 : const auto alpha_lam = kappa / (rho * cp);
147 261 : const auto Pr = cp * mu / kappa;
148 261 : const auto Pr_t = _Pr_t(elem_arg, state);
149 : const auto Pr_ratio = Pr / Pr_t;
150 : const auto jayatilleke_P =
151 1566 : 9.24 * (pow(Pr_ratio, 0.75) - 1.0) * (1.0 + 0.28 * exp(-0.007 * Pr_ratio));
152 : const auto wall_scaling =
153 522 : 1.0 / NS::von_karman_constant * log(NS::E_turb_constant * y_plus) + jayatilleke_P;
154 261 : const auto alpha_turb = u_tau * wall_dist / (Pr_t * wall_scaling);
155 261 : const auto blending_function = (y_plus - 5.0) / 25.0;
156 522 : alpha = blending_function * alpha_turb + (1.0 - blending_function) * alpha_lam;
157 : }
158 :
159 : // To make sure new derivatives are not introduced as the solve progresses
160 16984 : if (_newton_solve)
161 33672 : alpha += 0 * kappa * (rho * cp) + 0 * u_tau * _Pr_t(elem_arg, state) * y_plus;
162 :
163 16984 : const auto face_arg = singleSidedFaceArg();
164 50952 : return -rho * cp * alpha * (_T_w(face_arg, state) - _var(elem_arg, state)) / wall_dist;
165 : }
|
