fix: set true calculation of reynolds number

calculate reynolds number exact value in indvel.m
use its exact value when possible, if not use an approximation in ElemPerf.m

src/classes/@ElemPerf/ElemPerf.m

@@ -56,7 +56,6 @@ classdef ElemPerf < handle
         Rot (1, 1) Rotor % Handle of the rotor linked to this specific ElemPerf instance
         Op (1, 1)  Oper  % Handle of the Operating point linked to this specific ElemPerf instance
 
-        reynolds (1, :) double % Reynolds number, [-]
         upstreamVelAx (1, :) double % Upstream axial velocity above rotor disk, [m/s]
         upstreamVelTg (1, :) double % Upstream tangential velocity above rotor disk, [m/s]
 
@@ -85,6 +84,7 @@ classdef ElemPerf < handle
     % from a single set method otherwise.
     properties (Dependent)
         alpha    (1, :) double % Angle of attack, [rad]
+        reynolds (1, :) double % Reynolds number, [-]
         indVelAx (1, :) double % Induced axial velocity at rotor disk, [m/s]
         indVelTg (1, :) double % Induced tangential velocity at rotor disk, [m/s]
         inflowRat (1, :) double % Inflow ratio, [-]
@@ -96,6 +96,7 @@ classdef ElemPerf < handle
     % Cache for dependent properties to avoid excessive recalculation
     properties (Access = private, Hidden)
         alpha_    (1, :) double % Angle of attack, [rad]
+        reynolds_    (1, :) double % Reynolds number, [-]
         indVelAx_ (1, :) double % Induced axial velocity at rotor disk, [m/s]
         indVelTg_ (1, :) double % Induced tangential velocity at rotor disk, [m/s]
         inflowRat_ (1, :) double % Inflow ratio, [-]
@@ -128,10 +129,6 @@ classdef ElemPerf < handle
                 % In-plane velocities
                 self.tgSpeed = Op.omega .* Rot.Bl.y;
 
-                % Reynolds
-                relVel = sqrt(self.tgSpeed.^2 + Op.speed.^2);
-                self.reynolds = Flow.reynolds(relVel, Rot.Bl.chord, Op.Flow.mu, Op.Flow.rho);
-
                 % Get alpha_0 from reynolds
                 self.alpha0 = zeros(size(self.pitch));
                 for i = 1:length(Rot.Af)
@@ -162,6 +159,18 @@ classdef ElemPerf < handle
             alpha = self.alpha_;
         end
 
+        function reynolds = get.reynolds(self)
+            % if true reynolds is unknown, we use an approximated value
+            % with the upstream velocity instead of relative velocity at
+            % the rotor disk
+            if isempty(self.reynolds_)
+                relVel = sqrt(self.tgSpeed.^2 + self.upstreamVelAx.^2);
+                reynolds = Flow.reynolds(relVel, Rot.Bl.chord, Op.Flow.mu, Op.Flow.rho);
+            else
+                reynolds = self.reynolds_;
+            end
+        end
+
         function indVelAx = get.indVelAx(self)
             indVelAx = self.indVelAx_;
         end
@@ -197,6 +206,11 @@ classdef ElemPerf < handle
             [self.cl, self.cd] = getclcd(self, val);
         end
 
+        function set.reynolds(self, val)
+            % Cache value
+            self.reynolds_ = val;
+        end
+
         function set.indVelAx(self, val)
             % Cache value
             self.indVelAx_ = val;

src/solvers/indvel.m

@@ -68,15 +68,17 @@ function indvel(OpRot, Mod)
         % Local mass flow rate
         dmdot = 2 * pi * OpRot.Op.Flow.rho * OpRot.Rot.Bl.y * OpRot.Rot.Bl.dy .* v;
 
-        % Compute new estimates for the velocity magnitude and flow angle
+        % Compute new estimates for the velocity magnitude and Reynolds
+        % number
         relVel = sqrt(v.^2 + u.^2);
+        reynolds = Flow.reynolds(relVel, OpRot.Rot.Bl.chord, OpRot.Op.Flow.mu, OpRot.Op.Flow.rho);
 
         % Angles
         phi = atan2(v, u);  % Inflow angle
         alpha = OpRot.ElPerf.truePitch - phi;
 
         % Get new estimates for the coefficients
-        [cl, cd] = OpRot.ElPerf.getclcd(alpha);
+        [cl, cd] = OpRot.ElPerf.getclcd(alpha,reynolds);
 
         % Loss factor (separate swirl and axial components)
         lossFact = prandtlloss(OpRot.Rot.nBlades, OpRot.Rot.Bl.r, OpRot.Rot.r0, phi, ...
@@ -118,6 +120,7 @@ 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