From fd4b21279574b8f02ee286e9148987559a07476a Mon Sep 17 00:00:00 2001
From: Thomas Lambert <t.lambert@uliege.be>
Date: Fri, 28 Jul 2023 16:28:02 +0200
Subject: [PATCH] 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 | 1 -
src/solvers/indfact.m | 4 +++-
src/solvers/indvel.m | 3 +--
src/solvers/leishman.m | 8 ++++++++
src/solvers/stahlhut.m | 39 ++++++++++++++++++++++++++++++++++++-
5 files changed, 50 insertions(+), 5 deletions(-)
diff --git a/src/configs/caradonna1981.m b/src/configs/caradonna1981.m
index 017a0cb..e96cace 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 c5e19dc..54ed6d4 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 2ac1ac1..2eff70d 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 b783528..0160198 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 5a8a689..e79e1cb 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)
--
GitLab