From af5459537b25c15bab8a284676422cd41091a554 Mon Sep 17 00:00:00 2001
From: Rakotondratsimba Lyraie <liz_tp@proton.me>
Date: Tue, 19 Sep 2023 14:08:07 +0000
Subject: [PATCH] feat: add SMST coaxial model
- Adapt upstream velocity definition
- Proper calculation of inflow velocity for second rotor
- Modify Results class to store multiple rotors
- Adapt BEMT indvel
---
src/classes/@Blade/Blade.m | 2 +
src/classes/@ElemPerf/ElemPerf.m | 12 +++-
src/classes/@ElemPerf/calcforces.m | 2 +-
src/classes/@ElemPerf/updateupstreamvel.m | 85 +++++++++++++++++++++++
src/classes/@OperRotor/OperRotor.m | 6 ++
src/classes/@Result/Result.m | 23 +++---
src/classes/@Result/filter.m | 4 +-
src/configs/caradonna1981.m | 1 +
src/configs/knight1937.m | 1 +
src/configs/template.m | 1 +
src/configs/templatecoax.m | 1 +
src/rotare.m | 2 +-
src/solvers/bemt.m | 8 +--
src/solvers/indvel.m | 2 +
14 files changed, 131 insertions(+), 19 deletions(-)
create mode 100644 src/classes/@ElemPerf/updateupstreamvel.m
diff --git a/src/classes/@Blade/Blade.m b/src/classes/@Blade/Blade.m
index 19c5462..e8509d7 100644
--- a/src/classes/@Blade/Blade.m
+++ b/src/classes/@Blade/Blade.m
@@ -44,6 +44,7 @@ classdef Blade < handle
% twist : Twist of the guide stations, [deg] (same size as rad)
% iAf : Index of the airfoil for the guide stations, [-] (same size as rad)
% nElem : Number of elements to mesh the whole blade, [-]
+ % spinDir : Rotational direction of the current rotor, [-]
%
% See also: Rotor, rotare, template, af_tools.Airfoil.
%
@@ -69,6 +70,7 @@ classdef Blade < handle
chord (1, :) double {mustBePositive} % Element chord, [m]
twist (1, :) double {mustBeFinite} % Element twist (also called stagger angle), [rad]
iAf (1, :) double {mustBePositive} % Index of airfoil to use for each element, [-]
+ spinDir (1, 1) % Rotational direction of the current rotor, [-]
end
properties (SetAccess = private, Dependent)
diff --git a/src/classes/@ElemPerf/ElemPerf.m b/src/classes/@ElemPerf/ElemPerf.m
index e30aafb..68a75e7 100644
--- a/src/classes/@ElemPerf/ElemPerf.m
+++ b/src/classes/@ElemPerf/ElemPerf.m
@@ -85,6 +85,7 @@ classdef ElemPerf < handle
properties (Dependent)
alpha (1, :) double % Angle of attack, [rad]
reynolds (1, :) double % Reynolds number, [-]
+ massflow (1, :) double % Mass flow rate through the element, [kg/s]
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,7 +97,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, [-]
+ massflow_ (1, :) double % Mass flow rate through the element, [kg/s]
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, [-]
@@ -176,6 +177,10 @@ classdef ElemPerf < handle
end
+ function massflow = get.massflow(self)
+ massflow = self.massflow_;
+ end
+
function indVelAx = get.indVelAx(self)
indVelAx = self.indVelAx_;
end
@@ -211,6 +216,10 @@ classdef ElemPerf < handle
[self.cl, self.cd] = getclcd(self, val);
end
+ function set.massflow(self, val)
+ self.massflow_ = val;
+ end
+
function set.indVelAx(self, val)
% Cache value
self.indVelAx_ = val;
@@ -271,6 +280,7 @@ classdef ElemPerf < handle
calcforces(self) % Calculate forces, torque and power
[cl, cd] = getclcd(self, aoaVect, reyVect, i) % Get values of cl and cd
plotveltriangles(self, nTriangles, varargin)
+ updateupstreamvel(self, PrevOpRot)
end
end
diff --git a/src/classes/@ElemPerf/calcforces.m b/src/classes/@ElemPerf/calcforces.m
index f6700dc..6d41711 100644
--- a/src/classes/@ElemPerf/calcforces.m
+++ b/src/classes/@ElemPerf/calcforces.m
@@ -39,7 +39,7 @@ function calcforces(self)
if ~isempty(self.cl)
axVel = self.upstreamVelAx + self.indVelAx;
- tgVel = self.tgSpeed - self.indVelTg;
+ tgVel = self.tgSpeed - self.Rot.Bl.spinDir * (self.upstreamVelTg + self.indVelTg);
relVel = sqrt(axVel.^2 + tgVel.^2);
% Non-dimensionalization factor
diff --git a/src/classes/@ElemPerf/updateupstreamvel.m b/src/classes/@ElemPerf/updateupstreamvel.m
new file mode 100644
index 0000000..cea2781
--- /dev/null
+++ b/src/classes/@ElemPerf/updateupstreamvel.m
@@ -0,0 +1,85 @@
+function updateupstreamvel(self, PrevOpRot)
+ % UPDATEUPSTREAMVEL Update axial and tangential upstream velocities.
+ % This function is used to update upstream velocities of the current
+ % rotor based on the previous rotor inforations.
+ %
+ % Syntax:
+ % ElPerf.updateupstreamvel(PrevOpRot) updates upstream velocities in an ElPerf
+ % object according to the performance of the previous Operating Rotor `PrevOpRot`.
+ %
+ % Inputs:
+ % - prevOpRot : previous Operating Rotor object
+ %
+ % Outputs:
+ % Updates the following properties of the ElemPerf object:
+ % - self.upstreamVelAx : Upstream axial velocity above rotor disk, [m/s]
+ % - self.upstreamVelTg : Upstream tangential velocity above rotor disk, [m/s]
+ %
+ % It is assumed that :
+ % - a simplified multiple stream tube model is adopted
+ % - there is no contraction between upstream and rotor levels for the current rotor
+ % - the swirl is taken into account
+ %
+ % See also: ElemPerf.
+ %
+ % <a href="https:/gitlab.uliege.be/rotare/documentation">Complete documentation (online)</a>
+
+ % ----------------------------------------------------------------------------------------------
+ % (c) Copyright 2022-2023 University of Liege
+ % Author: Thomas Lambert <t.lambert@uliege.be>
+ % ULiege - Aeroelasticity and Experimental Aerodynamics
+ % MIT License
+ % Repo: https://gitlab.uliege.be/rotare/rotare
+ % Docs: https://gitlab.uliege.be/rotare/documentation
+ % Issues: https://gitlab.uliege.be/rotare/rotare/-/issues
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ % axial velocity at the previous rotor disk
+ prevRotVelAx = PrevOpRot.ElPerf.upstreamVelAx + PrevOpRot.ElPerf.indVelAx;
+
+ % axial downstream velocity for previous rotor
+ prevDownVelAx = PrevOpRot.ElPerf.upstreamVelAx + 2 * PrevOpRot.ElPerf.indVelAx;
+
+ % contraction variables
+ % radius ratio of of previous rotor contraction
+ vec_areaRatio = prevRotVelAx ./ prevDownVelAx;
+ % area ratio of previous rotor contraction
+ vec_radiusRatio = sqrt(vec_areaRatio);
+
+ vcRadius = PrevOpRot.Rot.radius * vec_radiusRatio(1, end); % vena contracta radius
+ vcRadIdx = find(self.Rot.Bl.y <= vcRadius, 1, 'last');
+
+ % tangential velocity at the previous rotor disk
+ prevDownVelTg = PrevOpRot.ElPerf.dQ ./ (PrevOpRot.Rot.Bl.y .* ...
+ PrevOpRot.ElPerf.massflow) .* (1 ./ vec_areaRatio);
+
+ % Mapping of velocities from initial upper rotor splitting to contracted splitting
+ % contracted previous rotor velocities at current rotor elements
+ yqueryAx = interp1(PrevOpRot.Rot.Bl.y .* vec_radiusRatio, prevDownVelAx, ...
+ self.Rot.Bl.y, 'linear', 0);
+ yqueryTg = interp1(PrevOpRot.Rot.Bl.y .* vec_radiusRatio, prevDownVelTg, ...
+ self.Rot.Bl.y, 'linear', 0);
+
+ % Update upstream axial velocity with previous rotor wake
+ self.upstreamVelAx(1, 1:vcRadIdx) = yqueryAx(1, 1:vcRadIdx);
+
+ % Update upstream tangential velocity with previous rotor wake
+ self.upstreamVelTg(1, 1:vcRadIdx) = yqueryTg(1, 1:vcRadIdx);
+
+ % Verification of mass flow conservation
+ % total flow mass at previous rotor level
+ debit_u2 = trapz(2 * pi * PrevOpRot.Op.Flow.rho * PrevOpRot.Rot.Bl.dy * ...
+ PrevOpRot.Rot.Bl.y .* prevRotVelAx);
+ % total flow mass at downstream of previous rotor - There is contraction of wake
+ area_contracted = PrevOpRot.Rot.Bl.dy * PrevOpRot.Rot.Bl.y .* vec_areaRatio;
+ debit_u3 = trapz(2 * pi * PrevOpRot.Op.Flow.rho * area_contracted .* prevDownVelAx);
+ % total flow mass at upstream of current rotor
+ debit_d1 = trapz(2 * pi * self.Op.Flow.rho * ...
+ self.Rot.Bl.dy * self.Rot.Bl.y(1, 1:vcRadIdx) .* ...
+ self.upstreamVelAx(1, 1:vcRadIdx));
+ % comparison of the flow masses in percentage
+ % Remove the below semicolons to display the results before running the code :
+ diff = 100 * (debit_u2 - debit_u3) / debit_u2;
+ diff2 = 100 * (debit_u3 - debit_d1) / debit_u3;
+
+end
diff --git a/src/classes/@OperRotor/OperRotor.m b/src/classes/@OperRotor/OperRotor.m
index 4c5842c..bf62161 100644
--- a/src/classes/@OperRotor/OperRotor.m
+++ b/src/classes/@OperRotor/OperRotor.m
@@ -64,6 +64,8 @@ classdef OperRotor < handle
cQ (1, 1) double % Rotor torque coefficient, [-]
cP (1, 1) double % Rotor power coefficient, [-]
+ totalmassflow (1, 1) double % Mass flow through the rotor, [kg/s]
+
nonDim (1, 2) char % Non-dimensionalization factor ('US', 'EU')
end
@@ -119,6 +121,10 @@ classdef OperRotor < handle
cP = calccoeff(self, 'power');
end
+ function totalmassflow = get.totalmassflow(self)
+ totalmassflow = trapz(self.ElPerf.massflow);
+ end
+
function nonDim = get.nonDim(self)
nonDim = self.nonDim_cached;
end
diff --git a/src/classes/@Result/Result.m b/src/classes/@Result/Result.m
index 5230a0e..e37170a 100644
--- a/src/classes/@Result/Result.m
+++ b/src/classes/@Result/Result.m
@@ -43,10 +43,10 @@ classdef Result < handle
% Base properties of the elements
properties
- leishman (1, :) OperRotor % Results obtained with the leishman solver
- indfact (1, :) OperRotor % Results obtained with the indFact solver
- indvel (1, :) OperRotor % Results obtained with the indVel solver
- stahlhut (1, :) OperRotor % Results obtained with the stalhlut solver
+ leishman (:, :) OperRotor % Results obtained with the leishman solver
+ indfact (:, :) OperRotor % Results obtained with the indFact solver
+ indvel (:, :) OperRotor % Results obtained with the indVel solver
+ stahlhut (:, :) OperRotor % Results obtained with the stalhlut solver
operPts table % Table of operating points for each individual result
ModExt (1, 1) struct % Modelling enstension options (losses, etc)
end
@@ -115,14 +115,17 @@ function resTable = maketable(self, prop)
resTable.Properties.VariableNames = propName;
- for iOp = 1:size(self.operPts, 1)
+ for iRotor = 1:size(self.indvel, 2)
for iSolv = 1:length(self.SOLVERS)
- solver = self.SOLVERS{iSolv};
+ for iOp = 1:size(self.operPts, 1)
- if ~isempty(self.(solver))
- resTable.(iSolv)(iOp) = self.(solver)(iOp).(prop);
- else
- resTable.(iSolv)(iOp) = NaN;
+ solver = self.SOLVERS{iSolv};
+
+ if ~isempty(self.(solver))
+ resTable.(iSolv)(iOp, iRotor) = self.(solver)(iOp, iRotor).(prop);
+ else
+ resTable.(iSolv)(iOp, iRotor) = NaN;
+ end
end
end
end
diff --git a/src/classes/@Result/filter.m b/src/classes/@Result/filter.m
index 436858f..3b87b5c 100644
--- a/src/classes/@Result/filter.m
+++ b/src/classes/@Result/filter.m
@@ -6,7 +6,7 @@ function [filtered, idx] = filter(self, alt, speed, rpm, coll)
% -----
%
% Syntax:
- % [filtered, idx] = Result.filter(alt, rpm, coll, speed) Return the Results of the four
+ % [filtered, idx] = Result.filter(alt, speed, rpm, coll) Return the Results of the four
% solvers for the operating point defined by `alt`, `rpm`, `coll`, `speed`.
%
% Inputs:
@@ -61,7 +61,7 @@ function [filtered, idx] = filter(self, alt, speed, rpm, coll)
for iSolv = 1:length(DEF.ALLOWED_SOLVERS)
thissolv = DEF.ALLOWED_SOLVERS{iSolv};
if ~isempty(self.(thissolv))
- filtered.(thissolv) = self.(thissolv)(1, idx);
+ filtered.(thissolv) = self.(thissolv)(idx, :);
else
filtered.(thissolv) = self.(thissolv);
end
diff --git a/src/configs/caradonna1981.m b/src/configs/caradonna1981.m
index e96cace..48410f2 100644
--- a/src/configs/caradonna1981.m
+++ b/src/configs/caradonna1981.m
@@ -115,6 +115,7 @@ Blade.nElem = 100; % Number of blade elements, [-]
% Rotor base position
Blade.hubPos = [0, 0, 0]; % Rotor center position (used for coaxial rotors), [m]
+Blade.spinDir = 1; % Rotor spin direction (used for coaxial rotors; 1 = cw, -1 = ccw)
% ==================================================================================================
% ======================================= Overwrites ===============================================
diff --git a/src/configs/knight1937.m b/src/configs/knight1937.m
index d28fce5..7e6ee22 100644
--- a/src/configs/knight1937.m
+++ b/src/configs/knight1937.m
@@ -110,3 +110,4 @@ Blade.nElem = 100; % Number of blade elements, [-]
% Rotor base position
Blade.hubPos = [0, 0, 0]; % Rotor center position (used for coaxial rotors), [m]
+Blade.spinDir = 1; % Rotor spin direction (used for coaxial rotors; 1 = cw, -1 = ccw)
diff --git a/src/configs/template.m b/src/configs/template.m
index dbd803e..4f6de7f 100644
--- a/src/configs/template.m
+++ b/src/configs/template.m
@@ -162,3 +162,4 @@ Blade.nElem = 100; % Number of blade elements, [-]
% Rotor base position
Blade.hubPos = [0, 0, 0]; % Rotor center position (used for coaxial rotors), [m]
+Blade.spinDir = 1; % Rotor spin direction (used for coaxial rotors; 1 = cw, -1 = ccw)
diff --git a/src/configs/templatecoax.m b/src/configs/templatecoax.m
index ccfa168..764d984 100644
--- a/src/configs/templatecoax.m
+++ b/src/configs/templatecoax.m
@@ -34,6 +34,7 @@ else
rotorShift = [-1, 0, 0];
end
Blade(2).hubPos = Blade(1).hubPos + rotorShift; % Shift second rotor w.r.t the application type
+Blade(2).spinDir = -Blade(1).spinDir; % Spin in opposite direction from first rotor
% ==================================================================================================
% ===================================== Operating points ===========================================
diff --git a/src/rotare.m b/src/rotare.m
index 0a80671..efdf786 100644
--- a/src/rotare.m
+++ b/src/rotare.m
@@ -95,7 +95,7 @@ function [Results] = rotare(configFile)
if Sim.Out.show3D
% TODO: Allow plot for multiple rotors
- Rot.plot('all', Sim.Out.hubType);
+ % Rot.plot('all', Sim.Out.hubType);
end
% (Following comment is so miss-hit metric checker stays happy)
diff --git a/src/solvers/bemt.m b/src/solvers/bemt.m
index 1f08f66..abb3deb 100644
--- a/src/solvers/bemt.m
+++ b/src/solvers/bemt.m
@@ -36,13 +36,13 @@ function bemt(OpRot, Mod)
% Issues: https://gitlab.uliege.be/rotare/rotare/-/issues
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % For each rotor, calculate the external velocity and solve the BEMT equations
+ % Always calculate first rotor as if isolated
+ bemtsinglerot(OpRot(1), Mod);
- for i = 1:length(OpRot)
-
- % TODO: Coaxial: Determine proper external velocity to pass to the rotor here
+ for i = 2:length(OpRot)
% Solves BEMT equations for the rotor using the true distribution of external velocity
+ OpRot(i).ElPerf.updateupstreamvel(OpRot(i - 1));
bemtsinglerot(OpRot(i), Mod);
end
diff --git a/src/solvers/indvel.m b/src/solvers/indvel.m
index 580da9a..e2b7592 100644
--- a/src/solvers/indvel.m
+++ b/src/solvers/indvel.m
@@ -119,6 +119,8 @@ function indvel(OpRot, Mod)
OpRot.ElPerf.alpha = alpha;
OpRot.ElPerf.indVelAx = v - OpRot.ElPerf.upstreamVelAx;
OpRot.ElPerf.indVelTg = uw / 2;
+ OpRot.ElPerf.massflow = dmdot;
+
end
function force = coeff2force(coeff, OpRot, speed)
--
GitLab