From 4781193e314c831708cdf42356e354bf26faec44 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Mon, 8 Jan 2018 10:33:56 +0000
Subject: [PATCH] ENH: nutUBlended wall function - added protection for uTau. 
 See #696

---
 ...utUBlendedWallFunctionFvPatchScalarField.C | 23 +++++++++++--------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C
index 81f1d10a259..47f25c6f0de 100644
--- a/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C
+++ b/src/TurbulenceModels/turbulenceModels/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUBlendedWallFunction/nutUBlendedWallFunctionFvPatchScalarField.C
@@ -93,17 +93,20 @@ Foam::nutUBlendedWallFunctionFvPatchScalarField::calcUTau
     forAll(uTaup, facei)
     {
         scalar ut = sqrt((nutw[facei] + nuw[facei])*magGradU[facei]);
-        scalar error = GREAT;
-        label iter = 0;
-        while (iter++ < 10 && error > 0.001)
+        if (mag(ut) > ROOTVSMALL)
         {
-            scalar yPlus = y[facei]*ut/nuw[facei];
-            scalar uTauVis = magUp[facei]/yPlus;
-            scalar uTauLog = kappa_*magUp[facei]/log(E_*yPlus);
-
-            scalar utNew = pow(pow(uTauVis, n_) + pow(uTauLog, n_), 1.0/n_);
-            error = mag(ut - utNew)/(ut + ROOTVSMALL);
-            ut = 0.5*(ut + utNew);
+            scalar error = GREAT;
+            label iter = 0;
+            while (iter++ < 10 && error > 0.001)
+            {
+                scalar yPlus = y[facei]*ut/nuw[facei];
+                scalar uTauVis = magUp[facei]/yPlus;
+                scalar uTauLog = kappa_*magUp[facei]/log(E_*yPlus);
+
+                scalar utNew = pow(pow(uTauVis, n_) + pow(uTauLog, n_), 1.0/n_);
+                error = mag(ut - utNew)/(ut + ROOTVSMALL);
+                ut = 0.5*(ut + utNew);
+            }
         }
         uTaup[facei] = ut;
     }
-- 
GitLab