diff --git a/src/configs/caradonna1981.m b/src/configs/caradonna1981.m
index 017a0cb677021b9d66375d35fee41c2f0bd368a2..e96cacec84a4321630bff298023a2b8a82540f25 100644
--- a/src/configs/caradonna1981.m
+++ b/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')
diff --git a/src/solvers/indfact.m b/src/solvers/indfact.m
index c5e19dc200533f8206df4a7f5fa806376f269788..54ed6d400804e9e8ea036f367c6b6cd2bdaa250c 100644
--- a/src/solvers/indfact.m
+++ b/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
diff --git a/src/solvers/indvel.m b/src/solvers/indvel.m
index 2ac1ac1ba4128880847ab1f2d6eaed4c08f27172..2eff70d92bfdb1491b8bcedc6e8befce7ae3d932 100644
--- a/src/solvers/indvel.m
+++ b/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);
diff --git a/src/solvers/leishman.m b/src/solvers/leishman.m
index b7835280b3ae2ddfbdbdbd66074466396a6eff16..0160198410f7b1c5b69b92344c107fdb788ec294 100644
--- a/src/solvers/leishman.m
+++ b/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
diff --git a/src/solvers/stahlhut.m b/src/solvers/stahlhut.m
index 5a8a6898c751e3fd5c2bf6499926a2566c11c608..e79e1cb2a335537016ddda1d5d08b43447821ef2 100644
--- a/src/solvers/stahlhut.m
+++ b/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)