fix(solvers): true Reynolds calculation in all solvers (closes #13)

This completes the work started in !4 by implementing the true
Reynolds recalculation for the other three solvers.

Note that special care had to be implemented in Stahlhut solver in order
to ensure proper numerical resolution of the solution.

src/configs/caradonna1981.m

@@ -49,7 +49,6 @@ Sim.Misc.appli  = 'heli'; % Type of application ('helicopter', 'propeller', 'win
 
 % Solvers
 Mod.solvers = {'all'};   % BEMT Solver ('leishman', 'indfact', 'indvel', 'stahlhut', 'all')
-% Mod.solvers = {'leishman'};   % BEMT Solver ('leishman', 'indfact', 'indvel', 'stahlhut', 'all')
 
 % Extensions/corrections
 Mod.Ext.losses = 'all';  % Include losses using Prandtl formula ('none', 'hub', 'tip', 'both')

src/solvers/indfact.m

@@ -77,13 +77,14 @@ function indfact(OpRot, Mod)
         % [DEBUG]: SQRT(PI) GIVES BETTER MATCH?
         % relVel = sqrt(sqrt(pi)) * sqrt(v_ax.^2 + v_ang.^2);
         relVel = sqrt(v_ax.^2 + v_ang.^2);
+        reynolds = Flow.reynolds(relVel, OpRot.Rot.Bl.chord, OpRot.Op.Flow.mu, OpRot.Op.Flow.rho);
 
         % Angles
         phi = atan2(v_ax, v_ang);  % Inflow angle
         alpha = OpRot.ElPerf.truePitch - phi;  % Angle of attack
 
         % Get new estimates for the coefficients
-        [cl, cd] = OpRot.ElPerf.getclcd(alpha);
+        [cl, cd] = OpRot.ElPerf.getclcd(alpha, reynolds);
 
         % Calculate loss factor
         lossFact = prandtlloss(OpRot.Rot.nBlades, OpRot.Rot.Bl.r, OpRot.Rot.r0, phi, ...
@@ -133,5 +134,6 @@ function indfact(OpRot, Mod)
     OpRot.ElPerf.indVelTg = b .* OpRot.ElPerf.tgSpeed;
     OpRot.ElPerf.inflAngle = phi;
     OpRot.ElPerf.alpha = alpha;
+    OpRot.ElPerf.reynolds = reynolds;
 
 end

src/solvers/indvel.m

@@ -68,8 +68,7 @@ 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 Reynolds
-        % number
+        % 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);
 

src/solvers/leishman.m

@@ -68,6 +68,13 @@ function leishman(OpRot, Mod)
         while ~converged
             lambda_old = lambda;
 
+            % Update Reynolds and then clSlope w.r.t. the new inflow velocity
+            axialVel = OpRot.Op.speed + lambda * OpRot.tgTipSpeed;
+            relVel = sqrt(axialVel.^2 + OpRot.ElPerf.tgSpeed.^2);
+            reynolds = Flow.reynolds(relVel, ...
+                                     OpRot.Rot.Bl.chord, OpRot.Op.Flow.mu, OpRot.Op.Flow.rho);
+            clSlope = getclslope(OpRot.Rot.Af, OpRot.Rot.Bl.iAf, reynolds);
+
             % Calculate loss factor
             lossFact = prandtlloss(OpRot.Rot.nBlades, OpRot.Rot.Bl.r, OpRot.Rot.r0, ...
                                    phi, Mod.Ext.losses);
@@ -93,6 +100,7 @@ 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
 

src/solvers/stahlhut.m

@@ -40,6 +40,8 @@ function stahlhut(OpRot, Mod)
     % Ref: Stahlhut and Leishman, "Aerodynamic design optimization of proprotors for convertible-
     %      rotor concepts", In American Helicopter Society 68th Annual Forum. Fort Worth, TX, USA.
     % ----------------------------------------------------------------------------------------------
+    % FIXME: (L) Find a way to solve this for the whole vector at a time to improve CPU time.
+    % ----------------------------------------------------------------------------------------------
     % (c) Copyright 2022-2023 University of Liege
     % Author: Thomas Lambert <t.lambert@uliege.be>
     % ULiege - Aeroelasticity and Experimental Aerodynamics
@@ -122,6 +124,12 @@ 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)
@@ -129,16 +137,40 @@ function [gphi, b1phi, b2phi] = stahlhuteq(i, phi, OpRot, lossType)
 
     % Defaults and constants
     SIGNUM_ZERO = 1; % Value of the signum function in 0;
+    MAX_REYNOLDS_DIFF = 50; % Maximum difference between the true the approximate Reynolds, [%]
+
+    persistent elemIdx indVelAx indVelTg
 
     % Abbreviations
     r = OpRot.Rot.Bl.r(i); % Relative position
     sol = OpRot.Rot.solidity(i); % Solidity
 
-    reynolds = OpRot.ElPerf.reynolds(i);
+    reynolds = OpRot.ElPerf.reynolds(i); % Use approximate Reynolds by default
     tgSpeed = OpRot.ElPerf.tgSpeed(i);
     pitch = OpRot.ElPerf.truePitch(i);
     airspeed = OpRot.Op.speed;
 
+    % Recalculate Reynolds based on actual inflow velocity
+    %   The first time we calculate it for a given element, we use the approximate value that
+    %   neglects the induced velocity. Afterwards, the reynolds is calculated based on the
+    %   previously found inflow velocity
+    if isempty(elemIdx) || elemIdx ~= i
+        elemIdx = i;
+    else
+        relVel = sqrt((airspeed + indVelAx)^2 + (tgSpeed - indVelTg)^2);
+        tmpReynolds = Flow.reynolds(relVel, ...
+                                    OpRot.Rot.Bl.chord(i), OpRot.Op.Flow.mu, OpRot.Op.Flow.rho);
+
+        % Stabilize first steps of numeric resolution
+        %   The main contirbutions to the Reynolds are the external velocity and the rotation speed
+        %   If the calculated value deviates too much to the approximate one, reject it and keep the
+        %   appximation. Once the solution is a bit more stable, the true Reynolds can start to be
+        %   used.
+        if abs(tmpReynolds - reynolds) / reynolds * 100 <= MAX_REYNOLDS_DIFF
+            reynolds = tmpReynolds;
+        end
+    end
+
     % Loss factor (separate swirl and axial components)
     lossFact = prandtlloss(OpRot.Rot.nBlades, r, OpRot.Rot.r0, phi, lossType);
     K_T = 1 - (1 - lossFact) .* cos(phi);
@@ -159,6 +191,11 @@ function [gphi, b1phi, b2phi] = stahlhuteq(i, phi, OpRot, lossType)
            sgn(phi, SIGNUM_ZERO) * 1 / (8 * r) * sol * cl * sec(gamma) * ...
            (tgSpeed ./ K_T .* cos(phi + gamma) + airspeed ./ K_P .* sin(phi + gamma));
 
+    % Current value for the induced velocities (needed for true Reynolds)
+    xi = r .* cos(phi) ./ b2phi;
+    lambda = xi  .* tan(phi);
+    indVelAx = lambda .* OpRot.tgTipSpeed;
+    indVelTg = tgSpeed - xi * OpRot.tgTipSpeed;
 end
 
 function plotgfunc(i, allphi, allg, unk0)