Merge branch 'coax_LR_new' into 'dev_coax'

feat: add SMST coaxial model
See merge request rotare/rotare!6

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)

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

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

src/classes/@ElemPerf/updateupstreamvel.m

+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

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

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

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

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 ===============================================

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)

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)

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 ===========================================

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)

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

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)