From dc1fc8f434235617564408368a0314ca46205e98 Mon Sep 17 00:00:00 2001
From: Thomas Lambert <t.lambert@uliege.be>
Date: Mon, 25 Sep 2023 17:33:56 +0200
Subject: [PATCH] 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 | 4 +-
src/classes/@ElemPerf/updateupstreamvel.m | 247 ++++++++++++++++------
src/utils/preproc/validateconfig.m | 2 +-
3 files changed, 188 insertions(+), 65 deletions(-)
diff --git a/README.md b/README.md
index 60a40a7..bfd21d9 100644
--- a/README.md
+++ b/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.
diff --git a/src/classes/@ElemPerf/updateupstreamvel.m b/src/classes/@ElemPerf/updateupstreamvel.m
index 778d6a9..00506a0 100644
--- a/src/classes/@ElemPerf/updateupstreamvel.m
+++ b/src/classes/@ElemPerf/updateupstreamvel.m
@@ -1,24 +1,23 @@
-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
diff --git a/src/utils/preproc/validateconfig.m b/src/utils/preproc/validateconfig.m
index 413ac30..4471240 100644
--- a/src/utils/preproc/validateconfig.m
+++ b/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
--
GitLab