diff --git a/README.md b/README.md
index 60a40a70bf43f5014480409a7c43c4890b334653..bfd21d98097c1112f6a5362d9bf0ce65e813747f 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 19c54624b332570b371bbeaebe6240bcd1a89e14..423c88e1b46e8ff3a25886564bf98df03467f143 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 e30aafb034b3a8d549a52b1ea4c5b0c9fa8b7178..4c4d7dfef27d3c5c92162167ac82b45421989124 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 f6700dc2c52f3b09dcccde00f86f1546c1291731..6d417116a17c9c45638371f31df84591e01844bb 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 0000000000000000000000000000000000000000..00506a05bebcc424ef9a0535e781061c2e60f806
--- /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 4c5842c1d8185585f47250747927356ebd3add71..948fccc8951cd3b9e6cd7d9d1074c3e9828ec97a 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 5230a0ed9f4f4e13d9206e7c9aab9051bcacca10..0e71c09ab0971a88a42f19ce24807092f6c7bb54 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 436858f386967841f1254a20931c041a697af4a4..c563ed008028425536d681f92757a053fd866b15 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 0749dbee49d65704cdbe7b39c15d5183b84753b4..fc789ec1a009acb6fc52c97394a4da5a7e60060d 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 e96cacec84a4321630bff298023a2b8a82540f25..bc082d9173ee18fb59120ec54773c90b22655502 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 d28fce505d3331f63a8c25243a16c34d3fbc09f4..a328fa7c36bb43d6ec8bb5fab2295d5f1d05c59d 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 dbd803ec28c88c378ef3a1fc983671fbea5fcd71..6df35d0a42c169f7b51cafe69b8c65d1a95a7b91 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 ccfa168eb3ad1cd95e8a94d51016e2b139e4b5de..aae4f71b77149ea2c5c54ae7de10eacb2d9ba678 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 0a806711d1e0127109290091b4c88f6982a26c6a..7ef0451e2afc361db4579acf9f108a506fbf04f5 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 1f08f66ef2da7b200eb007bb2c9d32c42ed20aeb..8924ea437b70ab862021e0b3b5f92050ea72e4c6 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 11ac1eaebd9adcab4a0affc684f98c5daab04b9d..f635429b4788b97e0a7badf2f667ea57d0569c0b 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 580da9a8d6ca8a02a4a90f8ebbbf4a296f1765bb..a0def66b266631dbafc87a74cb5a572b649e55ad 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 54b69783b31fb333092d8d9bda57a48ce366613b..aa51f1eada3985c5213a10b75cb58f3c0ac9ac8a 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 15b8dce0850024808cf6a2eca1201c854011d064..a6997be81d2780ea3bc606f5fda700249a401258 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 48696340293827117629d22e62efbdcf681c5470..44712405c2996fe1a90884dd3f8e185640e49929 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