From a0ada711e86286a601ccf09b38936b03dd12f7f8 Mon Sep 17 00:00:00 2001
From: Thomas Lambert <t.lambert@uliege.be>
Date: Mon, 25 Sep 2023 17:47:15 +0200
Subject: [PATCH] feat(coax): add coaxial models
Squashed commit of the following:
commit dc1fc8f434235617564408368a0314ca46205e98
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Mon Sep 25 17:33:56 2023 +0200
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.
commit 422dc4575c735a560e787fd70be73e46956f0192
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Fri Sep 22 14:53:42 2023 +0200
feat(coax,WIP): add options for contraction
commit 8ad2362531d4a211090223439c6d80d22da414d2
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Wed Sep 20 15:07:12 2023 +0200
fix(ElemPerf): issue with contraction
The contraction should be calculated based on the **last element**
position and not the actual rotor radius.
In some edge cases, an element of the aft rotor can fall in between the
contraction calculated from the previous rotor's last element and the one
calculated from its tip radius. When this happens, the velocity is
poorly interpolated and brought to 0 for this element on the aft rotor.
This not only skews the results, but also lead to issues with the
current implementation of the indfact solver, which only looks into
Op.speed~=0 and into the airspeed for each individual elements.
commit ccf3cfb40c04d0dc79c6c12254fa9c86ae80a56e
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Wed Sep 20 15:06:05 2023 +0200
fix(solvers): use proper upstreamvel in all solvers (closes #12)
commit eeb8cc2f655d06e2841a20f3c6fcb232bfd9aeb1
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Wed Sep 20 13:13:20 2023 +0200
fix(Result): display op pts for multi-rotors (closes #15)
commit c8525bcb8ebb8ab8f38819ee7d9e6ab832336e07
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Wed Sep 20 11:38:50 2023 +0200
feat(solvers, WIP): coax for all solvers
commit 67225be0a6c074aec849d2a26028729723be1cec
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Wed Sep 20 10:59:46 2023 +0200
refact(config): more generic coax template
commit bd6c21b5653aa35833066c1fa8e8dc4cec46d4b9
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Wed Sep 20 10:56:18 2023 +0200
refact(result): improve result tables
commit 4b24a542f0dadc77970239ddf18213b517fd52d5
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Wed Sep 20 09:48:46 2023 +0200
chore(bemt): add warning for coax solvers
commit ab5bc9f78c0723c85173f7eb34a1192ebe6a1652
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Tue Sep 19 16:50:30 2023 +0200
doc(Blade): spindir definition
commit d639096e05c4a10e94b699895d267ba72de7ca5a
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Tue Sep 19 16:46:31 2023 +0200
chore(template): minor adjustments
commit 76f8791de99310ed11c0a15f41b5029722ccbf3f
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Tue Sep 19 16:21:33 2023 +0200
chore(valid): add checks for spindir and hubpos
commit 510808d01574a9a46a73d6f72cbfca3e011af66d
Merge: 83427a0 af54595
Author: Thomas Lambert <t.lambert@uliege.be>
Date: Tue Sep 19 14:08:07 2023 +0000
Merge branch 'coax_LR_new' into 'dev_coax'
feat: add SMST coaxial model
See merge request rotare/rotare!6
commit af5459537b25c15bab8a284676422cd41091a554
Author: Rakotondratsimba Lyraie <liz_tp@proton.me>
Date: Tue Sep 19 14:08:07 2023 +0000
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
---
README.md | 4 +-
src/classes/@Blade/Blade.m | 21 ++-
src/classes/@ElemPerf/ElemPerf.m | 28 ++-
src/classes/@ElemPerf/calcforces.m | 2 +-
src/classes/@ElemPerf/updateupstreamvel.m | 206 ++++++++++++++++++++++
src/classes/@OperRotor/OperRotor.m | 8 +
src/classes/@Result/Result.m | 59 ++++---
src/classes/@Result/filter.m | 6 +-
src/classes/@Result/plotperf.m | 3 +-
src/configs/caradonna1981.m | 6 +
src/configs/knight1937.m | 6 +
src/configs/template.m | 8 +-
src/configs/templatecoax.m | 48 +++--
src/rotare.m | 5 +-
src/solvers/bemt.m | 26 ++-
src/solvers/indfact.m | 11 +-
src/solvers/indvel.m | 1 +
src/solvers/leishman.m | 6 +-
src/solvers/stahlhut.m | 2 +-
src/utils/preproc/validateconfig.m | 31 +++-
20 files changed, 397 insertions(+), 90 deletions(-)
create mode 100644 src/classes/@ElemPerf/updateupstreamvel.m
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/@Blade/Blade.m b/src/classes/@Blade/Blade.m
index 19c5462..423c88e 100644
--- a/src/classes/@Blade/Blade.m
+++ b/src/classes/@Blade/Blade.m
@@ -19,15 +19,16 @@ classdef Blade < handle
% -----
%
% Blade properties:
- % nElem - Number of elements, [-]
- % dy - Element span, [m]
- % y - Element absolute radial position, [m]
- % r - Element relative radial position (relative to the total radius), [-]
- % area - Element area, [m^2]
- % chord - Element chord, [m]
- % twist - Element twist (also called stagger angle), [rad]
- % iAf - Index of airfoil to use for each element, [-]
- % sol - Element solidity, [-]
+ % nElem - Number of elements, [-]
+ % dy - Element span, [m]
+ % y - Element absolute radial position, [m]
+ % r - Element relative radial position (relative to the total radius), [-]
+ % area - Element area, [m^2]
+ % chord - Element chord, [m]
+ % twist - Element twist (also called stagger angle), [rad]
+ % iAf - Index of airfoil to use for each element, [-]
+ % sol - Element solidity, [-]
+ % spinDir - Rotor spin direction (used for coaxial rotors; 1 = cw, -1 = ccw), [-]
%
% Blade methods:
%
@@ -44,6 +45,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 : Rotor spin direction (used for coaxial rotors; 1 = cw, -1 = ccw), [-]
%
% See also: Rotor, rotare, template, af_tools.Airfoil.
%
@@ -69,6 +71,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) = 1 % Rotor spin direction (used for coaxial rotors; 1 = cw, -1 = ccw), [-]
end
properties (SetAccess = private, Dependent)
diff --git a/src/classes/@ElemPerf/ElemPerf.m b/src/classes/@ElemPerf/ElemPerf.m
index e30aafb..4c4d7df 100644
--- a/src/classes/@ElemPerf/ElemPerf.m
+++ b/src/classes/@ElemPerf/ElemPerf.m
@@ -36,6 +36,8 @@ classdef ElemPerf < handle
%
% <a href="https:/gitlab.uliege.be/rotare/documentation">Complete documentation (online)</a>
+ % ----------------------------------------------------------------------------------------------
+ % FIXME: Documentation
% ----------------------------------------------------------------------------------------------
% Implementation:
% - alpha is implemented as a dependent property to prevent Matlab complaining when we set cl
@@ -85,6 +87,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 +99,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 +179,10 @@ classdef ElemPerf < handle
end
+ function massflow = get.massflow(self)
+ massflow = self.massflow_;
+ end
+
function indVelAx = get.indVelAx(self)
indVelAx = self.indVelAx_;
end
@@ -211,13 +218,21 @@ 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;
% Calculate inflow ratios
- self.inflowRat_ = (self.Op.speed + val) ./ (self.Op.omega * self.Rot.radius);
+ self.inflowRat_ = (self.upstreamVelAx + val) ./ (self.Op.omega * self.Rot.radius);
self.indInflowRat_ = val ./ (self.Op.omega * self.Rot.radius);
+
+ % Calculate massflow rate
+ axialVel = self.upstreamVelAx + val;
+ self.massflow_ = 2 * pi * self.Op.Flow.rho * self.Rot.Bl.y * self.Rot.Bl.dy .* axialVel;
end
function set.indVelTg(self, val)
@@ -234,8 +249,12 @@ classdef ElemPerf < handle
self.inflowRat_ = val;
% Calculate induced velocity and ratio
- self.indVelAx_ = val .* self.Op.omega * self.Rot.radius - self.Op.speed;
+ self.indVelAx_ = val .* self.Op.omega * self.Rot.radius - self.upstreamVelAx;
self.indInflowRat_ = self.indVelAx ./ (self.Op.omega * self.Rot.radius);
+
+ % Calculate massflow rate
+ axialVel = self.upstreamVelAx + self.indVelAx_;
+ self.massflow_ = 2 * pi * self.Op.Flow.rho * self.Rot.Bl.y * self.Rot.Bl.dy .* axialVel;
end
function set.swirlRat(self, val)
@@ -252,7 +271,7 @@ classdef ElemPerf < handle
self.indInflowRat_ = val;
% Calculate inflow ratios
- self.inflowRat_ = (self.Op.speed) ./ (self.Op.omega * self.Rot.radius) + val;
+ self.inflowRat_ = (self.upstreamVelAx) ./ (self.Op.omega * self.Rot.radius) + val;
self.indVelAx_ = val .* (self.Op.omega * self.Rot.radius);
end
@@ -271,6 +290,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, Mod)
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..00506a0
--- /dev/null
+++ b/src/classes/@ElemPerf/updateupstreamvel.m
@@ -0,0 +1,206 @@
+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 ElemPerf object
+ % according to the performance of the previous Operating Rotor `PrevOpRot`.
+ %
+ % Inputs:
+ % prevOpRot : previous Operating Rotor object
+ %
+ % Outputs:
+ % Updated values for upstreamVelAx and upstreamVelTg for the current ElemPerf object
+ %
+ % See also: ElemPerf.
+ %
+ % <a href="https:/gitlab.uliege.be/rotare/documentation">Complete documentation (online)</a>
+
+ % ----------------------------------------------------------------------------------------------
+ % (c) Copyright 2022-2023 University of Liege
+ % 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
+ % Docs: https://gitlab.uliege.be/rotare/documentation
+ % Issues: https://gitlab.uliege.be/rotare/rotare/-/issues
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ % 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/classes/@OperRotor/OperRotor.m b/src/classes/@OperRotor/OperRotor.m
index 4c5842c..948fccc 100644
--- a/src/classes/@OperRotor/OperRotor.m
+++ b/src/classes/@OperRotor/OperRotor.m
@@ -14,6 +14,8 @@ classdef OperRotor < handle
% ElPerf - Performance of the elements (angles, velocities, forces, etc)
% tgTipSpeed - Tangential tip speed, [m/s]
% relTipSpeed - Relative tip speed, [m/s]
+ % totalmassflow - Mass flow through the rotor, [kg/s]
+
%
% OperRotor methods:
% calcperf - Calculate operating rotor performance
@@ -64,6 +66,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 +123,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..0e71c09 100644
--- a/src/classes/@Result/Result.m
+++ b/src/classes/@Result/Result.m
@@ -5,7 +5,7 @@ classdef Result < handle
%
% Notes:
% - Almost all properties are private to prevent accidental overwriting.
- %
+ % - Most properties are structured as [nRotor, nOperPoint]
% -----
%
% Result properties:
@@ -14,11 +14,16 @@ classdef Result < handle
% indvel - Results obtained with the indVel solver
% stahlhut - Results obtained with the stalhlut solver
% operPts - Table with the operating points for each individual result
- % ModExt - Modelling enstension options (losses, etc)
- %
+ % ModExt - Modelling extension options (losses, etc)
+ % cT - Table with all Thrust coefficients
+ % cP - Table with all Power coefficients
+ % cQ - Table with all Torque coefficients
%
% Result methods:
- % save - Save the Result object to file
+ % save - Save the Result object to file
+ % summarize - Output a summary table with the operating points and the performance
+ % filter - Filter for a specific operating point
+ % plotperf - Plot and compare rotor performances
%
% Result constructor:
% Result = Result() creates an empty object.
@@ -43,10 +48,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
@@ -73,9 +78,16 @@ classdef Result < handle
% ---------------------------------------------
% Set methods
function set.operPts(self, val)
- % Note: A bug in Matlab prevents the user of variable names when size checks are added
- % in the property definition. https://stackoverflow.com/q/48423003
+ % Note: A bug in Matlab prevents the use of variable names when size checks are added
+ % in the property definition. See: https://stackoverflow.com/q/48423003
self.operPts = val;
+
+ % Properly name and group operating points for multi rotors
+ % Must merge in that order to have correct indexes
+ nRots = (size(val, 2) - 2) / 2;
+ self.operPts = mergevars(self.operPts, 3 + nRots:size(val, 2));
+ self.operPts = mergevars(self.operPts, 3:3 + nRots - 1);
+
self.operPts.Properties.VariableNames = {'altitude', 'speed', 'rpm', 'collective'};
end
@@ -105,26 +117,21 @@ end
function resTable = maketable(self, prop)
% MAKETABLE Gather properties from all solvers and all operating points in a table
- resTable = table('Size', [size(self.operPts, 1), 4], ...
- 'VariableTypes', {'double', 'double', 'double', 'double'});
+ resTable = table('Size', [size(self.operPts, 1), 0]);
- propName = cell(1, 4);
- for i = 1:length(self.SOLVERS)
- propName{i} = [prop, '_', self.SOLVERS{i}];
- end
+ for iSolv = 1:length(self.SOLVERS)
+ thissolver = self.SOLVERS{iSolv};
+ tmpTable = [];
- resTable.Properties.VariableNames = propName;
-
- for iOp = 1:size(self.operPts, 1)
- for iSolv = 1:length(self.SOLVERS)
- solver = self.SOLVERS{iSolv};
-
- if ~isempty(self.(solver))
- resTable.(iSolv)(iOp) = self.(solver)(iOp).(prop);
- else
- resTable.(iSolv)(iOp) = NaN;
+ for iRotor = size(self.(thissolver), 2):-1:1
+ for iOp = size(self.operPts, 1):-1:1
+ tmpTable(iOp, iRotor) = self.(thissolver)(iOp, iRotor).(prop);
end
end
+
+ if ~isempty(tmpTable)
+ resTable = addvars(resTable, tmpTable, 'NewVariableNames', [prop, '_', thissolver]);
+ end
end
end
diff --git a/src/classes/@Result/filter.m b/src/classes/@Result/filter.m
index 436858f..c563ed0 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:
@@ -27,6 +27,8 @@ function [filtered, idx] = filter(self, alt, speed, rpm, coll)
%
% <a href="https://gitlab.uliege.be/rotare/documentation">Complete documentation (online)</a>
+ % ----------------------------------------------------------------------------------------------
+ % FIXME: Needs adjustments for multiple rotors
% ----------------------------------------------------------------------------------------------
% (c) Copyright 2022-2023 University of Liege
% Author: Thomas Lambert <t.lambert@uliege.be>
@@ -61,7 +63,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/classes/@Result/plotperf.m b/src/classes/@Result/plotperf.m
index 0749dbe..fc789ec 100644
--- a/src/classes/@Result/plotperf.m
+++ b/src/classes/@Result/plotperf.m
@@ -33,6 +33,8 @@ function plotperf(self, varargin)
%
% <a href="https://gitlab.uliege.be/rotare/documentation">Complete documentation (online)</a>
+ % ----------------------------------------------------------------------------------------------
+ % FIXME: Needs adjustments for multiple rotors
% ----------------------------------------------------------------------------------------------
% (c) Copyright 2022-2023 University of Liege
% Author: Thomas Lambert <t.lambert@uliege.be>
@@ -72,7 +74,6 @@ function plotperf(self, varargin)
operVect = self.operPts(idx, :).(oper);
end
dataVect = vectorize(filtered.(thissolv), perf);
-
plot(operVect, dataVect, lineSpec.(thissolv){:}, 'DisplayName', thissolv);
hold on;
end
diff --git a/src/configs/caradonna1981.m b/src/configs/caradonna1981.m
index e96cace..bc082d9 100644
--- a/src/configs/caradonna1981.m
+++ b/src/configs/caradonna1981.m
@@ -53,6 +53,11 @@ Mod.solvers = {'all'}; % BEMT Solver ('leishman', 'indfact', 'indvel', 'stahlh
% Extensions/corrections
Mod.Ext.losses = 'all'; % Include losses using Prandtl formula ('none', 'hub', 'tip', 'both')
+% [COAXIAL ONLY] Contraction model
+Mod.coax.model = 'mst'; % Type of contraction ('none', 'sst', 'smst', 'mst')
+Mod.coax.contraType = 'farfield'; % (opt) Type of contraction ('farfield', 'value')
+Mod.coax.contraVal = 1 / 2; % (opt) Contraction ratio between upper disk and lower rotor, [-]
+
% Numerical parameters
Mod.Num.convCrit = 1e-4; % Convergence criterion value,
Mod.Num.maxIter = 500; % Maximum number of iterations for convergence calculations
@@ -115,6 +120,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..a328fa7 100644
--- a/src/configs/knight1937.m
+++ b/src/configs/knight1937.m
@@ -48,6 +48,11 @@ Mod.solvers = {'indvel'}; % BEMT Solver ('leishman', 'indfact', 'indvel', 'sta
% Extensions/corrections
Mod.Ext.losses = 'all'; % Include losses using Prandtl formula ('none', 'hub', 'tip', 'both')
+% [COAXIAL ONLY] Contraction model
+Mod.coax.model = 'mst'; % Type of contraction ('none', 'sst', 'smst', 'mst')
+Mod.coax.contraType = 'farfield'; % (opt) Type of contraction ('farfield', 'value')
+Mod.coax.contraVal = 1 / 2; % (opt) Contraction ratio between upper disk and lower rotor, [-]
+
% Numerical parameters
Mod.Num.convCrit = 1e-4; % Convergence criterion value,
Mod.Num.maxIter = 500; % Maximum number of iterations for convergence calculations
@@ -110,3 +115,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..6df35d0 100644
--- a/src/configs/template.m
+++ b/src/configs/template.m
@@ -63,7 +63,7 @@ Sim.Save.filename = 'template'; % File name of the saved result
% Outputs
Sim.Out.showPlots = true; % Show all plots (forces, angles, speed, ...)
-Sim.Out.show3D = false; % Show the 3D view of the whole rotor
+Sim.Out.show3D = true; % Show the 3D view of the whole rotor
Sim.Out.hubType = 'tangent_ogive'; % Hub (nose cone) type on the 3D view (see docs for list)
Sim.Out.console = true; % Print the final results in console
Sim.Out.verbosity = 'min'; % Verbosity level of the console output ('min', 'all')
@@ -85,6 +85,11 @@ Mod.solvers = 'all'; % BEMT Solver ('leishman', 'indfact', 'indvel', 'stahlhut'
% Extensions/corrections
Mod.Ext.losses = 'none'; % Include losses using Prandtl formula ('none', 'hub', 'tip', 'both')
+% [COAXIAL ONLY] Contraction model
+Mod.coax.model = 'mst'; % Type of contraction ('none', 'sst', 'smst', 'mst')
+Mod.coax.contraType = 'farfield'; % Type of contraction ('farfield', 'value')
+Mod.coax.contraVal = 1 / 2; % (opt) Contraction ratio between upper disk and lower rotor, [-]
+
% Numerical parameters
Mod.Num.convCrit = 1e-4; % Convergence criterion value,
Mod.Num.maxIter = 500; % Maximum number of iterations for convergence calculations
@@ -162,3 +167,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..aae4f71 100644
--- a/src/configs/templatecoax.m
+++ b/src/configs/templatecoax.m
@@ -1,5 +1,5 @@
% TEMPLATECOAX Template configuration file for Coaxial rotors simulations with Rotare.
-% This file extends the regular template configuration by adding a second rotor geometry.
+% This file extends the regular template configuration by adding N-1 coaxial rotors.
% -----
%
% See also: template, rotare, validconfig.
@@ -16,30 +16,42 @@
% Issues: https://gitlab.uliege.be/rotare/rotare/-/issues
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+N_ROTORS = 2; % Total number of rotors in coaxial configuration, [-]
+
% ==================================================================================================
% ================================== Baseline configuration ========================================
% ==================================================================================================
-template; % Just load exisiting template
+% Just load exisiting template
+template;
+
+Sim.Out.show3D = false; % Disable 3D as it is not supported for coaxial (yet)
% ==================================================================================================
-% ================================= Blade and rotor geometry =======================================
+% ====================================== COAXIAL SYSTEM ===========================================
% ==================================================================================================
-% Adds a second element to the blade structure to represent the second rotor
+% Adds a as many rotors as required, based on the initial blade geometry
-Blade(2) = Blade(1); % Make second rotor identical to first one
+% First rotor settings for collective and RPM
+init_coll = Op.collective;
+init_rpm = Op.rpm;
+Blade(N_ROTORS - 1) = Blade(1); % Pre-alloc
-if strcmp(Sim.Misc.appli, 'heli')
- rotorShift = [0, 0, -1];
-else
- rotorShift = [-1, 0, 0];
-end
-Blade(2).hubPos = Blade(1).hubPos + rotorShift; % Shift second rotor w.r.t the application type
+for i = 2:N_ROTORS
-% ==================================================================================================
-% ===================================== Operating points ===========================================
-% ==================================================================================================
+ % ================================ Blade and rotor geometry ====================================
+ Blade(i) = Blade(1); % Make second rotor identical to first one
+ Blade(i).spinDir = -1 * Blade(i - 1).spinDir; % Invert spin direction for each rotor
+ Blade(i).nElem = round(Blade(i - 1).nElem * 1);
+
+ % Shift second rotor w.r.t the application type
+ if strcmp(Sim.Misc.appli, 'heli')
+ rotorShift = [0, 0, -1] * Blade(i - 1).radius(end);
+ else
+ rotorShift = [-1, 0, 0] * Blade(i - 1).radius(end);
+ end
+ Blade(i).hubPos = Blade(i - 1).hubPos + rotorShift;
-% Operating points must be set for the two rotors
-% We simply extend the collective and rpm fields by adding values for the new rotor
-Op.collective = [Op.collective; 1, 4, 7];
-Op.rpm = [Op.rpm; 900, 1200];
+ % =================================== Operating points =========================================
+ Op.collective = [Op.collective; init_coll + 5];
+ Op.rpm = [Op.rpm; init_rpm * 1.5];
+end
diff --git a/src/rotare.m b/src/rotare.m
index 0a80671..7ef0451 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(1).plot('all', Sim.Out.hubType);
end
% (Following comment is so miss-hit metric checker stays happy)
@@ -141,7 +141,8 @@ function [Results] = rotare(configFile)
% Table that identifies the operating point for future reuse
operPoints(iOperPoint, :) = [Uop.altitude(iAlt), Uop.speed(iSpeed), ...
- Uop.rpm(i, iRpm), Uop.collective(i, iColl)];
+ [Uop.rpm(:, iRpm)'], ...
+ [Uop.collective(:, iColl)']];
iOperPoint = iOperPoint + 1;
diff --git a/src/solvers/bemt.m b/src/solvers/bemt.m
index 1f08f66..8924ea4 100644
--- a/src/solvers/bemt.m
+++ b/src/solvers/bemt.m
@@ -9,22 +9,21 @@ function bemt(OpRot, Mod)
% -----
%
% Syntax:
- % Rot = BEMT(Rot, Mod, Flow)
+ % BEMT(OpRot, Mod) Calculates the performances of each OpRot object with the Blade Element
+ % Momentum Theory, using the solvers and models defined in Mod
%
% Inputs:
- % Rot : Rotor object fully defined
- % Flow : Flow variables
- % Mod : Model and numerical parameters
+ % OpRot : Operating rotor object(s) fully defined
+ % Mod : Model and numerical parameters
%
% Output:
- % Rot : Rotor object with calculated performances
+ % OperRotor object(s) with calculated performances
%
- % See also: rotare, leishman, stahlhut, propsolv, turbsolv.
+ % See also: rotare, leishman, indfact, indvel, stahlhut.
%
% <a href="https://gitlab.uliege.be/rotare/documentation">Complete documentation (online)</a>
% ----------------------------------------------------------------------------------------------
- % TODO: Implement coaxial rotors
% TODO: Implement oblique flows
% ----------------------------------------------------------------------------------------------
% (c) Copyright 2022-2023 University of Liege
@@ -36,15 +35,12 @@ function bemt(OpRot, Mod)
% Issues: https://gitlab.uliege.be/rotare/rotare/-/issues
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % For each rotor, calculate the external velocity and solve the BEMT equations
-
- for i = 1:length(OpRot)
-
- % TODO: Coaxial: Determine proper external velocity to pass to the rotor here
-
- % Solves BEMT equations for the rotor using the true distribution of external velocity
- bemtsinglerot(OpRot(i), Mod);
+ % First rotor always calculated as if isolated
+ bemtsinglerot(OpRot(1), Mod);
+ for i = 2:length(OpRot)
+ OpRot(i).ElPerf.updateupstreamvel(OpRot(i - 1), Mod); % Update upstream velocity
+ bemtsinglerot(OpRot(i), Mod); % Calculate BEMT as usual, with new upstream vel distribution
end
end
diff --git a/src/solvers/indfact.m b/src/solvers/indfact.m
index 11ac1ea..f635429 100644
--- a/src/solvers/indfact.m
+++ b/src/solvers/indfact.m
@@ -54,12 +54,12 @@ function indfact(OpRot, Mod)
INIT_VI = 1; % Initial guess for the axial induced velocity (for when Airspeed=0), [m/s]
% Abbreviations
- airspeed = OpRot.Op.speed;
+ airspeed = OpRot.ElPerf.upstreamVelAx;
% Initial guesses
if OpRot.Op.speed ~= 0
a = INIT_AX_FACT * ones(1, OpRot.Rot.Bl.nElem);
- v_ax = (1 + a) * airspeed;
+ v_ax = (1 + a) .* airspeed;
else
v_ax = INIT_VI * ones(1, OpRot.Rot.Bl.nElem);
end
@@ -93,7 +93,7 @@ function indfact(OpRot, Mod)
K_P = 1 - (1 - lossFact) .* sin(phi);
% Analytical solution to the momentum and blade element equations.
- if airspeed == 0
+ if OpRot.Op.speed == 0
v_ax = sqrt(((relVel.^2 * OpRot.Rot.nBlades .* OpRot.Rot.Bl.chord) .* ...
(cl .* cos(phi) - cd .* sin(phi))) ./ ...
(8 * pi * OpRot.Rot.Bl.y .* K_T));
@@ -109,8 +109,11 @@ function indfact(OpRot, Mod)
b = ((relVel.^2 * OpRot.Rot.nBlades .* OpRot.Rot.Bl.chord) .* ...
(cl .* sin(phi) + cd .* cos(phi))) ./ ...
(8 * pi * OpRot.Rot.Bl.y.^2 .* airspeed .* (1 + a) .* OpRot.Op.omega .* K_P);
- v_ax = (1 + a) * airspeed;
+ v_ax = (1 + a) .* airspeed;
end
+ % FIXME: We should pay attention for the case where individual element has a 0 airspeed due
+ % to upstreamvelocity (highly improbable). In such situation, this should be recomputed
+ % separately to prevent issues. -> Refactor the coefficient calculation into a function
% Use relaxation to faciliate convergence of nonlinear system
v_ax = v_ax_old * (1 - Mod.Num.relax) + v_ax * Mod.Num.relax;
diff --git a/src/solvers/indvel.m b/src/solvers/indvel.m
index 580da9a..a0def66 100644
--- a/src/solvers/indvel.m
+++ b/src/solvers/indvel.m
@@ -119,6 +119,7 @@ function indvel(OpRot, Mod)
OpRot.ElPerf.alpha = alpha;
OpRot.ElPerf.indVelAx = v - OpRot.ElPerf.upstreamVelAx;
OpRot.ElPerf.indVelTg = uw / 2;
+
end
function force = coeff2force(coeff, OpRot, speed)
diff --git a/src/solvers/leishman.m b/src/solvers/leishman.m
index 54b6978..aa51f1e 100644
--- a/src/solvers/leishman.m
+++ b/src/solvers/leishman.m
@@ -69,7 +69,7 @@ function leishman(OpRot, Mod)
lambda_old = lambda;
% Update Reynolds and then clSlope w.r.t. the new inflow velocity
- axialVel = OpRot.Op.speed + lambda * OpRot.tgTipSpeed;
+ axialVel = OpRot.ElPerf.upstreamVelAx + 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);
@@ -98,7 +98,7 @@ function leishman(OpRot, Mod)
% Update values in ElemPerf
OpRot.ElPerf.alpha = alpha;
OpRot.ElPerf.inflAngle = phi;
- OpRot.ElPerf.indVelAx = lambda * OpRot.tgTipSpeed - OpRot.Op.speed;
+ OpRot.ElPerf.indVelAx = lambda * OpRot.tgTipSpeed - OpRot.ElPerf.upstreamVelAx;
OpRot.ElPerf.indVelTg = OpRot.ElPerf.tgSpeed - xi * OpRot.tgTipSpeed;
end
@@ -113,7 +113,7 @@ function [lambda, xi, phi, alpha] = leishmaneq(OpRot, clSlope, lossFact)
% alpha : Angle of attack, [rad]
% Abbreviations
- lambda_c = OpRot.Op.speed / OpRot.tgTipSpeed;
+ lambda_c = OpRot.ElPerf.upstreamVelAx / OpRot.tgTipSpeed;
sol = OpRot.Rot.solidity; % [FIXME], normally it should be local solidity (see Carroll)
r = OpRot.Rot.Bl.r;
pitch = OpRot.ElPerf.truePitch;
diff --git a/src/solvers/stahlhut.m b/src/solvers/stahlhut.m
index 15b8dce..a6997be 100644
--- a/src/solvers/stahlhut.m
+++ b/src/solvers/stahlhut.m
@@ -142,7 +142,7 @@ function [gphi, b1phi, b2phi] = stahlhuteq(i, phi, OpRot, lossType)
reynolds = OpRot.ElPerf.reynolds(i); % Use approximate Reynolds by default
tgSpeed = OpRot.ElPerf.tgSpeed(i);
pitch = OpRot.ElPerf.truePitch(i);
- airspeed = OpRot.Op.speed;
+ airspeed = OpRot.ElPerf.upstreamVelAx(i);
% Recalculate Reynolds based on actual inflow velocity
% The first time we calculate it for a given element, we use the approximate value that
diff --git a/src/utils/preproc/validateconfig.m b/src/utils/preproc/validateconfig.m
index 4869634..4471240 100644
--- a/src/utils/preproc/validateconfig.m
+++ b/src/utils/preproc/validateconfig.m
@@ -74,6 +74,8 @@ function [Sim, Mod, Flow, Op, Airfoil, Blade] = validateconfig(configFile)
DEF.NONDIM = {'US', 'EU'};
DEF.VERBOSITY = {'min', 'all'};
DEF.SOLVERS = {'leishman', 'indfact', 'indvel', 'stahlhut', 'all'};
+ DEF.CONTRACTION_TYPE = {'farfield', 'value'};
+ DEF.CONTRACTION_MODEL = {'none', 'sst', 'smst', 'mst'};
DEF.FLUID = {'air', 'seawater', 'freshwater'};
DEF.LOSSES = {'none', 'hub', 'tip', 'both', 'all'};
DEF.POLARS = {'polynomial', 'file'};
@@ -88,12 +90,12 @@ function [Sim, Mod, Flow, Op, Airfoil, Blade] = validateconfig(configFile)
% Input checks
Sim = checksim(Sim, configFile, DEF); % Simulation parameters
- Mod = checkmod(Mod, configFile, DEF); % Models, solvers, etc.
Flow.fluid = validatestring(Flow.fluid, DEF.FLUID, configFile, 'Flow.fluid'); % Fluid
Airfoil = checkairfoil(Airfoil, configFile, DEF); % Airfoils data
Blade = checkblade(Blade, configFile, length(Airfoil), DEF); % Blade
Op = checkop(Op, configFile, DEF, numel(Blade), Flow.fluid); % Operating points
+ Mod = checkmod(Mod, configFile, numel(Blade), DEF); % Models, solvers, etc.
% ===========================================
% Extra warnings
@@ -169,7 +171,7 @@ end
% ==================================================================================================
% ==================================== Models and solvers ==========================================
% ==================================================================================================
-function Mod = checkmod(Mod, configFile, DEF)
+function Mod = checkmod(Mod, configFile, nRotors, DEF)
% Solvers
Mod.solvers = cellstr(Mod.solvers);
@@ -195,6 +197,21 @@ function Mod = checkmod(Mod, configFile, DEF)
{'scalar', 'positive'}, ...
configFile, 'Mod.Num.maxIter');
+ if nRotors > 1
+ Mod.coax.model = validatestring(Mod.coax.model, DEF.CONTRACTION_MODEL, ...
+ configFile, 'Mod.coax.model');
+ if ~strcmpi(Mod.coax.model, 'none')
+ Mod.coax.contraType = validatestring(Mod.coax.contraType, DEF.CONTRACTION_TYPE, ...
+ configFile, 'Mod.coax.contraType');
+ end
+ if strcmp(Mod.coax.contraType, 'value')
+ validateattributes(Mod.coax.contraVal, ...
+ {'double'}, ...
+ {'scalar', 'positive', '<=', 1}, ...
+ configFile, 'Mod.Num.convCrit');
+ end
+ end
+
end
% ==================================================================================================
@@ -345,6 +362,16 @@ function Blade = checkblade(Blade, configFile, nAirfoils, DEF)
'<=', nAirfoils}, ...
configFile, 'Blade.iAirfoil');
+ validateattributes(Blade(i).hubPos, ...
+ {'numeric'}, ...
+ {'vector', 'finite', 'nonnan', 'numel', 3}, ...
+ configFile, 'Blade.hubPos');
+ validateattributes(Blade(i).spinDir, ...
+ {'numeric'}, ...
+ {'scalar', 'integer', 'nonzero', ...
+ '<=', 1, '>=', -1}, ...
+ configFile, 'Blade.spinDir');
+
end
% Discretization
--
GitLab