diff --git a/src/classes/@ElemPerf/ElemPerf.m b/src/classes/@ElemPerf/ElemPerf.m
index 08f4fe3db6fbe553bdff9826fd6e86442ef06c22..c700032a2e7541daa1f93bc22e085275d738295e 100644
--- a/src/classes/@ElemPerf/ElemPerf.m
+++ b/src/classes/@ElemPerf/ElemPerf.m
@@ -227,7 +227,7 @@ classdef ElemPerf < handle
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
@@ -249,7 +249,7 @@ 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
@@ -271,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
diff --git a/src/classes/@OperRotor/OperRotor.m b/src/classes/@OperRotor/OperRotor.m
index bf621617b13e883867b8d6a36f37b27299b0d19f..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
diff --git a/src/classes/@Result/filter.m b/src/classes/@Result/filter.m
index 3b87b5c9e20ab6ec267a2b193217afd6a07db96e..c563ed008028425536d681f92757a053fd866b15 100644
--- a/src/classes/@Result/filter.m
+++ b/src/classes/@Result/filter.m
@@ -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>
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/templatecoax.m b/src/configs/templatecoax.m
index 5615acb23d50ee7f1d5d7a85711964d69f1e03b6..978dfcda4d1393630728421b8512221aa3afc07a 100644
--- a/src/configs/templatecoax.m
+++ b/src/configs/templatecoax.m
@@ -16,7 +16,7 @@
% Issues: https://gitlab.uliege.be/rotare/rotare/-/issues
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-N_ROTORS = 4; % Total number of rotors in coaxial configuration, [-]
+N_ROTORS = 2; % Total number of rotors in coaxial configuration, [-]
% ==================================================================================================
% ================================== Baseline configuration ========================================
@@ -26,9 +26,6 @@ template;
Sim.Out.show3D = false; % Disable 3D as it is not supported for coaxial (yet)
-% FIXME: Currently only indvel is supported for coaxial
-Mod.solvers = 'indvel'; % BEMT Solver ('leishman', 'indfact', 'indvel', 'stahlhut', 'all')
-
% ==================================================================================================
% ====================================== COAXIAL SYSTEM ===========================================
% ==================================================================================================
@@ -54,6 +51,6 @@ for i = 2:N_ROTORS
Blade(i).hubPos = Blade(i - 1).hubPos + rotorShift;
% =================================== Operating points =========================================
- Op.collective = [Op.collective; init_coll];
- Op.rpm = [Op.rpm; init_rpm];
+ Op.collective = [Op.collective; init_coll + 5];
+ Op.rpm = [Op.rpm; init_rpm * 1.5];
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/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