From 93480ba3fce295c822814d812c8d293899e65e6d Mon Sep 17 00:00:00 2001 From: Freile Date: Wed, 12 Jun 2024 14:45:10 -0600 Subject: [PATCH] Fix findYPlus routine. #27800 --- .../navier_stokes/src/base/NavierStokesMethods.C | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/modules/navier_stokes/src/base/NavierStokesMethods.C b/modules/navier_stokes/src/base/NavierStokesMethods.C index 8a61a259d7cb..42f1cef24170 100644 --- a/modules/navier_stokes/src/base/NavierStokesMethods.C +++ b/modules/navier_stokes/src/base/NavierStokesMethods.C @@ -69,13 +69,20 @@ findUStar(const ADReal & mu, const ADReal & rho, const ADReal & u, const Real di const ADReal nu = mu / rho; - ADReal u_star = std::sqrt(nu * u / dist); + // Wall-function linearized guess + const Real a_c = 1 / NS::von_karman_constant; + const ADReal b_c = + 1 / NS::von_karman_constant * (std::log(NS::E_turb_constant * dist / mu) + 1.0); + const ADReal c_c = u; + + /// This is important to reduce the number of nonlinear iterations + ADReal u_star = (-b_c + std::sqrt(std::pow(b_c, 2) + 4.0 * a_c * c_c)) / (2.0 * a_c); // Newton-Raphson method to solve for u_star (friction velocity). for (int i = 0; i < MAX_ITERS; ++i) { - ADReal residual = u_star / NS::von_karman_constant * std::log(u_star * dist / (0.111 * nu)) - u; - ADReal deriv = (1 + std::log(u_star * dist / (0.111 * nu))) / NS::von_karman_constant; + ADReal residual = u_star / NS::von_karman_constant * std::log(NS::E_turb_constant * u_star * dist / nu) - u; + ADReal deriv = (1 + std::log(NS::E_turb_constant * u_star * dist / nu)) / NS::von_karman_constant; ADReal new_u_star = std::max(1e-20, u_star - residual / deriv); Real rel_err = std::abs((new_u_star.value() - u_star.value()) / new_u_star.value()); @@ -114,7 +121,7 @@ findyPlus(const ADReal & mu, const ADReal & rho, const ADReal & u, const Real di { yPlusLast = yPlus; yPlus = (kappa_time_Re + yPlus) / (1.0 + std::log(NS::E_turb_constant * yPlus)); - } while (rev_yPlusLam * (yPlus - yPlusLast) > REL_TOLERANCE && ++iters < MAX_ITERS); + } while ( std::abs(rev_yPlusLam * (yPlus - yPlusLast)) > REL_TOLERANCE && ++iters < MAX_ITERS); return std::max(0.0, yPlus); }