diff --git a/src/classes/@ElemPerf/ElemPerf.m b/src/classes/@ElemPerf/ElemPerf.m
index d17cb070c33dfa316f7191ef7db9b6037e000668..e30aafb034b3a8d549a52b1ea4c5b0c9fa8b7178 100644
--- a/src/classes/@ElemPerf/ElemPerf.m
+++ b/src/classes/@ElemPerf/ElemPerf.m
@@ -160,15 +160,20 @@ classdef ElemPerf < handle
         end
 
         function reynolds = get.reynolds(self)
-            % If Reynolds in unknown, we use an approximate solution
-            % neglecting induced velocities.
-            if isempty(self.reynolds_)
-                relVel = sqrt(self.tgSpeed.^2 + self.upstreamVelAx.^2);
-                reynolds = Flow.reynolds(relVel, self.Rot.Bl.chord, ...
-                                         self.Op.Flow.mu, self.Op.Flow.rho);
-            else
-                reynolds = self.reynolds_;
+            % Always recalculate Reynolds based on the most accurate data
+            velAx = self.upstreamVelAx;
+            velTg = self.tgSpeed;
+
+            if ~isempty(self.indVelAx)
+                velAx = velAx + self.indVelAx;
+            end
+            if ~isempty(self.indVelTg)
+                velTg =  velTg - self.indVelTg;
             end
+            relVel = sqrt(velAx.^2 + velTg.^2);
+            reynolds = Flow.reynolds(relVel, self.Rot.Bl.chord, ...
+                                     self.Op.Flow.mu, self.Op.Flow.rho);
+
         end
 
         function indVelAx = get.indVelAx(self)
@@ -206,10 +211,6 @@ classdef ElemPerf < handle
             [self.cl, self.cd] = getclcd(self, val);
         end
 
-        function set.reynolds(self, val)
-            self.reynolds_ = val;
-        end
-
         function set.indVelAx(self, val)
             % Cache value
             self.indVelAx_ = val;
diff --git a/src/solvers/indfact.m b/src/solvers/indfact.m
index 54ed6d400804e9e8ea036f367c6b6cd2bdaa250c..11ac1eaebd9adcab4a0affc684f98c5daab04b9d 100644
--- a/src/solvers/indfact.m
+++ b/src/solvers/indfact.m
@@ -134,6 +134,4 @@ function indfact(OpRot, Mod)
     OpRot.ElPerf.indVelTg = b .* OpRot.ElPerf.tgSpeed;
     OpRot.ElPerf.inflAngle = phi;
     OpRot.ElPerf.alpha = alpha;
-    OpRot.ElPerf.reynolds = reynolds;
-
 end
diff --git a/src/solvers/indvel.m b/src/solvers/indvel.m
index 2eff70d92bfdb1491b8bcedc6e8befce7ae3d932..580da9a8d6ca8a02a4a90f8ebbbf4a296f1765bb 100644
--- a/src/solvers/indvel.m
+++ b/src/solvers/indvel.m
@@ -119,8 +119,6 @@ function indvel(OpRot, Mod)
     OpRot.ElPerf.alpha = alpha;
     OpRot.ElPerf.indVelAx = v - OpRot.ElPerf.upstreamVelAx;
     OpRot.ElPerf.indVelTg = uw / 2;
-    OpRot.ElPerf.reynolds = reynolds;
-
 end
 
 function force = coeff2force(coeff, OpRot, speed)
diff --git a/src/solvers/leishman.m b/src/solvers/leishman.m
index 0160198410f7b1c5b69b92344c107fdb788ec294..54b69783b31fb333092d8d9bda57a48ce366613b 100644
--- a/src/solvers/leishman.m
+++ b/src/solvers/leishman.m
@@ -100,8 +100,6 @@ function leishman(OpRot, Mod)
     OpRot.ElPerf.inflAngle = phi;
     OpRot.ElPerf.indVelAx = lambda * OpRot.tgTipSpeed - OpRot.Op.speed;
     OpRot.ElPerf.indVelTg = OpRot.ElPerf.tgSpeed - xi * OpRot.tgTipSpeed;
-    OpRot.ElPerf.reynolds = reynolds;
-
 end
 
 % --------------------------------------------------------------------------------------------------
diff --git a/src/solvers/stahlhut.m b/src/solvers/stahlhut.m
index e79e1cb2a335537016ddda1d5d08b43447821ef2..15b8dce0850024808cf6a2eca1201c854011d064 100644
--- a/src/solvers/stahlhut.m
+++ b/src/solvers/stahlhut.m
@@ -124,12 +124,6 @@ function stahlhut(OpRot, Mod)
     OpRot.ElPerf.swirlRat = OpRot.Rot.Bl.r .* cos(OpRot.ElPerf.inflAngle) ./ b2phi;
     OpRot.ElPerf.inflowRat = OpRot.ElPerf.swirlRat .* tan(OpRot.ElPerf.inflAngle);
 
-    indVelAx = OpRot.ElPerf.inflowRat .* OpRot.tgTipSpeed;
-    indVelTg = OpRot.ElPerf.tgSpeed - OpRot.ElPerf.swirlRat .* OpRot.tgTipSpeed;
-    relVel = sqrt((OpRot.Op.speed + indVelAx).^2 + (OpRot.ElPerf.tgSpeed - indVelTg).^2);
-    OpRot.ElPerf.reynolds = Flow.reynolds(relVel, ...
-                                          OpRot.Rot.Bl.chord, OpRot.Op.Flow.mu, OpRot.Op.Flow.rho);
-
 end
 
 function [gphi, b1phi, b2phi] = stahlhuteq(i, phi, OpRot, lossType)