Merge branch 'coax_LR_new' into 'main'

Add upstream velocities

See merge request rotare/rotare!3

src/classes/@ElemPerf/ElemPerf.m

@@ -57,6 +57,8 @@ classdef ElemPerf < handle
         Op (1, 1)  Oper  % Handle of the Operating point linked to this specific ElemPerf instance
 
         reynolds (1, :) double % Reynolds number, [-]
+        upstreamVelAx (1, :) double % Upstream axial velocity above rotor disk, [m/s]
+        upstreamVelTg (1, :) double % Upstream tangential velocity above rotor disk, [m/s]
 
         tgSpeed   (1, :) double % Tangential speed, [m/s]
         pitch     (1, :) double % Pitch (twist + collective), [rad]
@@ -117,6 +119,12 @@ classdef ElemPerf < handle
                 % Blade pitch
                 self.pitch = Rot.Bl.twist + Op.coll;
 
+                % Upstream axial velocity
+                self.upstreamVelAx = ones(size(self.pitch)) * Op.speed;
+
+                % Upstream tangential velocity
+                self.upstreamVelTg = zeros(size(self.pitch));
+
                 % In-plane velocities
                 self.tgSpeed = Op.omega .* Rot.Bl.y;
 

src/classes/@ElemPerf/calcforces.m

@@ -38,7 +38,7 @@ function calcforces(self)
 
     if ~isempty(self.cl)
 
-        axVel = self.Op.speed + self.indVelAx;
+        axVel = self.upstreamVelAx + self.indVelAx;
         tgVel = self.tgSpeed - self.indVelTg;
         relVel = sqrt(axVel.^2 + tgVel.^2);
 

src/classes/@ElemPerf/plotveltriangles.m

@@ -67,13 +67,13 @@ function plotveltriangles(self, nTriangles, varargin)
 
     % Velocity triangles
     % FIXME: Needs adaptation for coaxial
-    vAx_up = ones(size(self.tgSpeed)) * self.Op.speed;
+    vAx_up = ones(size(self.tgSpeed)) .* self.upstreamVelAx;
     vTg_up = self.tgSpeed;
 
-    vAx_rot = self.Op.speed + self.indVelAx;
+    vAx_rot = self.upstreamVelAx + self.indVelAx;
     vTg_rot = self.tgSpeed - self.indVelTg;
 
-    vAx_down = self.Op.speed + 2 * self.indVelAx;
+    vAx_down = self.upstreamVelAx + 2 * self.indVelAx;
     vTg_down = self.tgSpeed - 2 * self.indVelTg;
 
     % Determine sections to be plotted

src/classes/@OperRotor/OperRotor.m

@@ -100,7 +100,7 @@ classdef OperRotor < handle
         end
 
         function relSpeed = get.relTipSpeed(self)
-            relSpeed = sqrt(self.tgTipSpeed^2 + self.Op.speed^2);
+            relSpeed = sqrt(self.tgTipSpeed.^2 + self.upstreamVelAx.^2);
         end
 
         function advRat = get.advanceRatio(self)

src/solvers/indvel.m

@@ -49,7 +49,7 @@ function indvel(OpRot, Mod)
     % FSOLVE_OPTS = optimoptions('fsolve', 'Display', 'none');
 
     % DEFAULTS AND CONSTANTS
-    vw = ones(1, OpRot.Rot.Bl.nElem) * OpRot.Op.speed + 1; % Axial vel in the slipstream
+    vw = OpRot.ElPerf.upstreamVelAx + 1; % Axial vel in the slipstream
     uw = zeros(size(vw)); % Tangential vel in the slipstream
 
     % Solve iteratively
@@ -62,7 +62,7 @@ function indvel(OpRot, Mod)
         uw_old = uw;
 
         % Compute the velocity components at the propeller disk
-        v = (OpRot.Op.speed + vw) / 2;
+        v = (OpRot.ElPerf.upstreamVelAx + vw) / 2;
         u = OpRot.ElPerf.tgSpeed - uw / 2;
 
         % Local mass flow rate
@@ -93,7 +93,7 @@ function indvel(OpRot, Mod)
         dFu = OpRot.Rot.nBlades * (dL .* sin(phi) + dD .* cos(phi));
 
         % Final system to solve for v and u
-        vw = OpRot.Op.speed + dFa ./ dmdot ./ K_T;
+        vw = OpRot.ElPerf.upstreamVelAx + dFa ./ dmdot ./ K_T;
         uw = dFu ./ dmdot ./ K_P;
 
         % Use relaxation to faciliate convergence of nonlinear system
@@ -116,7 +116,7 @@ function indvel(OpRot, Mod)
     % Update values in ElemPerf
     OpRot.ElPerf.inflAngle = phi;
     OpRot.ElPerf.alpha = alpha;
-    OpRot.ElPerf.indVelAx = v - OpRot.Op.speed;
+    OpRot.ElPerf.indVelAx = v - OpRot.ElPerf.upstreamVelAx;
     OpRot.ElPerf.indVelTg = uw / 2;
 
 end