From 204b1295b4cb238e4a8c775ba3578601fec2d69d Mon Sep 17 00:00:00 2001
From: Rakotondratsimba Lyraie <liz_tp@proton.me>
Date: Wed, 26 Jul 2023 14:41:40 +0000
Subject: [PATCH] Add upstream velocities

---
 src/classes/@ElemPerf/ElemPerf.m         | 8 ++++++++
 src/classes/@ElemPerf/calcforces.m       | 2 +-
 src/classes/@ElemPerf/plotveltriangles.m | 6 +++---
 src/classes/@OperRotor/OperRotor.m       | 2 +-
 src/solvers/indvel.m                     | 8 ++++----
 5 files changed, 17 insertions(+), 9 deletions(-)

diff --git a/src/classes/@ElemPerf/ElemPerf.m b/src/classes/@ElemPerf/ElemPerf.m
index beaaacd..78b33cf 100644
--- a/src/classes/@ElemPerf/ElemPerf.m
+++ b/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;
 
diff --git a/src/classes/@ElemPerf/calcforces.m b/src/classes/@ElemPerf/calcforces.m
index 62300ea..f6700dc 100644
--- a/src/classes/@ElemPerf/calcforces.m
+++ b/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);
 
diff --git a/src/classes/@ElemPerf/plotveltriangles.m b/src/classes/@ElemPerf/plotveltriangles.m
index 243b4f8..bd5a558 100644
--- a/src/classes/@ElemPerf/plotveltriangles.m
+++ b/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
diff --git a/src/classes/@OperRotor/OperRotor.m b/src/classes/@OperRotor/OperRotor.m
index d749984..4c5842c 100644
--- a/src/classes/@OperRotor/OperRotor.m
+++ b/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)
diff --git a/src/solvers/indvel.m b/src/solvers/indvel.m
index 2f5c172..0a340c2 100644
--- a/src/solvers/indvel.m
+++ b/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
-- 
GitLab