From 0f0a012a6b6992b22c62a37372598f5dcf3bd045 Mon Sep 17 00:00:00 2001
From: Thomas Lambert <t.lambert@uliege.be>
Date: Fri, 28 Jul 2023 16:58:02 +0200
Subject: [PATCH] chore(reynolds): more accurate Reynolds recalculation

---
 src/classes/@ElemPerf/ElemPerf.m | 25 +++++++++++++------------
 src/solvers/indfact.m            |  2 --
 src/solvers/indvel.m             |  2 --
 src/solvers/leishman.m           |  2 --
 src/solvers/stahlhut.m           |  6 ------
 5 files changed, 13 insertions(+), 24 deletions(-)

diff --git a/src/classes/@ElemPerf/ElemPerf.m b/src/classes/@ElemPerf/ElemPerf.m
index d17cb07..e30aafb 100644
--- a/src/classes/@ElemPerf/ElemPerf.m
+++ b/src/classes/@ElemPerf/ElemPerf.m
@@ -160,15 +160,20 @@ classdef ElemPerf < handle
         end
 
         function reynolds = get.reynolds(self)
-            % If Reynolds in unknown, we use an approximate solution
-            % neglecting induced velocities.
-            if isempty(self.reynolds_)
-                relVel = sqrt(self.tgSpeed.^2 + self.upstreamVelAx.^2);
-                reynolds = Flow.reynolds(relVel, self.Rot.Bl.chord, ...
-                                         self.Op.Flow.mu, self.Op.Flow.rho);
-            else
-                reynolds = self.reynolds_;
+            % Always recalculate Reynolds based on the most accurate data
+            velAx = self.upstreamVelAx;
+            velTg = self.tgSpeed;
+
+            if ~isempty(self.indVelAx)
+                velAx = velAx + self.indVelAx;
+            end
+            if ~isempty(self.indVelTg)
+                velTg =  velTg - self.indVelTg;
             end
+            relVel = sqrt(velAx.^2 + velTg.^2);
+            reynolds = Flow.reynolds(relVel, self.Rot.Bl.chord, ...
+                                     self.Op.Flow.mu, self.Op.Flow.rho);
+
         end
 
         function indVelAx = get.indVelAx(self)
@@ -206,10 +211,6 @@ classdef ElemPerf < handle
             [self.cl, self.cd] = getclcd(self, val);
         end
 
-        function set.reynolds(self, val)
-            self.reynolds_ = val;
-        end
-
         function set.indVelAx(self, val)
             % Cache value
             self.indVelAx_ = val;
diff --git a/src/solvers/indfact.m b/src/solvers/indfact.m
index 54ed6d4..11ac1ea 100644
--- a/src/solvers/indfact.m
+++ b/src/solvers/indfact.m
@@ -134,6 +134,4 @@ function indfact(OpRot, Mod)
     OpRot.ElPerf.indVelTg = b .* OpRot.ElPerf.tgSpeed;
     OpRot.ElPerf.inflAngle = phi;
     OpRot.ElPerf.alpha = alpha;
-    OpRot.ElPerf.reynolds = reynolds;
-
 end
diff --git a/src/solvers/indvel.m b/src/solvers/indvel.m
index 2eff70d..580da9a 100644
--- a/src/solvers/indvel.m
+++ b/src/solvers/indvel.m
@@ -119,8 +119,6 @@ function indvel(OpRot, Mod)
     OpRot.ElPerf.alpha = alpha;
     OpRot.ElPerf.indVelAx = v - OpRot.ElPerf.upstreamVelAx;
     OpRot.ElPerf.indVelTg = uw / 2;
-    OpRot.ElPerf.reynolds = reynolds;
-
 end
 
 function force = coeff2force(coeff, OpRot, speed)
diff --git a/src/solvers/leishman.m b/src/solvers/leishman.m
index 0160198..54b6978 100644
--- a/src/solvers/leishman.m
+++ b/src/solvers/leishman.m
@@ -100,8 +100,6 @@ function leishman(OpRot, Mod)
     OpRot.ElPerf.inflAngle = phi;
     OpRot.ElPerf.indVelAx = lambda * OpRot.tgTipSpeed - OpRot.Op.speed;
     OpRot.ElPerf.indVelTg = OpRot.ElPerf.tgSpeed - xi * OpRot.tgTipSpeed;
-    OpRot.ElPerf.reynolds = reynolds;
-
 end
 
 % --------------------------------------------------------------------------------------------------
diff --git a/src/solvers/stahlhut.m b/src/solvers/stahlhut.m
index e79e1cb..15b8dce 100644
--- a/src/solvers/stahlhut.m
+++ b/src/solvers/stahlhut.m
@@ -124,12 +124,6 @@ function stahlhut(OpRot, Mod)
     OpRot.ElPerf.swirlRat = OpRot.Rot.Bl.r .* cos(OpRot.ElPerf.inflAngle) ./ b2phi;
     OpRot.ElPerf.inflowRat = OpRot.ElPerf.swirlRat .* tan(OpRot.ElPerf.inflAngle);
 
-    indVelAx = OpRot.ElPerf.inflowRat .* OpRot.tgTipSpeed;
-    indVelTg = OpRot.ElPerf.tgSpeed - OpRot.ElPerf.swirlRat .* OpRot.tgTipSpeed;
-    relVel = sqrt((OpRot.Op.speed + indVelAx).^2 + (OpRot.ElPerf.tgSpeed - indVelTg).^2);
-    OpRot.ElPerf.reynolds = Flow.reynolds(relVel, ...
-                                          OpRot.Rot.Bl.chord, OpRot.Op.Flow.mu, OpRot.Op.Flow.rho);
-
 end
 
 function [gphi, b1phi, b2phi] = stahlhuteq(i, phi, OpRot, lossType)
-- 
GitLab