From 47f2b3bca4ce151dd9fc1d6baea5a8b3f4bbd077 Mon Sep 17 00:00:00 2001 From: Lynna Date: Thu, 17 Jul 2025 16:16:07 -0400 Subject: [PATCH 1/9] Added Jensen and TurbOPark wake models to imports upon initialization --- mitwindfarm/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mitwindfarm/__init__.py b/mitwindfarm/__init__.py index 0fe7af7..cbfee3e 100644 --- a/mitwindfarm/__init__.py +++ b/mitwindfarm/__init__.py @@ -2,6 +2,6 @@ from .Rotor import RotorSolution, AD, UnifiedAD, BEM, CosineRotor from .RotorGrid import Point, Line, Area from .Superposition import Linear, Niayifar, Quadratic, Dominant -from .Wake import WakeModel, GaussianWakeModel, GaussianWake, VariableKwGaussianWakeModel +from .Wake import WakeModel, GaussianWakeModel, GaussianWake, VariableKwGaussianWakeModel, JensenWake, JensenWakeModel, TurbOParkWake, TurbOParkWakeModel from .windfarm import WindfarmSolution, PartialWindfarmSolution, Windfarm, CosineWindfarm from .Windfield import Uniform, PowerLaw, Superimposed From 2c1657324dac929645f73bfb07a46694230bf1dc Mon Sep 17 00:00:00 2001 From: Lynna Date: Thu, 17 Jul 2025 16:19:33 -0400 Subject: [PATCH 2/9] Added Niayifar quadratic wake superposition option --- mitwindfarm/Superposition.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/mitwindfarm/Superposition.py b/mitwindfarm/Superposition.py index 5418854..7834b33 100644 --- a/mitwindfarm/Superposition.py +++ b/mitwindfarm/Superposition.py @@ -25,3 +25,7 @@ def __call__(self, base_windfield: Windfield, wakes: list[Wake]) -> Windfield: class Dominant(Superposition): def __call__(self, base_windfield: Windfield, wakes: list[Wake]) -> Windfield: return Superimposed(base_windfield, wakes, method="dominant") + +class NQuadratic(Superposition): + def __call__(self, base_windfield: Windfield, wakes: list[Wake]) -> Windfield: + return Superimposed(base_windfield, wakes, method="nquadratic") \ No newline at end of file From e97e2e59bbe6d4c381f8bc3a3f86ad8d7b0de375 Mon Sep 17 00:00:00 2001 From: Lynna Date: Thu, 17 Jul 2025 16:20:47 -0400 Subject: [PATCH 3/9] Added Jensen, top hat TurbOPark, Gaussian TurbOPark wake models --- mitwindfarm/Wake.py | 707 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 705 insertions(+), 2 deletions(-) diff --git a/mitwindfarm/Wake.py b/mitwindfarm/Wake.py index a82abe6..383e2fe 100644 --- a/mitwindfarm/Wake.py +++ b/mitwindfarm/Wake.py @@ -38,6 +38,709 @@ def __call__(self, x, y, z, rotor_sol: "RotorSolution") -> Wake: ... +class JensenWake(Wake): + """ + Attributes: + x: x-position of rotor in global coordinate frame + y: y-position of rotor in global coordinate frame + z: z-position of rotor in global coordinate frame + rotor_sol: Rotor solution + sigma: Proportionality constant for wake diameter used in Gaussian + kw: Constant coefficient for wake growth + xmax: Maximum x value evaluated + dx: Interval of x values evaluated + TIamb: Ambient turbulence intensity + + """ + def __init__( + self, + x: float, + y: float, + z: float, + rotor_sol: "RotorSolution", + sigma: float = None, + kw: float = 0.07, + TIamb: float = None, + xmax: float = 100.0, + dx: float = 0.05 + ): + self.x, self.y, self.z = x, y, z + self.rotor_sol = rotor_sol + self.sigma, self.kw = sigma, kw + self.TIamb = TIamb or 0.0 + self.xmax = xmax + self.dx = dx + + # Precompute centerline far downstream + self.x_centerline, self.y_centerline = self._centerline(xmax, dx) + + def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0) -> ArrayLike: + + # Into rotor coordinate frame + x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z + yc = self.centerline(x_glob) - self.y + r = np.sqrt((y - yc) ** 2 + z ** 2) + + # Calculate wake diameter for each streamwise coordinate + d = self._wake_diameter(x) + + # Calculate du for each streamwise coordinate + du = self._du(x, wake_diameter = d) + + # Set deficit to be only inside cone of rotor + deficit = np.zeros(np.shape(y)) + + for ii, _ in np.ndenumerate(deficit): + if (np.sqrt(r[ii] ** 2) <= 1/2*d[ii]) and (x[ii] >= 0): + deficit[ii] = du[ii] + else: + # Outside of wake, there is no deficit (velocity is the same as u_inf) + deficit[ii] = 0 + + return deficit + + + def _wake_diameter(self, x: ArrayLike) -> ArrayLike: + """ + Solves the normalized far-wake diameter + """ + + return 1 + 2 * self.kw * x + + def _du(self, x: ArrayLike, wake_diameter: Optional[float] = None) -> ArrayLike: + """ + Solves for nondimensionalized wake deficit velocity + """ + + d = self._wake_diameter(x) if wake_diameter is None else wake_diameter + + du = (1 - self.rotor_sol.u4 / self.rotor_sol.REWS) / d**2 + + return du + + def wake_added_turbulence(self, x_glob, y_glob, z_glob): + + # Placeholder of zeroes for now + return np.zeros(np.shape(x_glob)) + + def centerline(self, x_glob: ArrayLike) -> ArrayLike: + """ + Interpolates Eq. 6 from Shapiro, Gayme, and Meneveau, 2018 (same as for + GaussianWake) for centerline y position in global coordinates + """ + x = x_glob - self.x + + yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) + + return yc_temp * self.rotor_sol.v4 / self.rotor_sol.REWS + self.y + + def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: + """ + Based on principle from Shapiro, Gayme, and Meneveau, 2018, the + transverse velocity wake recovery should mirror the axial velocity + wake recovery. The centerline y position in global coordinates is + then computed numerically using the transverse velocity deficit, + based on Equation 9 in Shapiro, Gayme, and Meneveau, 2018. + """ + + _x = np.arange(0, max(xmax, 2 * dx), dx) + d = self._wake_diameter(_x) + + # Handle d = 0 -> centerline is self.y when d = 0 + d_mask = d > 0 + dv = np.zeros(np.shape(d)) + + dv[d_mask] = (1 - self.rotor_sol.v4 / self.rotor_sol.REWS) / d[d_mask]**2 + _yc = cumulative_trapezoid(-dv, dx=dx, initial=0) + + _yc[d < 0] = self.y + + return _x, _yc + + def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0) -> ArrayLike: + + # Calculate deficit + deficit = self.deficit(x_glob, y_glob, z_glob) + + return self.rotor_sol.REWS * deficit + + +class JensenWakeModel(WakeModel): + def __init__( + self, + sigma: float = None, + kw = 0.07, + xmax: float = 100.0, + dx: float = 0.05 + ): + self.sigma = sigma + self.kw = kw + self.xmax = xmax + self.dx = dx + + def __call__( + self, + x: ArrayLike, + y: ArrayLike, + z: ArrayLike, + rotor_sol: "RotorSolution", + TIamb: float = None + ) -> JensenWake: + return JensenWake( + x, + y, + z, + rotor_sol, + sigma = self.sigma, + kw = self.kw, + TIamb = TIamb, + xmax = self.xmax, + dx = self.dx + ) + + +class TurbOParkWake(Wake): + def __init__( + self, + x: float, + y: float, + z: float, + rotor_sol: "RotorSolution", + TIamb: float = 0.1, + WATI_Iw_multiplier: float = 0.04, + c_1: float = 1.5, + c_2: float = 0.8, + xmax: float = 100.0, + dx: float = 0.05 + ): + self.x, self.y, self.z = x, y, z + self.rotor_sol = rotor_sol + self.TIamb = TIamb + self.WATI_Iw_multiplier = WATI_Iw_multiplier + self.c_1 = c_1 + self.c_2 = c_2 + + # Precompute centerline far downstream + self.x_centerline, self.y_centerline = self._centerline(xmax, dx) + + def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): + """ + Solves Eq. 6 in Pedersen et al. 2022 for wake deficit profile + """ + + x = x_glob - self.x + y = y_glob - self.y + z = z_glob - self.z + yc = self.centerline(x_glob) - self.y + r = np.sqrt((y - yc) ** 2 + z ** 2) + + sigma = self._char_wake_diameter(x) + + # Only compute behind rotor (x > 0) + x_mask = x > 0 + + # Handle sigma = 0, as this gives 0 in denominator of gaussian's exponent + # -> gaussian, du, C = 0 when sigma = 0 + sigma_mask = sigma > 0 + + comb_mask = x_mask & sigma_mask + + gaussian = np.zeros(np.shape(sigma)) + du = np.zeros(np.shape(sigma)) + + du[comb_mask] = self._du(x[comb_mask], sigma[comb_mask]) + + # Gaussian deficit formulation proposed by Bastankhah and Porte-Agel + gaussian[comb_mask] = np.exp(- r[comb_mask] ** 2 / + (2 * sigma[comb_mask] ** 2)) + + return gaussian * du + + def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): + """ + Computes wake deficit profile, weighted by rotor effective wind speed, + for use in Niayifar superposition method + """ + + x = x_glob - self.x + y = y_glob - self.y + z = z_glob - self.z + yc = self.centerline(x_glob) - self.y + r = np.sqrt((y - yc) ** 2 + z ** 2) + + sigma = self._char_wake_diameter(x) + + # Only compute behind rotor (x > 0) + x_mask = x > 0 + + # Handle sigma = 0, as this gives 0 in denominator of gaussian's exponent + # -> gaussian, du, C = 0 when sigma = 0 + sigma_mask = sigma > 0 + comb_mask = x_mask & sigma_mask + + gaussian = np.zeros(np.shape(sigma)) + du = np.zeros(np.shape(sigma)) + + du[comb_mask] = self._du(x[comb_mask], sigma[comb_mask]) + + # Gaussian deficit formulation proposed by Bastankhah and Porte-Agel + gaussian[comb_mask] = np.exp(- r[comb_mask] ** 2 / + (2 * sigma[comb_mask] ** 2)) + + return gaussian * du * self.rotor_sol.REWS + + def _centerline_wake_added_turb(self, x: ArrayLike): + """ + Solves Eq. 3 in Pedersen et al. 2022 for I_w(x), turbulence due to + wake at the centerline + """ + + I_w = 1 / (self.c_1 + self.c_2 * x/(np.sqrt(self.rotor_sol.Ct / self.rotor_sol.REWS ** 2))) + I_w[x < 0] = 0.0 + + return I_w + + def _tot_turb_intensity(self, x: ArrayLike): + """ + Solves Eq. 2 in Pedersen et al. 2022 for I(x), total turbulence + intensity + """ + + I_0 = self.TIamb + I_w = self._centerline_wake_added_turb(x) + I_total = np.sqrt(I_0 ** 2 + I_w ** 2) + + return I_total + + def _char_wake_diameter(self, x: ArrayLike) -> ArrayLike: + """ + Solves Eq. 4 in Pedersen et al. 2022 for nondimensionalized characteristic + wake diameter sigma. Expression was derived through analytical integration. + """ + + epsilon = self._wake_diameter_prop_const() + + # alpha and beta are parameters introduced by Pedersen et al. for + # simpler expression evaluation + alpha = self.c_1 * self.TIamb + beta = self.c_2 * self.TIamb / np.sqrt(self.rotor_sol.Ct / self.rotor_sol.REWS ** 2) + + # Only compute behind rotor (x > 0) + x_mask = x > 0 + + # Handle log(x) undefined when x <= 0 + # -> sigma = 0 when log expression undefined + log_expr = ((np.sqrt((alpha + beta * x) ** 2 + 1) + 1) * alpha) / (( + np.sqrt(1 + alpha ** 2) + 1) * (alpha + beta * x)) + log_mask = log_expr > 0 + + comb_mask = x_mask & log_mask + + sigma = np.zeros(np.shape(x)) + + sigma[comb_mask] = epsilon + self.WATI_Iw_multiplier * self.TIamb / beta * ( + np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1) + - np.sqrt(1 + alpha ** 2) + - np.log(log_expr[comb_mask]) + ) + + return sigma + + def _wake_diameter_prop_const(self) -> np.float64: + """ + Solves Eq. 5 in Pedersen et al. 2022 for epsilon, wake diameter + at rotor + """ + + Ct = self.rotor_sol.Ct / self.rotor_sol.REWS ** 2 + + epsilon = 0.25 * np.sqrt( + (1 + np.sqrt(1 - Ct)) + / (2 * np.sqrt(1 - Ct))) + + return epsilon + + def _wake_diameter(self, x) -> ArrayLike: + """ + Compute wake diameter (d) from characteristic wake diameter (sigma) + using the relation sigma(x) = epsilon * d(x) + """ + + sigma = self._char_wake_diameter(x) + epsilon = self._wake_diameter_prop_const() + + d = sigma / epsilon + + return d + + def _du(self, x, sigma: Optional[float] = None) -> ArrayLike: + """ + Solves Eq. 7 in Pedersen et al. for C(x), axial deficit at centerline + """ + + sigma = self._char_wake_diameter(x) if sigma is None else sigma + + # Only compute behind rotor (x > 0) + x_mask = x > 0 + + # Handle sigma = 0 nonphysical (there can't be a deficit without a char_wake_diameter) + # -> du = 0 when sigma <= 0 + sigma_mask = sigma > 0 + + # Handle sqrt does not exist -> du does not exist when sqrt does not exist + sqrt_expr = np.zeros(np.shape(sigma)) + sqrt_expr[sigma_mask] = 1 - (self.rotor_sol.Ct / self.rotor_sol.REWS ** 2) / (8 * sigma[sigma_mask] ** 2) + sqrt_mask = sqrt_expr >= 0 + + comb_mask = x_mask & sigma_mask & sqrt_mask + + du = np.zeros(np.shape(sigma)) + + du[comb_mask] = 1 - np.sqrt(sqrt_expr[comb_mask]) + + return du + + def wake_added_turbulence( + self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob=0 + ) -> ArrayLike: + """ + Returns wake added turbulence intensity caused by a wake at particular + points in space. Laterally smeared with the gaussian twice as wide as + the wake deficit model. As recommended by Niayifar and Porte-Agel 2016 + """ + x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z + yc = self.centerline(x_glob) - self.y + + # Only compute behind rotor (x > 0) + x_mask = x > 0 + + WATI = np.zeros(np.shape(x)) + sigma = np.zeros(np.shape(x)) + d = np.zeros(np.shape(x)) + _gaussian = np.zeros(np.shape(x)) + + WATI[x_mask] = self._centerline_wake_added_turb(x[x_mask]) + sigma[x_mask] = self._char_wake_diameter(x[x_mask]) + d[x_mask] = self._wake_diameter(x[x_mask]) + + # Handle sigma must be nonzero as this gives 0 in denominator of gaussian's exponent + # -> gaussian = 0 when sigma = 0 + sigma_mask = sigma != 0 + sigma = np.where(sigma==0, np.nan, sigma) + + comb_mask = x_mask & sigma_mask + + _gaussian[comb_mask] = ( + 1 + / (8 * (self.WATI_Iw_multiplier * sigma[comb_mask]) ** 2) + * np.exp( + -( + ((y[comb_mask] - yc[comb_mask]) ** 2 + z[comb_mask]**2) + / (2 * (self.WATI_Iw_multiplier * sigma[comb_mask]) ** 2 * d[comb_mask]**2) + ) + ) + ) + + return _gaussian * np.nan_to_num(WATI) + + # [***FLAG***] Is this same model by Niayifar and Porte Agel applicable for TurbOPark? + + def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: + """ + Based on principle from Shapiro, Gayme, and Meneveau, 2018, the + transverse velocity wake recovery should mirror the axial velocity + wake recovery. The centerline y position in global coordinates is + then computed numerically using the transverse velocity deficit, + based on Equation 9 in Shapiro, Gayme, and Meneveau, 2018. + """ + + _x = np.arange(0, max(xmax, 2 * dx), dx) + sigma = self._char_wake_diameter(_x) + + # Only compute behind rotor (x > 0) + x_mask = _x > 0 + + # Handle sigma = 0 nonphysical (there can't be a deficit without a char_wake_diameter) + # -> du = 0 when sigma <= 0 + sigma_mask = sigma > 0 + + # Handle sqrt does not exist -> du does not exist when sqrt does not exist + sqrt_expr = np.zeros(np.shape(sigma)) + sqrt_expr[sigma_mask] = 1 - (self.rotor_sol.Ct / self.rotor_sol.REWS ** 2) / (8 * sigma[sigma_mask] ** 2) + sqrt_mask = sqrt_expr >= 0 + + comb_mask = x_mask & sigma_mask & sqrt_mask + + dv = np.zeros(np.shape(sigma)) + + dv[comb_mask] = 1 - np.sqrt(sqrt_expr[comb_mask]) + _yc = cumulative_trapezoid(-dv, dx=dx, initial=0) + + _yc[sigma < 0] = self.y + + return _x, _yc + + # def _centerline_dv(self, ) + + def centerline(self, x_glob: ArrayLike) -> ArrayLike: + """ + Interpolates Eq. 6 from Shapiro, Gayme, and Meneveau, 2018 (same as for + GaussianWake) for centerline y position in global coordinates + """ + x = x_glob - self.x + + yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) + + return yc_temp * self.rotor_sol.extra.v4 + self.y + + +class TopHatTurbOParkWake(Wake): + def __init__( + self, + x: float, + y: float, + z: float, + rotor_sol: "RotorSolution", + TIamb: float = 0.1, + WATI_Iw_multiplier: float = 0.04, + c_1: float = 1.5, + c_2: float = 0.8, + xmax: float = 100.0, + dx: float = 0.05 + ): + self.x, self.y, self.z = x, y, z + self.rotor_sol = rotor_sol + self.TIamb = TIamb + self.WATI_Iw_multiplier = WATI_Iw_multiplier + self.c_1 = c_1 + self.c_2 = c_2 + + # Precompute centerline far downstream + self.x_centerline, self.y_centerline = self._centerline(xmax, dx) + + def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): + """ + Solves Eq. 6 in Pedersen et al. 2022 for wake deficit profile + """ + + x = x_glob - self.x + y = y_glob - self.y + z = z_glob - self.z + yc = self.centerline(x_glob) - self.y + r = np.sqrt((y - yc) ** 2 + z ** 2) + + d = self._wake_diameter(x) + + # Only compute for x > 0 + x_mask = x > 0 + + du = np.zeros(np.shape(d)) + + du[x_mask] = self._du(x[x_mask], d[x_mask]) + + # Set deficit to be only inside cone of rotor + deficit = np.zeros(np.shape(y)) + + for ii, _ in np.ndenumerate(deficit): + if (np.sqrt(r[ii] ** 2) <= 1/2*d[ii]) and (x[ii] >= 0): + deficit[ii] = du[ii] + else: + # Outside of wake, there is no deficit (velocity is the same as u_inf) + deficit[ii] = 0 + + return deficit + + def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): + """ + Computes wake deficit profile, weighted by rotor effective wind speed, + for use in Niayifar superposition method + """ + + # Calculate deficit + deficit = self.deficit(x_glob, y_glob, z_glob) + + return deficit * self.rotor_sol.REWS + + def _centerline_wake_added_turb(self, x: ArrayLike): + """ + Solves Eq. 3 in Pedersen et al. 2022 for I_w(x), turbulence due to + wake at the centerline + """ + + I_w = 1 / (self.c_1 + self.c_2 * x/(np.sqrt(self.rotor_sol.Ct / self.rotor_sol.REWS ** 2))) + I_w[x < 0] = 0.0 + + return I_w + + def _tot_turb_intensity(self, x: ArrayLike): + """ + Solves Eq. 2 in Pedersen et al. 2022 for I(x), total turbulence + intensity + """ + + I_0 = self.TIamb + I_w = self._centerline_wake_added_turb(x) + I_total = np.sqrt(I_0 ** 2 + I_w ** 2) + + return I_total + + def _wake_diameter(self, x) -> ArrayLike: + """ + Compute wake diameter (d) from characteristic wake diameter (sigma) + using the relation sigma(x) = epsilon * d(x) + """ + + # alpha and beta are parameters introduced by Nygaard et al. for + # simpler expression evaluation + alpha = self.c_1 * self.TIamb + # beta = self.c_2 * self.TIamb / np.sqrt(self.rotor_sol.Ct) + beta = self.c_2 * self.TIamb / np.sqrt(self.rotor_sol.Ct / self.rotor_sol.REWS ** 2) + + # Only compute behind rotor (x > 0) + x_mask = x > 0 + + # Handle log(x) undefined when x <= 0 + # -> d = 0 when log expression undefined + log_expr = ((np.sqrt((alpha + beta * x) ** 2 + 1) + 1) * alpha) / (( + np.sqrt(1 + alpha ** 2) + 1) * (alpha + beta * x)) + log_mask = log_expr > 0 + + comb_mask = x_mask & log_mask + + d = np.zeros(np.shape(x)) + + d[comb_mask] = 1 + self.WATI_Iw_multiplier * self.TIamb / beta * ( + np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1) + - np.sqrt(1 + alpha ** 2) + - np.log(log_expr[comb_mask]) + ) + + return d + + def _du(self, x, wake_diameter: Optional[float] = None) -> ArrayLike: + """ + Solves Eq. 7 in Pedersen et al. for C(x), deficit at centerline + """ + + d = self._wake_diameter(x) if wake_diameter is None else wake_diameter + + # Only compute for x > 0 + x_mask = x > 0 + + # Handle d = 0 -> du = 0 when sigma = 0 + d_mask = d > 0 + + comb_mask = x_mask & d_mask + + du = np.zeros(np.shape(d)) + + # du[comb_mask] = (1 - self.rotor_sol.REWS * np.sqrt(1 - self.rotor_sol.Ct)) / (d[comb_mask]**2) + du[comb_mask] = (1 - self.rotor_sol.REWS * np.sqrt(1 - self.rotor_sol.Ct / self.rotor_sol.REWS ** 2)) / (d[comb_mask]**2) + + return du + + def wake_added_turbulence( + self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob=0 + ) -> ArrayLike: + + # Placeholder of zeroes for now + return np.zeros(np.shape(x_glob)) + + + def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: + """ + Based on principle from Shapiro, Gayme, and Meneveau, 2018, the + transverse velocity wake recovery should mirror the axial velocity + wake recovery. The centerline y position in global coordinates is + then computed numerically using the transverse velocity deficit, + based on Equation 9 in Shapiro, Gayme, and Meneveau, 2018. + """ + + _x = np.arange(0, max(xmax, 2 * dx), dx) + d = self._wake_diameter(_x) + + # Handle d = 0 -> centerline is self.y when d = 0 + d_mask = d > 0 + dv = np.zeros(np.shape(d)) + + # [***FLAG***] I'm not sure if this would be the right way to handle the proportionality + # between deficit in v and deficit in u + dv[d_mask] = (1 - np.tan(self.rotor_sol.yaw) * np.sqrt(1 - self.rotor_sol.Ct / self.rotor_sol.REWS ** 2)) / (d[d_mask]**2) + _yc = cumulative_trapezoid(-dv, dx=dx, initial=0) + + _yc[d < 0] = self.y + + return _x, _yc + + def centerline(self, x_glob: ArrayLike) -> ArrayLike: + """ + Interpolates Eq. 6 from Shapiro, Gayme, and Meneveau, 2018 (same as for + GaussianWake) for centerline y position in global coordinates + """ + x = x_glob - self.x + + yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) + + return yc_temp * self.rotor_sol.extra.v4 + self.y + + +class TurbOParkWakeModel(WakeModel): + def __init__( + self, + WATI_Iw_multiplier: float = 0.04, + c_1: float = 1.5, + c_2: float = 0.8, + xmax: float = 100.0, + dx: float = 0.05, + gaussian_profile: bool = True + ): + self.xmax = xmax + self.dx = dx + self.WATI_Iw_multiplier = WATI_Iw_multiplier + self.c_1 = c_1 + self.c_2 = c_2 + self.gaussian_profile = gaussian_profile + + def __call__( + self, + x, + y, + z, + rotor_sol: "RotorSolution", + TIamb: float = 0.06 + ) -> TurbOParkWake: + + # Return Gaussian TurbOPark by default (Pedersen et al. 2022) + if self.gaussian_profile == True: + return TurbOParkWake( + x, + y, + z, + rotor_sol, + TIamb = TIamb, + xmax = self.xmax, + dx = self.dx, + WATI_Iw_multiplier = self.WATI_Iw_multiplier, + c_1 = self.c_1, + c_2 = self.c_2, + ) + + # Return Top Hat TurbOPark if desired (Nygaard et al. 2020) + else: + return TopHatTurbOParkWake( + x, + y, + z, + rotor_sol, + TIamb = TIamb, + xmax = self.xmax, + dx = self.dx, + WATI_Iw_multiplier = self.WATI_Iw_multiplier, + c_1 = self.c_1, + c_2 = self.c_2, + ) + + class GaussianWake(Wake): def __init__( self, @@ -85,7 +788,7 @@ def centerline(self, x_glob: ArrayLike) -> ArrayLike: yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) - return yc_temp * self.rotor_sol.v4 + self.y + return yc_temp * self.rotor_sol.extra.v4 + self.y def centerline_wake_added_turb(self, x: ArrayLike) -> ArrayLike: """ @@ -161,7 +864,7 @@ def wake_added_turbulence( """ Returns wake added turbulence intensity caused by a wake at particular points in space. Laterally smeared with the gaussian twice as wide as - the wake deficit model. as recommended by Niayifar and Porte-Agel 2016 + the wake deficit model. As recommended by Niayifar and Porte-Agel 2016 """ x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z d = self._wake_diameter(x) From dd617df925004826fd0a127bf325aa0b194b7509 Mon Sep 17 00:00:00 2001 From: Lynna Date: Thu, 17 Jul 2025 16:30:26 -0400 Subject: [PATCH 4/9] Added memory-friendly superimposed wind speed computation --- mitwindfarm/Windfield.py | 59 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 57 insertions(+), 2 deletions(-) diff --git a/mitwindfarm/Windfield.py b/mitwindfarm/Windfield.py index a139f56..c47708f 100644 --- a/mitwindfarm/Windfield.py +++ b/mitwindfarm/Windfield.py @@ -213,7 +213,7 @@ def wsp(self, x: ArrayLike, y: ArrayLike, z: ArrayLike) -> ArrayLike: wsp_base = self.base_windfield.wsp(x, y, z) deficits = [] for wake in self.wakes: - if self.method == "niayifar": + if (self.method == "niayifar") | (self.method == "nquadratic"): deficits.append(wake.niayifar_deficit(x, y, z)) else: deficits.append(wake.deficit(x, y, z)) @@ -223,7 +223,7 @@ def wsp(self, x: ArrayLike, y: ArrayLike, z: ArrayLike) -> ArrayLike: if (self.method == "linear") | (self.method == "niayifar"): wsp_out = wsp_base - np.sum(deficits, axis=0) - elif self.method == "quadratic": + elif (self.method == "quadratic") | (self.method == "nquadratic"): wsp_out = wsp_base - np.sqrt(np.sum(np.array(deficits)**2, axis=0)) elif self.method == "dominant": wsp_out = wsp_base - np.array(deficits).max(axis=0, initial=0) @@ -231,6 +231,61 @@ def wsp(self, x: ArrayLike, y: ArrayLike, z: ArrayLike) -> ArrayLike: raise NotImplementedError return wsp_out + + + def mem_eff_wsp(self, x: ArrayLike, y: ArrayLike, z: ArrayLike) -> ArrayLike: + """ + Returns an array of wind speed based on superposing the windfield's wakes. Runs without creating a + x by y by z by n_rotors array to conserve memory. + + Parameters: + - x: x-coordinates. + - y: y-coordinates. + - z: z-coordinates. + + Returns: + ArrayLike: Array of ind speed based on the windfield's wakes with the same shape as input coordinates. + """ + wsp_base = self.base_windfield.wsp(x, y, z) + tot_deficit = np.zeros_like(wsp_base) + deficit_count = 0 + + for wake in self.wakes: + if (self.method == "niayifar") | (self.method == "nquadratic"): + deficit = wake.niayifar_deficit(x, y, z) + else: + deficit = wake.deficit(x, y, z) + + if deficit is not None: + deficit_count += 1 + + # Compute linear methods by first totalling deficit fields + if (self.method == "linear") | (self.method == "niayifar"): + tot_deficit += deficit + # Compute quadratic methods by first summing squares of deficit fields + elif (self.method == "quadratic") | (self.method == "nquadratic"): + tot_deficit += deficit**2 + # Compute dominant method by first checking maximum of deficit fields + elif self.method == "dominant": + # Stack current deficit maximum and deficit + deficits = np.stack((tot_deficit, deficit), axis = 0) + # Store maximum of the two + tot_deficit = deficits.max(axis = 0, initial = 0) + else: + raise NotImplementedError + + # Handle no deficits + if deficit_count == 0: + wsp_out = wsp_base + # Output linear and dominant methods by subtracting total of deficits from base + elif (self.method == "linear") | (self.method == "niayifar") | (self.method == "dominant"): + wsp_out = wsp_base - tot_deficit + # Output quadratic methods by subtracting square root of tot_deficit from base + elif (self.method == "quadratic") | (self.method == "nquadratic"): + wsp_out = wsp_base - np.sqrt(tot_deficit) + + return wsp_out + def TI(self, x: ArrayLike, y: ArrayLike, z: ArrayLike) -> ArrayLike: """ From b5bc028b9a5d2705931b87c387d055596934f9e6 Mon Sep 17 00:00:00 2001 From: dengl Date: Fri, 8 Aug 2025 12:18:48 -0400 Subject: [PATCH 5/9] Renamed A parameter and added comments --- mitwindfarm/Wake.py | 70 ++++++++++++++++++++++++++++++--------------- 1 file changed, 47 insertions(+), 23 deletions(-) diff --git a/mitwindfarm/Wake.py b/mitwindfarm/Wake.py index 383e2fe..23f859f 100644 --- a/mitwindfarm/Wake.py +++ b/mitwindfarm/Wake.py @@ -45,12 +45,12 @@ class JensenWake(Wake): y: y-position of rotor in global coordinate frame z: z-position of rotor in global coordinate frame rotor_sol: Rotor solution - sigma: Proportionality constant for wake diameter used in Gaussian + sigma: Proportionality constant for wake diameter used in Gaussian and thus expected, + even though Jensen doesn't require a sigma kw: Constant coefficient for wake growth xmax: Maximum x value evaluated dx: Interval of x values evaluated TIamb: Ambient turbulence intensity - """ def __init__( self, @@ -120,7 +120,7 @@ def _du(self, x: ArrayLike, wake_diameter: Optional[float] = None) -> ArrayLike: def wake_added_turbulence(self, x_glob, y_glob, z_glob): - # Placeholder of zeroes for now + # Placeholder of zeroes return np.zeros(np.shape(x_glob)) def centerline(self, x_glob: ArrayLike) -> ArrayLike: @@ -162,6 +162,7 @@ def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0) -> # Calculate deficit deficit = self.deficit(x_glob, y_glob, z_glob) + # Multiply by REWS to get Niayifar deficit return self.rotor_sol.REWS * deficit @@ -200,6 +201,18 @@ def __call__( class TurbOParkWake(Wake): + """ + Attributes: + x: x-position of rotor in global coordinate frame + y: y-position of rotor in global coordinate frame + z: z-position of rotor in global coordinate frame + rotor_sol: Rotor solution + TIamb: Ambient turbulence intensity + A: Wake expansion calibration parameter + c_1, c_2: Empirically estimated constants in wake-added turbulence intensity + xmax: Maximum x value evaluated + dx: Interval of x values evaluated + """ def __init__( self, x: float, @@ -207,7 +220,7 @@ def __init__( z: float, rotor_sol: "RotorSolution", TIamb: float = 0.1, - WATI_Iw_multiplier: float = 0.04, + A: float = 0.04, c_1: float = 1.5, c_2: float = 0.8, xmax: float = 100.0, @@ -216,7 +229,7 @@ def __init__( self.x, self.y, self.z = x, y, z self.rotor_sol = rotor_sol self.TIamb = TIamb - self.WATI_Iw_multiplier = WATI_Iw_multiplier + self.A = A self.c_1 = c_1 self.c_2 = c_2 @@ -338,7 +351,7 @@ def _char_wake_diameter(self, x: ArrayLike) -> ArrayLike: sigma = np.zeros(np.shape(x)) - sigma[comb_mask] = epsilon + self.WATI_Iw_multiplier * self.TIamb / beta * ( + sigma[comb_mask] = epsilon + self.A * self.TIamb / beta * ( np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1) - np.sqrt(1 + alpha ** 2) - np.log(log_expr[comb_mask]) @@ -405,7 +418,7 @@ def wake_added_turbulence( ) -> ArrayLike: """ Returns wake added turbulence intensity caused by a wake at particular - points in space. Laterally smeared with the gaussian twice as wide as + points in space. Laterally smeared with the Gaussian twice as wide as the wake deficit model. As recommended by Niayifar and Porte-Agel 2016 """ x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z @@ -432,18 +445,19 @@ def wake_added_turbulence( _gaussian[comb_mask] = ( 1 - / (8 * (self.WATI_Iw_multiplier * sigma[comb_mask]) ** 2) + / (8 * (self.A * sigma[comb_mask]) ** 2) * np.exp( -( ((y[comb_mask] - yc[comb_mask]) ** 2 + z[comb_mask]**2) - / (2 * (self.WATI_Iw_multiplier * sigma[comb_mask]) ** 2 * d[comb_mask]**2) + / (2 * (self.A * sigma[comb_mask]) ** 2 * d[comb_mask]**2) ) ) ) return _gaussian * np.nan_to_num(WATI) - # [***FLAG***] Is this same model by Niayifar and Porte Agel applicable for TurbOPark? + # TODO: Not exactly what is given in TurbOPark, + # Check if same model by Niayifar and Porte Agel applicable def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: """ @@ -480,8 +494,6 @@ def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: return _x, _yc - # def _centerline_dv(self, ) - def centerline(self, x_glob: ArrayLike) -> ArrayLike: """ Interpolates Eq. 6 from Shapiro, Gayme, and Meneveau, 2018 (same as for @@ -491,10 +503,22 @@ def centerline(self, x_glob: ArrayLike) -> ArrayLike: yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) - return yc_temp * self.rotor_sol.extra.v4 + self.y + return yc_temp * self.rotor_sol.v4 / self.rotor_sol.REWS + self.y class TopHatTurbOParkWake(Wake): + """ + Attributes: + x: x-position of rotor in global coordinate frame + y: y-position of rotor in global coordinate frame + z: z-position of rotor in global coordinate frame + rotor_sol: Rotor solution + TIamb: Ambient turbulence intensity + A: Wake expansion calibration parameter + c_1, c_2: Empirically estimated constants in wake-added turbulence intensity + xmax: Maximum x value evaluated + dx: Interval of x values evaluated + """ def __init__( self, x: float, @@ -502,7 +526,7 @@ def __init__( z: float, rotor_sol: "RotorSolution", TIamb: float = 0.1, - WATI_Iw_multiplier: float = 0.04, + A: float = 0.04, c_1: float = 1.5, c_2: float = 0.8, xmax: float = 100.0, @@ -511,7 +535,7 @@ def __init__( self.x, self.y, self.z = x, y, z self.rotor_sol = rotor_sol self.TIamb = TIamb - self.WATI_Iw_multiplier = WATI_Iw_multiplier + self.A = A self.c_1 = c_1 self.c_2 = c_2 @@ -609,7 +633,7 @@ def _wake_diameter(self, x) -> ArrayLike: d = np.zeros(np.shape(x)) - d[comb_mask] = 1 + self.WATI_Iw_multiplier * self.TIamb / beta * ( + d[comb_mask] = 1 + self.A * self.TIamb / beta * ( np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1) - np.sqrt(1 + alpha ** 2) - np.log(log_expr[comb_mask]) @@ -663,7 +687,7 @@ def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: d_mask = d > 0 dv = np.zeros(np.shape(d)) - # [***FLAG***] I'm not sure if this would be the right way to handle the proportionality + # TODO: Not sure if this would be the right way to handle the proportionality # between deficit in v and deficit in u dv[d_mask] = (1 - np.tan(self.rotor_sol.yaw) * np.sqrt(1 - self.rotor_sol.Ct / self.rotor_sol.REWS ** 2)) / (d[d_mask]**2) _yc = cumulative_trapezoid(-dv, dx=dx, initial=0) @@ -681,13 +705,13 @@ def centerline(self, x_glob: ArrayLike) -> ArrayLike: yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) - return yc_temp * self.rotor_sol.extra.v4 + self.y + return yc_temp * self.rotor_sol.v4 / self.rotor_sol.REWS + self.y class TurbOParkWakeModel(WakeModel): def __init__( self, - WATI_Iw_multiplier: float = 0.04, + A: float = 0.04, c_1: float = 1.5, c_2: float = 0.8, xmax: float = 100.0, @@ -696,7 +720,7 @@ def __init__( ): self.xmax = xmax self.dx = dx - self.WATI_Iw_multiplier = WATI_Iw_multiplier + self.A = A self.c_1 = c_1 self.c_2 = c_2 self.gaussian_profile = gaussian_profile @@ -720,7 +744,7 @@ def __call__( TIamb = TIamb, xmax = self.xmax, dx = self.dx, - WATI_Iw_multiplier = self.WATI_Iw_multiplier, + A = self.A, c_1 = self.c_1, c_2 = self.c_2, ) @@ -735,7 +759,7 @@ def __call__( TIamb = TIamb, xmax = self.xmax, dx = self.dx, - WATI_Iw_multiplier = self.WATI_Iw_multiplier, + A = self.A, c_1 = self.c_1, c_2 = self.c_2, ) @@ -788,7 +812,7 @@ def centerline(self, x_glob: ArrayLike) -> ArrayLike: yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) - return yc_temp * self.rotor_sol.extra.v4 + self.y + return yc_temp * self.rotor_sol.v4 / self.rotor_sol.REWS + self.y def centerline_wake_added_turb(self, x: ArrayLike) -> ArrayLike: """ From 42e95851d7a02a6b51f65e310188fcf42225f83f Mon Sep 17 00:00:00 2001 From: dengl Date: Sun, 14 Sep 2025 11:54:28 -0400 Subject: [PATCH 6/9] Added tilt compatibility and renamed A parameter --- mitwindfarm/Wake.py | 474 +++++++++++++++++++++----------------------- 1 file changed, 225 insertions(+), 249 deletions(-) diff --git a/mitwindfarm/Wake.py b/mitwindfarm/Wake.py index 23f859f..ba5e00d 100644 --- a/mitwindfarm/Wake.py +++ b/mitwindfarm/Wake.py @@ -1,5 +1,6 @@ from abc import ABC, abstractmethod from typing import Optional, TYPE_CHECKING +import warnings import numpy as np from numpy.typing import ArrayLike @@ -45,12 +46,12 @@ class JensenWake(Wake): y: y-position of rotor in global coordinate frame z: z-position of rotor in global coordinate frame rotor_sol: Rotor solution - sigma: Proportionality constant for wake diameter used in Gaussian and thus expected, - even though Jensen doesn't require a sigma + sigma: Proportionality constant for wake diameter used in Gaussian kw: Constant coefficient for wake growth xmax: Maximum x value evaluated dx: Interval of x values evaluated TIamb: Ambient turbulence intensity + """ def __init__( self, @@ -72,14 +73,17 @@ def __init__( self.dx = dx # Precompute centerline far downstream - self.x_centerline, self.y_centerline = self._centerline(xmax, dx) + self.x_centerline, self.y_centerline, self.z_centerline = self._centerline(xmax, dx) def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0) -> ArrayLike: # Into rotor coordinate frame - x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z - yc = self.centerline(x_glob) - self.y - r = np.sqrt((y - yc) ** 2 + z ** 2) + x = np.float64(x_glob - self.x) + y = np.float64(y_glob - self.y) + z = np.float64(z_glob - self.z) + yc = self.centerline(x_glob)[0] - self.y + zc = self.centerline(x_glob)[1] - self.z + r = np.sqrt((y - yc) ** 2 + (z - zc) ** 2) # Calculate wake diameter for each streamwise coordinate d = self._wake_diameter(x) @@ -87,15 +91,28 @@ def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0) -> ArrayLike # Calculate du for each streamwise coordinate du = self._du(x, wake_diameter = d) - # Set deficit to be only inside cone of rotor - deficit = np.zeros(np.shape(y)) + # Set deficit to be only inside wake region + # Handle x is a constant + if np.shape(x) == (): + deficit = np.zeros(np.shape(y)) + + for ii, _ in np.ndenumerate(deficit): + if (np.sqrt(r[ii] ** 2) <= 1/2*d) and (x >= 0): + deficit[ii] = du + else: + # Outside of wake, there is no deficit (velocity is the same as u_inf) + deficit[ii] = 0 + + # Handle x is an array + else: + deficit = np.zeros(np.shape(x)) - for ii, _ in np.ndenumerate(deficit): - if (np.sqrt(r[ii] ** 2) <= 1/2*d[ii]) and (x[ii] >= 0): - deficit[ii] = du[ii] - else: - # Outside of wake, there is no deficit (velocity is the same as u_inf) - deficit[ii] = 0 + for ii, _ in np.ndenumerate(deficit): + if (np.sqrt(r[ii] ** 2) <= 1/2*d[ii]) and (x[ii] >= 0): + deficit[ii] = du[ii] + else: + # Outside of wake, there is no deficit (velocity is the same as u_inf) + deficit[ii] = 0 return deficit @@ -130,39 +147,48 @@ def centerline(self, x_glob: ArrayLike) -> ArrayLike: """ x = x_glob - self.x + # Interpolate previously computed centerlines yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) + zc_temp = np.interp(x, self.x_centerline, self.z_centerline, left=0) - return yc_temp * self.rotor_sol.v4 / self.rotor_sol.REWS + self.y + # Translate centerlines into global coordinates + yc = yc_temp + self.y + zc = zc_temp + self.z + + return yc, zc def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: """ + Returns centerline y and z position in rotor coordinates. + Based on principle from Shapiro, Gayme, and Meneveau, 2018, the - transverse velocity wake recovery should mirror the axial velocity + transverse velocity wake recovery should parallel the axial velocity wake recovery. The centerline y position in global coordinates is then computed numerically using the transverse velocity deficit, - based on Equation 9 in Shapiro, Gayme, and Meneveau, 2018. + based on Equation 3.9 in Shapiro, Gayme, and Meneveau, 2018. """ _x = np.arange(0, max(xmax, 2 * dx), dx) d = self._wake_diameter(_x) - # Handle d = 0 -> centerline is self.y when d = 0 + # Compute centerline y position when d != 0 d_mask = d > 0 dv = np.zeros(np.shape(d)) - - dv[d_mask] = (1 - self.rotor_sol.v4 / self.rotor_sol.REWS) / d[d_mask]**2 + dv[d_mask] = (- self.rotor_sol.v4 / self.rotor_sol.REWS) / d[d_mask]**2 _yc = cumulative_trapezoid(-dv, dx=dx, initial=0) - _yc[d < 0] = self.y + # Compute centerline z position when d != 0 + dw = np.zeros(np.shape(d)) + dw[d_mask] = (- self.rotor_sol.w4 / self.rotor_sol.REWS) / d[d_mask]**2 + _zc = cumulative_trapezoid(-dw, dx=dx, initial=0) - return _x, _yc + return _x, _yc, _zc def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0) -> ArrayLike: # Calculate deficit deficit = self.deficit(x_glob, y_glob, z_glob) - # Multiply by REWS to get Niayifar deficit return self.rotor_sol.REWS * deficit @@ -201,18 +227,6 @@ def __call__( class TurbOParkWake(Wake): - """ - Attributes: - x: x-position of rotor in global coordinate frame - y: y-position of rotor in global coordinate frame - z: z-position of rotor in global coordinate frame - rotor_sol: Rotor solution - TIamb: Ambient turbulence intensity - A: Wake expansion calibration parameter - c_1, c_2: Empirically estimated constants in wake-added turbulence intensity - xmax: Maximum x value evaluated - dx: Interval of x values evaluated - """ def __init__( self, x: float, @@ -220,7 +234,7 @@ def __init__( z: float, rotor_sol: "RotorSolution", TIamb: float = 0.1, - A: float = 0.04, + WATI_Iw_multiplier: float = 0.04, c_1: float = 1.5, c_2: float = 0.8, xmax: float = 100.0, @@ -229,23 +243,21 @@ def __init__( self.x, self.y, self.z = x, y, z self.rotor_sol = rotor_sol self.TIamb = TIamb - self.A = A + self.WATI_Iw_multiplier = WATI_Iw_multiplier self.c_1 = c_1 self.c_2 = c_2 - # Precompute centerline far downstream - self.x_centerline, self.y_centerline = self._centerline(xmax, dx) - def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): """ Solves Eq. 6 in Pedersen et al. 2022 for wake deficit profile """ - x = x_glob - self.x - y = y_glob - self.y - z = z_glob - self.z - yc = self.centerline(x_glob) - self.y - r = np.sqrt((y - yc) ** 2 + z ** 2) + x = np.float64(x_glob - self.x) + y = np.float64(y_glob - self.y) + z = np.float64(z_glob - self.z) + yc = self.centerline(x_glob)[0] - self.y + zc = self.centerline(x_glob)[1] - self.z + r = np.sqrt((y - yc) ** 2 + (z - zc) ** 2) sigma = self._char_wake_diameter(x) @@ -258,14 +270,30 @@ def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): comb_mask = x_mask & sigma_mask - gaussian = np.zeros(np.shape(sigma)) - du = np.zeros(np.shape(sigma)) + # Gaussian deficit formulation proposed by Bastankhah and Porte-Agel + # Handle x is a constant + if np.shape(x) == () and comb_mask: + + # Initialize Gaussian profile and peak deficit + gaussian = np.zeros(np.shape(y)) + du = np.zeros(np.shape(x)) - du[comb_mask] = self._du(x[comb_mask], sigma[comb_mask]) + du = self._du(x, sigma[comb_mask]) - # Gaussian deficit formulation proposed by Bastankhah and Porte-Agel - gaussian[comb_mask] = np.exp(- r[comb_mask] ** 2 / - (2 * sigma[comb_mask] ** 2)) + gaussian[comb_mask] = np.exp(- r[comb_mask] ** 2 / + (2 * sigma[comb_mask] ** 2)) + + # Handle x is an array + else: + + # Initialize Gaussian profile and peak deficit + gaussian = np.zeros(np.shape(x)) + du = np.zeros(np.shape(x)) + + du[comb_mask] = self._du(x[comb_mask], sigma[comb_mask]) + + gaussian[comb_mask] = np.exp(- r[comb_mask] ** 2 / + (2 * sigma[comb_mask] ** 2)) return gaussian * du @@ -275,32 +303,10 @@ def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): for use in Niayifar superposition method """ - x = x_glob - self.x - y = y_glob - self.y - z = z_glob - self.z - yc = self.centerline(x_glob) - self.y - r = np.sqrt((y - yc) ** 2 + z ** 2) - - sigma = self._char_wake_diameter(x) - - # Only compute behind rotor (x > 0) - x_mask = x > 0 - - # Handle sigma = 0, as this gives 0 in denominator of gaussian's exponent - # -> gaussian, du, C = 0 when sigma = 0 - sigma_mask = sigma > 0 - comb_mask = x_mask & sigma_mask - - gaussian = np.zeros(np.shape(sigma)) - du = np.zeros(np.shape(sigma)) - - du[comb_mask] = self._du(x[comb_mask], sigma[comb_mask]) - - # Gaussian deficit formulation proposed by Bastankhah and Porte-Agel - gaussian[comb_mask] = np.exp(- r[comb_mask] ** 2 / - (2 * sigma[comb_mask] ** 2)) + # Calculate deficit + deficit = self.deficit(x_glob, y_glob, z_glob) - return gaussian * du * self.rotor_sol.REWS + return self.rotor_sol.REWS * deficit def _centerline_wake_added_turb(self, x: ArrayLike): """ @@ -351,11 +357,11 @@ def _char_wake_diameter(self, x: ArrayLike) -> ArrayLike: sigma = np.zeros(np.shape(x)) - sigma[comb_mask] = epsilon + self.A * self.TIamb / beta * ( - np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1) - - np.sqrt(1 + alpha ** 2) - - np.log(log_expr[comb_mask]) - ) + sigma[comb_mask] = epsilon + self.WATI_Iw_multiplier * self.TIamb / beta * ( + np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1) + - np.sqrt(1 + alpha ** 2) + - np.log(log_expr[comb_mask]) + ) return sigma @@ -418,11 +424,12 @@ def wake_added_turbulence( ) -> ArrayLike: """ Returns wake added turbulence intensity caused by a wake at particular - points in space. Laterally smeared with the Gaussian twice as wide as + points in space. Laterally smeared with the gaussian twice as wide as the wake deficit model. As recommended by Niayifar and Porte-Agel 2016 """ x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z - yc = self.centerline(x_glob) - self.y + yc = self.centerline(x_glob)[0] - self.y + zc = self.centerline(x_glob)[1] - self.z # Only compute behind rotor (x > 0) x_mask = x > 0 @@ -445,80 +452,37 @@ def wake_added_turbulence( _gaussian[comb_mask] = ( 1 - / (8 * (self.A * sigma[comb_mask]) ** 2) + / (8 * (self.WATI_Iw_multiplier * sigma[comb_mask]) ** 2) * np.exp( - -( - ((y[comb_mask] - yc[comb_mask]) ** 2 + z[comb_mask]**2) - / (2 * (self.A * sigma[comb_mask]) ** 2 * d[comb_mask]**2) + -(((y[comb_mask] - yc[comb_mask]) ** 2 + (z[comb_mask] - zc[comb_mask])**2) + / (2 * (self.WATI_Iw_multiplier * sigma[comb_mask]) ** 2 * d[comb_mask]**2) ) ) ) return _gaussian * np.nan_to_num(WATI) - # TODO: Not exactly what is given in TurbOPark, - # Check if same model by Niayifar and Porte Agel applicable - - def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: - """ - Based on principle from Shapiro, Gayme, and Meneveau, 2018, the - transverse velocity wake recovery should mirror the axial velocity - wake recovery. The centerline y position in global coordinates is - then computed numerically using the transverse velocity deficit, - based on Equation 9 in Shapiro, Gayme, and Meneveau, 2018. - """ - - _x = np.arange(0, max(xmax, 2 * dx), dx) - sigma = self._char_wake_diameter(_x) - - # Only compute behind rotor (x > 0) - x_mask = _x > 0 - - # Handle sigma = 0 nonphysical (there can't be a deficit without a char_wake_diameter) - # -> du = 0 when sigma <= 0 - sigma_mask = sigma > 0 - - # Handle sqrt does not exist -> du does not exist when sqrt does not exist - sqrt_expr = np.zeros(np.shape(sigma)) - sqrt_expr[sigma_mask] = 1 - (self.rotor_sol.Ct / self.rotor_sol.REWS ** 2) / (8 * sigma[sigma_mask] ** 2) - sqrt_mask = sqrt_expr >= 0 - - comb_mask = x_mask & sigma_mask & sqrt_mask - - dv = np.zeros(np.shape(sigma)) - - dv[comb_mask] = 1 - np.sqrt(sqrt_expr[comb_mask]) - _yc = cumulative_trapezoid(-dv, dx=dx, initial=0) - - _yc[sigma < 0] = self.y - - return _x, _yc + # [***FLAG***] Is this same model by Niayifar and Porte Agel applicable for TurbOPark? def centerline(self, x_glob: ArrayLike) -> ArrayLike: """ - Interpolates Eq. 6 from Shapiro, Gayme, and Meneveau, 2018 (same as for - GaussianWake) for centerline y position in global coordinates + No wake deflection implemented for TurbOPark wake model, as + Shapiro et al. 2018 method is likely to overpredict deflection """ + x = x_glob - self.x - yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) + # Zero wake deflection for no yaw, no tilt case + if (self.rotor_sol.yaw != 0) or (self.rotor_sol.tilt != 0): + raise NotImplementedError("No wake deflection implemented for Gaussian TurbOPark wake model.") + else: + yc = np.full(np.shape(x), self.y) + zc = np.full(np.shape(x), self.z) - return yc_temp * self.rotor_sol.v4 / self.rotor_sol.REWS + self.y + return yc, zc class TopHatTurbOParkWake(Wake): - """ - Attributes: - x: x-position of rotor in global coordinate frame - y: y-position of rotor in global coordinate frame - z: z-position of rotor in global coordinate frame - rotor_sol: Rotor solution - TIamb: Ambient turbulence intensity - A: Wake expansion calibration parameter - c_1, c_2: Empirically estimated constants in wake-added turbulence intensity - xmax: Maximum x value evaluated - dx: Interval of x values evaluated - """ def __init__( self, x: float, @@ -526,7 +490,7 @@ def __init__( z: float, rotor_sol: "RotorSolution", TIamb: float = 0.1, - A: float = 0.04, + WATI_Iw_multiplier: float = 0.04, c_1: float = 1.5, c_2: float = 0.8, xmax: float = 100.0, @@ -535,23 +499,21 @@ def __init__( self.x, self.y, self.z = x, y, z self.rotor_sol = rotor_sol self.TIamb = TIamb - self.A = A + self.WATI_Iw_multiplier = WATI_Iw_multiplier self.c_1 = c_1 self.c_2 = c_2 - # Precompute centerline far downstream - self.x_centerline, self.y_centerline = self._centerline(xmax, dx) - def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): """ Solves Eq. 6 in Pedersen et al. 2022 for wake deficit profile """ - x = x_glob - self.x - y = y_glob - self.y - z = z_glob - self.z - yc = self.centerline(x_glob) - self.y - r = np.sqrt((y - yc) ** 2 + z ** 2) + x = np.float64(x_glob - self.x) + y = np.float64(y_glob - self.y) + z = np.float64(z_glob - self.z) + yc = self.centerline(x_glob)[0] - self.y + zc = self.centerline(x_glob)[1] - self.z + r = np.sqrt((y - yc) ** 2 + (z - zc) ** 2) d = self._wake_diameter(x) @@ -562,15 +524,28 @@ def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): du[x_mask] = self._du(x[x_mask], d[x_mask]) - # Set deficit to be only inside cone of rotor - deficit = np.zeros(np.shape(y)) + # Set deficit to be only inside wake region + # Handle x is a constant + if np.shape(x) == (): + deficit = np.zeros(np.shape(y)) - for ii, _ in np.ndenumerate(deficit): - if (np.sqrt(r[ii] ** 2) <= 1/2*d[ii]) and (x[ii] >= 0): - deficit[ii] = du[ii] - else: - # Outside of wake, there is no deficit (velocity is the same as u_inf) - deficit[ii] = 0 + for ii, _ in np.ndenumerate(deficit): + if (np.sqrt(r[ii] ** 2) <= 1/2*d) and (x >= 0): + deficit[ii] = du + else: + # Outside of wake, there is no deficit (velocity is the same as u_inf) + deficit[ii] = 0 + + # Handle x is an array + else: + deficit = np.zeros(np.shape(x)) + + for ii, _ in np.ndenumerate(deficit): + if (np.sqrt(r[ii] ** 2) <= 1/2*d[ii]) and (x[ii] >= 0): + deficit[ii] = du[ii] + else: + # Outside of wake, there is no deficit (velocity is the same as u_inf) + deficit[ii] = 0 return deficit @@ -633,7 +608,7 @@ def _wake_diameter(self, x) -> ArrayLike: d = np.zeros(np.shape(x)) - d[comb_mask] = 1 + self.A * self.TIamb / beta * ( + d[comb_mask] = 1 + self.WATI_Iw_multiplier * self.TIamb / beta * ( np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1) - np.sqrt(1 + alpha ** 2) - np.log(log_expr[comb_mask]) @@ -670,48 +645,28 @@ def wake_added_turbulence( # Placeholder of zeroes for now return np.zeros(np.shape(x_glob)) - - def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: - """ - Based on principle from Shapiro, Gayme, and Meneveau, 2018, the - transverse velocity wake recovery should mirror the axial velocity - wake recovery. The centerline y position in global coordinates is - then computed numerically using the transverse velocity deficit, - based on Equation 9 in Shapiro, Gayme, and Meneveau, 2018. - """ - - _x = np.arange(0, max(xmax, 2 * dx), dx) - d = self._wake_diameter(_x) - - # Handle d = 0 -> centerline is self.y when d = 0 - d_mask = d > 0 - dv = np.zeros(np.shape(d)) - - # TODO: Not sure if this would be the right way to handle the proportionality - # between deficit in v and deficit in u - dv[d_mask] = (1 - np.tan(self.rotor_sol.yaw) * np.sqrt(1 - self.rotor_sol.Ct / self.rotor_sol.REWS ** 2)) / (d[d_mask]**2) - _yc = cumulative_trapezoid(-dv, dx=dx, initial=0) - - _yc[d < 0] = self.y - - return _x, _yc - def centerline(self, x_glob: ArrayLike) -> ArrayLike: """ - Interpolates Eq. 6 from Shapiro, Gayme, and Meneveau, 2018 (same as for - GaussianWake) for centerline y position in global coordinates + No wake deflection implemented for TurbOPark wake model, as + Shapiro et al. 2018 method is likely to overpredict deflection """ + x = x_glob - self.x - yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) + # Zero wake deflection for no yaw, no tilt case + if (self.rotor_sol.yaw != 0) or (self.rotor_sol.tilt != 0): + raise NotImplementedError("No wake deflection implemented for top hat TurbOPark wake model.") + else: + yc = np.full(np.shape(x), self.y) + zc = np.full(np.shape(x), self.z) - return yc_temp * self.rotor_sol.v4 / self.rotor_sol.REWS + self.y + return yc, zc class TurbOParkWakeModel(WakeModel): def __init__( self, - A: float = 0.04, + WATI_Iw_multiplier: float = 0.04, c_1: float = 1.5, c_2: float = 0.8, xmax: float = 100.0, @@ -720,7 +675,7 @@ def __init__( ): self.xmax = xmax self.dx = dx - self.A = A + self.WATI_Iw_multiplier = WATI_Iw_multiplier self.c_1 = c_1 self.c_2 = c_2 self.gaussian_profile = gaussian_profile @@ -744,7 +699,7 @@ def __call__( TIamb = TIamb, xmax = self.xmax, dx = self.dx, - A = self.A, + WATI_Iw_multiplier = self.WATI_Iw_multiplier, c_1 = self.c_1, c_2 = self.c_2, ) @@ -759,13 +714,38 @@ def __call__( TIamb = TIamb, xmax = self.xmax, dx = self.dx, - A = self.A, + WATI_Iw_multiplier = self.WATI_Iw_multiplier, c_1 = self.c_1, c_2 = self.c_2, ) class GaussianWake(Wake): + """ + Gaussian wake model using the derivation from Shapiro et al. (2018) + with a wake added turbulence model from Crespo and Hernandez (1996). + + The model uses a constant wake spreading rate (kw) and a constant standard + deviation (sigma) for the Gaussian profile. + + __init__: + - Args + - sigma: float, standard deviation of the Gaussian profile (default: 0.25) + - kw: float, wake spreading rate (default: 0.07) + - WATI_sigma_multiplier: float, multiplier for the wake added + turbulence intensity (WATI) sigma (default: 1.0) + - xmax: float, maximum x value for the centerline calculation (default: 100.0) + + __call__: function to create a GaussianWake instance called by the wake model solver + - Args + - x: float, x-coordinate of the turbine + - y: float, y-coordinate of the turbine + - z: float, z-coordinate of the turbine + - rotor_sol: RotorSolution, solution object containing rotor parameters + - TIamb: float, ambient turbulence intensity (default: None) + - Returns: + - GaussianWake instance with the specified parameters. + """ def __init__( self, x: float, @@ -786,33 +766,39 @@ def __init__( self.TIamb = TIamb or 0.0 # precompute centerline far downstream - self.x_centerline, self.y_centerline = self._centerline(xmax, dx) + self.x_centerline, self.yz_centerline = self._centerline(xmax, dx) def __repr__(self): return f"GaussianWake(x={self.x}, y={self.y}, z={self.z}, sigma={self.sigma}, kw={self.kw})" def _centerline(self, xmax: float, dx: float = 0.05) -> ArrayLike: """ - Solves Eq. C4. Returns centerline y position in global coordinates. + Returns centerline y/z position in rotor coordinates without constant v4/w4 factor. + Eqs. C3 and C4 in Heck et al 2023 Appendix C / Eqs. 3.8-3.9 in Shapiro et al 2018. """ - _x = np.arange(0, max(xmax, 2 * dx), dx) d = self._wake_diameter(_x) + # v and w deficit values at _x, without constant factors of v4 and w4 (see equation C3) + dvw = -0.5 / d**2 * (1 + erf(_x / (np.sqrt(2) / 2))) + # y and z centerline without without constant factors of v4 and w4 (see integration in C4) + _yzc = cumulative_trapezoid(-dvw, dx=dx, initial=0) - dv = -0.5 / d**2 * (1 + erf(_x / (np.sqrt(2) / 2))) - _yc = cumulative_trapezoid(-dv, dx=dx, initial=0) - - return _x, _yc + return _x, _yzc def centerline(self, x_glob: ArrayLike) -> ArrayLike: """ - Solves Eq. C4. Returns centerline y position in global coordinates. + Returns centerline y position in global coordinates - scales and translates the results of the `_centerline` function. + Eq. C4. in Heck et al 2023 Appendix C / Eq. 3.9 in Shapiro et al 2018. """ x = x_glob - self.x - yc_temp = np.interp(x, self.x_centerline, self.y_centerline, left=0) - - return yc_temp * self.rotor_sol.v4 / self.rotor_sol.REWS + self.y + yzc_temp = np.interp(x, self.x_centerline, self.yz_centerline, left=0) + # scale and translate centerlines + # TODO: Should this be dimensionalized with respect to the freestream + # as opposed to by REWS? (Multiply by REWS) + yc = yzc_temp * self.rotor_sol.v4 + self.y + zc = yzc_temp * self.rotor_sol.w4 + self.z + return yc, zc def centerline_wake_added_turb(self, x: ArrayLike) -> ArrayLike: """ @@ -836,50 +822,54 @@ def centerline_wake_added_turb(self, x: ArrayLike) -> ArrayLike: def _wake_diameter(self, x: ArrayLike) -> ArrayLike: """ - Solves the normalized far-wake diameter (between C1 and C2) + Solves the normalized far-wake diameter. Note that this is non-dimensionalized by D and x is actually x/D. + Defined in text between Eq. C1 and C2 in Heck et al 2023 Appendix C / in text between Eq. 3.4-3.5 in Shapiro et al 2018. + + Note that this is non-dimensionalized by D and x is actually x/D. """ return 1 + self.kw * np.log(1 + np.exp(2 * (x - 1))) def _du(self, x: ArrayLike, wake_diameter: Optional[float] = None) -> ArrayLike: """ - Solves Eq. C2 + Calcualtes the streamwise velocity deficit. + Eq. C2. in Heck et al 2023 Appendix C / Eq. 3.7 in Shapiro et al 2018. """ d = self._wake_diameter(x) if wake_diameter is None else wake_diameter - du = 0.5 * (1 - self.rotor_sol.u4) / d**2 * (1 + erf(x / (np.sqrt(2) / 2))) + du = 0.5 * (1 - (self.rotor_sol.u4 / self.rotor_sol.REWS)) / d**2 * (1 + erf(x / (np.sqrt(2) / 2))) return du - def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob=0) -> ArrayLike: + def _gaussian(self, x_glob, y_glob, z_glob, sigma_multiplier = 1): """ - Solves Eq. C1 + Defines a Gaussian profile centered at the turbine. + Gaussian portion of Eq. C1 in Heck et al 2023 Appendix C / Eq. 3.5 in Shapiro et al 2018 (without streamwise velocity deficit multiplier). """ - x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z + # calculate the centerline + yc, zc = self.centerline(x_glob) + # transform coordinates to be in turbine/wake centerline frame of reference + x = x_glob - self.x + y = y_glob - yc + z = z_glob - zc + # find wake diameter d = self._wake_diameter(x) - yc = self.centerline(x_glob) - self.y - du = self._du(x, wake_diameter=d) + # calculate sigma + sigma = self.sigma * sigma_multiplier + # calculate gaussian gaussian_ = ( 1 - / (8 * self.sigma**2) - * np.exp(-(((y - yc) ** 2 + z**2) / (2 * self.sigma**2 * d**2))) + / (8 * sigma**2) + * np.exp(-((y** 2 + z**2) / (2 * sigma**2 * d**2))) ) + return gaussian_, x, d - return gaussian_ * du - - def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob=0) -> ArrayLike: + + def deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob=0) -> ArrayLike: """ - Solves Eq. C1 where the wake deficit is defined relative to the - incident rotor wind speed following Niayifar (2016) Energies. + Distributes average streamwise velocity deficit using a Gaussian profile centered at the turbine. + Eq. C1 in Heck et al 2023 Appendix C / Eq. 3.5 in Shapiro et al 2018. """ - x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z - d = self._wake_diameter(x) - yc = self.centerline(x_glob) - self.y - du = 0.5 * (self.rotor_sol.REWS - self.rotor_sol.u4) / d**2 * (1 + erf(x / (np.sqrt(2) / 2))) - gaussian_ = ( - 1 - / (8 * self.sigma**2) - * np.exp(-(((y - yc) ** 2 + z**2) / (2 * self.sigma**2 * d**2))) - ) - + gaussian_, x, d = self._gaussian(x_glob, y_glob, z_glob) + du = self._du(x, wake_diameter=d) return gaussian_ * du def wake_added_turbulence( @@ -890,32 +880,18 @@ def wake_added_turbulence( points in space. Laterally smeared with the gaussian twice as wide as the wake deficit model. As recommended by Niayifar and Porte-Agel 2016 """ - x, y, z = x_glob - self.x, y_glob - self.y, z_glob - self.z - d = self._wake_diameter(x) - yc = self.centerline(x_glob) - self.y + gaussian_, x, _ = self._gaussian(x_glob, y_glob, z_glob, sigma_multiplier = self.WATI_sigma_multiplier) WATI = self.centerline_wake_added_turb(x) + return gaussian_ * np.nan_to_num(WATI) - _gaussian = ( - 1 - / (8 * (self.WATI_sigma_multiplier * self.sigma) ** 2) - * np.exp( - -( - ((y - yc) ** 2 + z**2) - / (2 * (self.WATI_sigma_multiplier * self.sigma) ** 2 * d**2) - ) - ) - ) - - return _gaussian * np.nan_to_num(WATI) - - def line_deficit(self, x: np.array, y: np.array): + def line_deficit(self, x: np.array, y: np.array, z = 0): """ Returns the deficit at hub height averaged along a lateral line of length 1, centered at (x, y). """ - + warnings.warn("Line deficit is deprecated as it cannot handle a z-offset. Use another deficit function.", DeprecationWarning) d = self._wake_diameter(x) - yc = self.centerline(x) + yc, zc = self.centerline(x) du = self._du(x, wake_diameter=d) erf_plus = erf((y + 0.5 - yc) / (np.sqrt(2) * self.sigma * d)) @@ -993,4 +969,4 @@ def __call__( TIamb=TIamb, xmax=self.xmax, WATI_sigma_multiplier=self.WATI_sigma_multiplier, - ) + ) \ No newline at end of file From 7fc505556ae322075e0d7db67226dd5b01b2f5e6 Mon Sep 17 00:00:00 2001 From: dengl Date: Sun, 14 Sep 2025 13:06:38 -0400 Subject: [PATCH 7/9] Added comments and renamed A parameter --- mitwindfarm/Wake.py | 57 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 43 insertions(+), 14 deletions(-) diff --git a/mitwindfarm/Wake.py b/mitwindfarm/Wake.py index b492dd0..ae6400d 100644 --- a/mitwindfarm/Wake.py +++ b/mitwindfarm/Wake.py @@ -46,7 +46,8 @@ class JensenWake(Wake): y: y-position of rotor in global coordinate frame z: z-position of rotor in global coordinate frame rotor_sol: Rotor solution - sigma: Proportionality constant for wake diameter used in Gaussian + sigma: Proportionality constant for wake diameter used in Gaussian and thus expected + by WakeModel class, even though Jensen doesn't require a sigma kw: Constant coefficient for wake growth xmax: Maximum x value evaluated dx: Interval of x values evaluated @@ -189,6 +190,7 @@ def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0) -> # Calculate deficit deficit = self.deficit(x_glob, y_glob, z_glob) + # Multiply by REWS to get Niayifar deficit return self.rotor_sol.REWS * deficit @@ -227,6 +229,19 @@ def __call__( class TurbOParkWake(Wake): + """ + Attributes: + x: x-position of rotor in global coordinate frame + y: y-position of rotor in global coordinate frame + z: z-position of rotor in global coordinate frame + rotor_sol: Rotor solution + TIamb: Ambient turbulence intensity + A: Wake expansion calibration parameter + c_1, c_2: Empirically estimated constants in wake-added turbulence intensity + xmax: Maximum x value evaluated + dx: Interval of x values evaluated + """ + def __init__( self, x: float, @@ -234,7 +249,7 @@ def __init__( z: float, rotor_sol: "RotorSolution", TIamb: float = 0.1, - WATI_Iw_multiplier: float = 0.04, + A: float = 0.04, c_1: float = 1.5, c_2: float = 0.8, xmax: float = 100.0, @@ -243,7 +258,7 @@ def __init__( self.x, self.y, self.z = x, y, z self.rotor_sol = rotor_sol self.TIamb = TIamb - self.WATI_Iw_multiplier = WATI_Iw_multiplier + self.A = A self.c_1 = c_1 self.c_2 = c_2 @@ -306,6 +321,7 @@ def niayifar_deficit(self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob = 0): # Calculate deficit deficit = self.deficit(x_glob, y_glob, z_glob) + # Multiply by REWS to get Niayifar deficit return self.rotor_sol.REWS * deficit def _centerline_wake_added_turb(self, x: ArrayLike): @@ -357,7 +373,7 @@ def _char_wake_diameter(self, x: ArrayLike) -> ArrayLike: sigma = np.zeros(np.shape(x)) - sigma[comb_mask] = epsilon + self.WATI_Iw_multiplier * self.TIamb / beta * ( + sigma[comb_mask] = epsilon + self.A * self.TIamb / beta * ( np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1) - np.sqrt(1 + alpha ** 2) - np.log(log_expr[comb_mask]) @@ -452,10 +468,10 @@ def wake_added_turbulence( _gaussian[comb_mask] = ( 1 - / (8 * (self.WATI_Iw_multiplier * sigma[comb_mask]) ** 2) + / (8 * (self.A * sigma[comb_mask]) ** 2) * np.exp( -(((y[comb_mask] - yc[comb_mask]) ** 2 + (z[comb_mask] - zc[comb_mask])**2) - / (2 * (self.WATI_Iw_multiplier * sigma[comb_mask]) ** 2 * d[comb_mask]**2) + / (2 * (self.A * sigma[comb_mask]) ** 2 * d[comb_mask]**2) ) ) ) @@ -483,6 +499,19 @@ def centerline(self, x_glob: ArrayLike) -> ArrayLike: class TopHatTurbOParkWake(Wake): + """ + Attributes: + x: x-position of rotor in global coordinate frame + y: y-position of rotor in global coordinate frame + z: z-position of rotor in global coordinate frame + rotor_sol: Rotor solution + TIamb: Ambient turbulence intensity + A: Wake expansion calibration parameter + c_1, c_2: Empirically estimated constants in wake-added turbulence intensity + xmax: Maximum x value evaluated + dx: Interval of x values evaluated + """ + def __init__( self, x: float, @@ -490,7 +519,7 @@ def __init__( z: float, rotor_sol: "RotorSolution", TIamb: float = 0.1, - WATI_Iw_multiplier: float = 0.04, + A: float = 0.6, c_1: float = 1.5, c_2: float = 0.8, xmax: float = 100.0, @@ -499,7 +528,7 @@ def __init__( self.x, self.y, self.z = x, y, z self.rotor_sol = rotor_sol self.TIamb = TIamb - self.WATI_Iw_multiplier = WATI_Iw_multiplier + self.A = A self.c_1 = c_1 self.c_2 = c_2 @@ -608,7 +637,7 @@ def _wake_diameter(self, x) -> ArrayLike: d = np.zeros(np.shape(x)) - d[comb_mask] = 1 + self.WATI_Iw_multiplier * self.TIamb / beta * ( + d[comb_mask] = 1 + self.A * self.TIamb / beta * ( np.sqrt((alpha + beta * x[comb_mask]) ** 2 + 1) - np.sqrt(1 + alpha ** 2) - np.log(log_expr[comb_mask]) @@ -642,7 +671,7 @@ def wake_added_turbulence( self, x_glob: ArrayLike, y_glob: ArrayLike, z_glob=0 ) -> ArrayLike: - # Placeholder of zeroes for now + # Placeholder of zeroes return np.zeros(np.shape(x_glob)) def centerline(self, x_glob: ArrayLike) -> ArrayLike: @@ -666,7 +695,7 @@ def centerline(self, x_glob: ArrayLike) -> ArrayLike: class TurbOParkWakeModel(WakeModel): def __init__( self, - WATI_Iw_multiplier: float = 0.04, + A: float = 0.04, c_1: float = 1.5, c_2: float = 0.8, xmax: float = 100.0, @@ -675,7 +704,7 @@ def __init__( ): self.xmax = xmax self.dx = dx - self.WATI_Iw_multiplier = WATI_Iw_multiplier + self.A = A self.c_1 = c_1 self.c_2 = c_2 self.gaussian_profile = gaussian_profile @@ -699,7 +728,7 @@ def __call__( TIamb = TIamb, xmax = self.xmax, dx = self.dx, - WATI_Iw_multiplier = self.WATI_Iw_multiplier, + A = self.A, c_1 = self.c_1, c_2 = self.c_2, ) @@ -714,7 +743,7 @@ def __call__( TIamb = TIamb, xmax = self.xmax, dx = self.dx, - WATI_Iw_multiplier = self.WATI_Iw_multiplier, + A = self.A, c_1 = self.c_1, c_2 = self.c_2, ) From 768f6c3412a25ccef39993ec439bace925e48984 Mon Sep 17 00:00:00 2001 From: dengl Date: Sun, 14 Sep 2025 13:07:20 -0400 Subject: [PATCH 8/9] Added Jensen and TurbOPark wake model descriptions --- examples/MITWindfarm_quickstart.ipynb | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/examples/MITWindfarm_quickstart.ipynb b/examples/MITWindfarm_quickstart.ipynb index 43243ab..19ad61a 100644 --- a/examples/MITWindfarm_quickstart.ipynb +++ b/examples/MITWindfarm_quickstart.ipynb @@ -223,24 +223,31 @@ "source": [ "### Wake Model\n", "\n", - "Within the package, there is a `WakeModel` abstract class with two concrete types: the `GaussianWakeModel` and the `VariableKwGaussianWakeModel`. \n", + "Within the package, there is a `WakeModel` abstract class with four concrete types: `GaussianWakeModel`, `VariableKwGaussianWakeModel`, `JensenWakeModel`, and `TurbOParkWakeModel`. \n", "\n", "The `GaussianWakeModel` has the following optional arguments (and default values): `sigma` (0.25), `kw` (0.07), `WATI_sigma_multiplier` (1.0), and `xmax` (100.0).\n", "\n", - "On the other hand, the `VariableKwGaussianWakeModel` adjusts the wake spreading rate based on the\n", + "The `VariableKwGaussianWakeModel` adjusts the wake spreading rate based on the\n", "`Ctprime` and the `TI` experienced by the wake-generating turbine based on the linear relationship: $kw = a * TI + b * C_T' + c$. Therefore, `VariableKwGaussianWakeModel` has the following optional arguments (and default values): `a`, `b`, `c`, `sigma` ($\\frac{1}{\\sqrt{8}}$), `WATI_sigma_multiplier`(1.0), and `xmax` (100.0).\n", "\n", - "If neither of these is explicitly provided to the `WindFarm` model, the `GaussianWakeModel` is used." + "The `JensenWakeModel` has one optional argument `kw` (0.07).\n", + "\n", + "The `TurbOParkWakeModel` has two versions: the newer Gaussian TurbOPark model proposed in Pedersen 2022 (`TurbOParkWake`) and the older top-hat TurbOPark model proposed in Nygaard 2020 (`TopHatTurbOParkWake`). The Gaussian TurbOPark model is used by default. Both models have the following optional arguments and default values: `A` (0.04), `TIamb` (0.06), `c_1` (1.5), `c_2` (0.8), and `xmax` (100.0).\n", + "\n", + "If no wake model is explicitly provided to the `WindFarm` model, the `GaussianWakeModel` is used." ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gaussian_wake_model = GaussianWakeModel()\n", - "variable_gaussian_wake_model = VariableKwGaussianWakeModel(a = 1, b = 1, c = 1)" + "variable_gaussian_wake_model = VariableKwGaussianWakeModel(a = 1, b = 1, c = 1)\n", + "jensen_wake_model = JensenWakeModel()\n", + "gaussian_turbopark_wake_model = TurbOParkWakeModel()\n", + "tophat_turbopark_wake_model = TurbOParkWakeModel(gaussian_profile = False, A = 0.6)" ] }, { From 14decf1b2a00b81f5de85a2d9e6fd030a2885e9b Mon Sep 17 00:00:00 2001 From: dengl Date: Sun, 21 Sep 2025 13:47:49 -0400 Subject: [PATCH 9/9] Fixed top hat TurbOPark import in initialization --- mitwindfarm/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mitwindfarm/__init__.py b/mitwindfarm/__init__.py index cbfee3e..1d08f50 100644 --- a/mitwindfarm/__init__.py +++ b/mitwindfarm/__init__.py @@ -2,6 +2,6 @@ from .Rotor import RotorSolution, AD, UnifiedAD, BEM, CosineRotor from .RotorGrid import Point, Line, Area from .Superposition import Linear, Niayifar, Quadratic, Dominant -from .Wake import WakeModel, GaussianWakeModel, GaussianWake, VariableKwGaussianWakeModel, JensenWake, JensenWakeModel, TurbOParkWake, TurbOParkWakeModel +from .Wake import WakeModel, GaussianWakeModel, GaussianWake, VariableKwGaussianWakeModel, JensenWake, JensenWakeModel, TurbOParkWake, TopHatTurbOParkWake, TurbOParkWakeModel from .windfarm import WindfarmSolution, PartialWindfarmSolution, Windfarm, CosineWindfarm from .Windfield import Uniform, PowerLaw, Superimposed