diff --git a/src/classes/@ElemPerf/ElemPerf.m b/src/classes/@ElemPerf/ElemPerf.m
index beaaacd4ece9f3e5d957698fbebaefdf7d02f246..78b33cf5dcf62dff22bf80141212ec1449b6f309 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 62300eabd9ec94ba1c29598b5c6826a5c31bb86b..f6700dc2c52f3b09dcccde00f86f1546c1291731 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 243b4f839b26134be4535ecdc596e912ddabad14..bd5a558233982c1009868ad7a05c84ae82354cf6 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 d749984148068a3bb84c811fb422423a66c70f30..4c5842c1d8185585f47250747927356ebd3add71 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 2f5c172d52e85ed531e4727d29f29d5d429ecd04..0a340c28aea0ff7b94a07dbd818f3139670db5f0 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