feat(coax): add sst and mst coax models + refact

Complete refactoring of the updateupstreamvel function in order to
facilitate inclusion of new models.

Added Single StreamTube model and a more exact calculation of the
streamtube constraction using the MST model.

README.md

@@ -35,6 +35,7 @@ A more exhaustive list of features can be found in the [complete
 documentation][rotare-doc] that will be made available with version 1.0.0.
 
 ### Types of rotors
+
 - [x] Helicopters main or tail rotors
 - [x] Aircraft propellers
 - [ ] Wind/tidal turbines
@@ -128,7 +129,8 @@ contact me directly at tlambert@uliege.be.
 
 We would like to thank Theo Delvaux for his initial implementation of the
 oblique flight with the `indfact` solver. We would also like to thank Johan Le
-for his work on the implementation of coaxial rotors during his master thesis.
+and Lyraie Rakotondratsimba for their work on the implementation of coaxial
+rotors during his master thesis.
 
 Their excellent contributions facilitated considerably the final implementation
 of these features in Rotare.

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.
+function updateupstreamvel(self, PrevOpRot, Mod)
+    % UPDATEUPSTREAMVEL Update upstream velocities to account for previous rotor downwash.
+    %   This function updates the axial and tangential upstream velocity of the current elements by
+    %   taking into account the downwash of the previous rotor.
+    %
+    % Notes:
+    %   At the moment, only the Simplified Multiple Stream Tube (SMST) model is implemented. It
+    %   assumes:
+    %     - no pre-contraction caused by the second rotor
+    % -----
     %
     % Syntax:
-    %    ElPerf.updateupstreamvel(PrevOpRot) updates upstream velocities in an ElPerf
-    %    object according to the performance of the previous Operating Rotor `PrevOpRot`.
+    %   ElPerf.updateupstreamvel(PrevOpRot) updates upstream velocities in an ElemPerf object
+    %   according to the performance of the previous Operating Rotor `PrevOpRot`.
     %
     % Inputs:
-    %   - prevOpRot : previous Operating Rotor object
+    %   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
+    %   Updated values for upstreamVelAx and upstreamVelTg for the current ElemPerf object
     %
     % See also: ElemPerf.
     %
@@ -26,7 +25,8 @@ function updateupstreamvel(self, PrevOpRot)
 
     % ----------------------------------------------------------------------------------------------
     % (c) Copyright 2022-2023 University of Liege
-    % Author: Thomas Lambert <t.lambert@uliege.be>
+    % Author: Lyraie Rakotondratsimba
+    % Co-Author/Maintainer: Thomas Lambert <t.lambert@uliege.be>
     % ULiege - Aeroelasticity and Experimental Aerodynamics
     % MIT License
     % Repo: https://gitlab.uliege.be/rotare/rotare
@@ -34,52 +34,173 @@ function updateupstreamvel(self, PrevOpRot)
     % 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.Bl.y(end) * 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;
+    % Map the velocities according to the appropriate mode
+    switch Mod.coax.model
+        case 'none'
+            Mod.coax.contraType = 'none';
+            prevWake = calccontraction(PrevOpRot, Mod.coax);
+            mapvelocities(self, prevWake, length(self.Rot.Bl.y));
+        case 'sst'
+            sst(self, PrevOpRot, Mod.coax);
+        case 'mst'
+            mst(self, PrevOpRot, Mod.coax, false);
+        case 'smst'
+            mst(self, PrevOpRot, Mod.coax, true);
+    end
+
+end
+
+% ==================================== Helper functions ============================================
+
+function prevWake = calccontraction(PrevOpRot, ModCoax)
+    % CALCCONTRACTION Calculates the wake characteristics of the first rotor
+
+    % Velocity at the disk of previous rotor
+    prevRot.velAx = PrevOpRot.ElPerf.upstreamVelAx + PrevOpRot.ElPerf.indVelAx;
+    prevRot.velTg = PrevOpRot.ElPerf.indVelTg;
+    prevRot.upstreamVel = PrevOpRot.ElPerf.upstreamVelAx;
+    if strcmpi(ModCoax.model, 'sst')
+        prevRot.velAx = mean(prevRot.velAx);
+        prevRot.velTg = 0;
+        prevRot.upstreamVel = mean(PrevOpRot.ElPerf.upstreamVelAx);
+    end
+
+    % Velocities in the wake of previous rotor after contraction
+    switch ModCoax.contraType
+        case 'none'
+            prevWake.velAx = prevRot.velAx;
+            prevWake.velTg = prevRot.velTg;
+            prevWake.areaR = 1;
+        case 'farfield'
+            prevWake.velAx = 2 * prevRot.velAx - prevRot.upstreamVel;
+            prevWake.velTg = 2 * prevRot.velTg; % FIXME?
+            prevWake.areaR = prevRot.velAx ./ prevWake.velAx;
+            % prevWake.VelTg = PrevOpRot.ElPerf.dQ ./ (PrevOpRot.Rot.Bl.y .* ...
+            %                  PrevOpRot.ElPerf.massflow) .* (1 ./ areaRatio);
+        case 'value'
+            prevWake.areaR = Mod.coax.contraVal;
+            prevWake.velAx = prevRot.velAx ./ prevWake.areaR;
+            prevWake.velTg = prevRot.velTg ./ prevRot.velTg;
+    end
+
+    % This assumes annulus contraction, but not thickness contraction as well
+    % It must be recalculated in general MST model
+    prevWake.y = PrevOpRot.Rot.Bl.y .* sqrt(prevWake.areaR);
+    prevWake.dy = PrevOpRot.Rot.Bl.dy .* sqrt(prevWake.areaR);
+end
+
+function sst(self, PrevOpRot, ModCoax)
+    % SSTMODEL Implement the Simple Stream Tube model.
+    %   Only consider the mean axial velocity and the entire streamtube contraction
+
+    prevWake = calccontraction(PrevOpRot, ModCoax);
+    prevWake.spinDir = self.Rot.Bl.spinDir * PrevOpRot.Rot.Bl.spinDir;
+
+    % Resize properly and discard the swirl velocity
+    prevWake.velAx = ones(size(prevWake.y)) * prevWake.velAx;
+    prevWake.velTg = zeros(size(prevWake.y));
+
+    % Map velocities
+    vcRadIdx = calccradius(prevWake.y, self.Rot.Bl.y);
+    mapvelocities(self, prevWake, vcRadIdx);
+
+    % Check if massflows are conserved
+    massflowcheck(self, PrevOpRot, prevWake, vcRadIdx, ModCoax);
+
+end
+
+function mst(self, PrevOpRot, ModCoax, simple)
+    % MSTMODEL Implement the Multiple Stream Tube model.
+    %   Consider individual annulus contraction, as well as the effect on the swirl velocity
+
+    % Base calculation with only annulus contraction, no tube narrowing
+    prevWake = calccontraction(PrevOpRot, ModCoax);
+    prevWake.spinDir = self.Rot.Bl.spinDir * PrevOpRot.Rot.Bl.spinDir;
+
+    if ~simple
+        % Contract both annulus and tube thickness
+        yPrevExt = zeros(size(PrevOpRot.Rot.Bl.y)); % External radius of previous tube
+        yPrevExt(1) = PrevOpRot.Rot.Bl.y(1) - PrevOpRot.Rot.Bl.dy(1) / 2; % Hub
+
+        yw = zeros(size(PrevOpRot.Rot.Bl.y)); % Annulus center in wake
+        dyw = zeros(size(PrevOpRot.Rot.Bl.y)); % Annulus thickness in wake
+        for i = 1:numel(PrevOpRot.Rot.Bl.y)
+            if i > 1
+                yPrevExt(i) = yw(i - 1) + dyw(i - 1) / 2;
+            end
+
+            dyw(i) = -yPrevExt(i) + sqrt(yPrevExt(i)^2 + ...
+                                         2 * PrevOpRot.Rot.Bl.y(i) .* PrevOpRot.Rot.Bl.dy .* ...
+                                         prevWake.areaR(i));
+
+            yw(i) = yPrevExt(i) + dyw(i) / 2;
+        end
+        prevWake.y = yw;
+        prevWake.dy = dyw;
+
+    end
+
+    % Map velocities
+    vcRadIdx = calccradius(prevWake.y, self.Rot.Bl.y);
+    mapvelocities(self, prevWake, vcRadIdx);
+
+    % Check if massflows are conserved
+    massflowcheck(self, PrevOpRot, prevWake, vcRadIdx, ModCoax);
+end
+
+function vcIdxOnLower = calccradius(vcY, y)
+    % CALCVCRADIUS Calculate the radius of the vena contracta and its index on lower rotor
+
+    vcIdxOnLower = find(y <= vcY(end), 1, 'last');
+end
+
+function mapvelocities(self, prevWake, vcRadIdx)
+    % MAPVELOCITIES Map the velocities of the previous rotor wake onto the current rotor
+
+    self.upstreamVelAx(1:vcRadIdx) = interp1(prevWake.y, prevWake.velAx, ...
+                                             self.Rot.Bl.y(1:vcRadIdx), 'spline', 0);
+    self.upstreamVelTg(1:vcRadIdx) = interp1(prevWake.y, prevWake.spinDir * prevWake.velTg, ...
+                                             self.Rot.Bl.y(1:vcRadIdx), 'spline', 0);
+end
+
+function massflowcheck(self, PrevOpRot, prevWake, vcRadIdx, ModCoax)
+    % MASSFLOWCHECK Verify if the massflow is conserved after inteprolation in the vena contracta
+    %   The following three massflows should be roughly equal:
+    %     1. Mass flow accros entire previous rotor disk
+    %     2. Mass flow in wake after contraction of slipstream
+    %     3. Mass flow on lower rotor elements inside vena contracta after interpolation of wake
+
+    % 1. Entire previous rotor disk
+    if strcmpi(ModCoax.model, 'sst')
+        mf_prevRot = trapz(2 * pi * PrevOpRot.Rot.Bl.y .* PrevOpRot.Rot.Bl.dy .* ...
+                           PrevOpRot.Op.Flow.rho .* ...
+                           mean(PrevOpRot.ElPerf.upstreamVelAx + PrevOpRot.ElPerf.indVelAx));
+    else
+        mf_prevRot = PrevOpRot.totalmassflow;
+    end
+
+    % 2. Contracted wake
+    mf_prevWake = trapz(2 * pi * prevWake.y .* prevWake.dy .* PrevOpRot.Op.Flow.rho .* ...
+                        prevWake.velAx);
+    % 3. Lower rotor elements inside vena contracta
+    mf_curRot = trapz(2 * pi * self.Rot.Bl.y(1, 1:vcRadIdx) .* self.Rot.Bl.dy .* ...
+                      PrevOpRot.Op.Flow.rho .* self.upstreamVelAx(1:vcRadIdx));
+
+    % Error in mass flow conservation
+    err_prevRot = (mf_prevRot - mf_prevWake) / mf_prevRot;  % prevRot disk/wake
+    err_intrpCurr = (mf_prevWake - mf_curRot) / mf_prevWake; % prevRot wake/current vena contracta
+
+    %     close all;
+    %     figure
+    %     hold on
+    %     plot(PrevOpRot.Rot.Bl.y, PrevOpRot.ElPerf.upstreamVelAx + PrevOpRot.ElPerf.indVelAx)
+    %     plot(prevWake.y, prevWake.velAx)
+    %     plot(self.Rot.Bl.y, self.upstreamVelAx,'-o')
+    %     legend('upper', 'vc','lower')
+    %     pause
+
+    % Throw error if mass flow not conserved
+    assert(abs(err_prevRot) < 1e-12); % Should be 0 in theory
+    assert(abs(err_intrpCurr) < 2e-2); % Should be small
 
 end

src/utils/preproc/validateconfig.m

@@ -207,7 +207,7 @@ function Mod = checkmod(Mod, configFile, nRotors, DEF)
         if strcmp(Mod.coax.contraType, 'value')
             validateattributes(Mod.coax.contraVal, ...
                                {'double'}, ...
-                               {'scalar', 'positive'}, ...
+                               {'scalar', 'positive', '<=', 1}, ...
                                configFile, 'Mod.Num.convCrit');
         end
     end