fix(solvers): use proper upstreamvel in all solvers (closes #12)

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
 

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

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>

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

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

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;

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;

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