fix initial wake velocities

MITRotor/BEMSolver.py

@@ -41,10 +41,8 @@ class BEMSolution:
     geom: BEMGeometry = field(repr=False)
     converged: bool
     niter: int
-
-    def __post_init__(self):
-        sol = Heck()(self.Ctprime(), self.yaw)
-        self._u4, self._v4 = sol.u4, sol.v4
+    u4: float
+    v4: float
 
     def a(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.an, grid)
@@ -99,12 +97,6 @@ def Ctau_uncorr(self, grid: Literal["sector ", "annulus", "rotor"] = "rotor"):
 
     def F(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         return average(self.geom, self.aero_props.F, grid)
-
-    def u4(self):
-        return self._u4
-
-    def v4(self):
-        return self._v4
 
     def Cp(self, grid: Literal["sector", "annulus", "rotor"] = "rotor"):
         dCp = (
@@ -218,5 +210,7 @@ def post_process(self, result: FixedPointIterationResult, pitch, tsr, yaw, U=Non
         an, aprime = result.x
         aero_props = self.aerodynamic_model(an, aprime, pitch, tsr, yaw, self.rotor, self.geometry, U, wdir)
         aero_props.F = self.tiploss_model(aero_props, pitch, tsr, yaw, self.rotor, self.geometry)
+        avg_Ct = average(self.geometry, aero_props.C_x)
+        u4,v4 = self.momentum_model.compute_initial_wake_velocities(avg_Ct, yaw)
 
-        return BEMSolution(pitch, tsr, yaw, aero_props, self.geometry, result.converged, result.niter)
+        return BEMSolution(pitch, tsr, yaw, aero_props, self.geometry, result.converged, result.niter, u4, v4)

MITRotor/Momentum.py

@@ -21,10 +21,6 @@
 
 
 class MomentumModel(ABC):
-    @abstractmethod
-    def compute_induction(self, Cx: ArrayLike, yaw: float) -> ArrayLike:
-        ...
-
     @abstractmethod
     def __call__(
         self,
@@ -36,6 +32,16 @@ def __call__(
         geom: "BEMGeometry",
     ) -> ArrayLike:
         ...
+
+    @abstractmethod
+    def compute_induction(self, Cx: ArrayLike, yaw: float) -> ArrayLike:
+        ...
+
+    @abstractmethod
+    def compute_initial_wake_velocities(self, Ct: float, yaw: float) -> ArrayLike:
+        ...
+
+
 
     def _func_rotor(
         self,
@@ -109,9 +115,13 @@ class ConstantInduction(MomentumModel):
     def __init__(self, a = 1/3):
         self.a = a
 
-    def _func(self, aero_props, pitch, tsr, yaw, rotor, geom) -> ArrayLike:
+    def compute_induction(self, Cx, yaw) -> ArrayLike:
         return self.a * np.ones_like(yaw)
-
+
+    def compute_initial_wake_velocities(self, Ct: float, yaw: float) -> ArrayLike:
+        u4 = 1 - 2 * self.a
+        v4 = - (1/4) * Ct * np.sin(yaw)
+        return u4, v4
 
 
 class ClassicalMomentum(MomentumModel):
@@ -128,6 +138,11 @@ def __init__(self, averaging: Literal["sector", "annulus", "rotor"] = "rotor"):
 
     def compute_induction(self, Cx, yaw):
         return 0.5 * (1 - np.sqrt(1 - Cx))
+
+    def compute_initial_wake_velocities(self, Ct: float, yaw: float) -> ArrayLike:
+        u4 = np.sqrt(1 - Ct)
+        v4 = - (1/4) * Ct * np.sin(yaw)
+        return u4, v4
 
 
 
@@ -167,6 +182,12 @@ def compute_induction(self, Cx: ArrayLike, yaw: float) -> ArrayLike:
                 a[mask] = (Cx[mask] - Ctc) / slope + self.ac
 
         return a
+
+    def compute_initial_wake_velocities(self, Ct: float, yaw: float) -> ArrayLike:
+        a = self.compute_induction(Ct, yaw)
+        u4 = 1 - Ct /(2  * (1 - a))
+        v4 = - (1/4) * Ct * np.sin(yaw)
+        return u4, v4
 
 
 class UnifiedMomentum(MomentumModel):
@@ -192,6 +213,10 @@ def __init__(self, averaging: Literal["sector", "annulus", "rotor"] = "rotor", b
     def compute_induction(self, Cx: ArrayLike, yaw: float) -> ArrayLike:
         sol = self.model_Ct(Cx, yaw)
         return sol.an
+
+    def compute_initial_wake_velocities(self, Ct: float, yaw: float) -> ArrayLike:
+        sol = self.model_Ct(Ct, yaw)
+        return sol.u4[0], sol.v4[0]
 
 
 class MadsenMomentum(MomentumModel):
@@ -222,3 +247,8 @@ def compute_induction(self, Cx: ArrayLike, yaw: float) -> ArrayLike:
 
         an = Ct**3 * 0.0883 + Ct**2 * 0.0586 + Ct * 0.2460
         return an
+
+    def compute_initial_wake_velocities(self, Ct: float, yaw: float) -> ArrayLike:
+        u4 = np.sqrt(1 - Ct)
+        v4 = - (1/4) * Ct * np.sin(yaw)
+        return u4, v4

examples/example_01_basic_BEM.py

@@ -23,5 +23,5 @@
     print(f"Local thrust coefficient: {sol.Ctprime():2.2f}")
     print(f"Axial induction: {sol.a():2.2f}")
     print(f"Rotor-effective windspeed: {sol.U():2.2f}")
-    print(f"Far-wake streamwise velocity: {sol.u4():2.2f}")
-    print(f"Far-wake lateral velocity: {sol.v4():2.2f}")
+    print(f"Far-wake streamwise velocity: {sol.u4:2.2f}")
+    print(f"Far-wake lateral velocity: {sol.v4:2.2f}")