possible bug in vortex model
I was using the pywake implementation for a benchmark for an implementation. Some issues:
The deficit seems to not symmetric:
Code for generating the plot below:
from py_wake.deficit_models import VortexCylinder
if __name__ == '__main__':
import numpy as np
import matplotlib.pyplot as plt
import time
bm_PW = VortexCylinder()
WS, D_src, ct = 10, 75, 0.8
print("Across centerline")
dw = 2*D_src
cw_array = np.arange(-3*D_src, 3*D_src, 1)
deficit_PW = []
for cw in cw_array:
deficit_PW.append(bm_PW.calc_deficit(np.array([[[WS]]]),
np.array([[D_src]]),
np.array([[[[-dw]]]]),
np.array([[[[cw]]]]),
np.array([[[ct]]]))[0, 0, 0, 0])
plt.figure()
plt.plot(cw_array/D_src, deficit_PW, label='pywake')
plt.legend()
plt.xlabel("relative crosswind distance [-]")
plt.ylabel("wind speed deficit [m/s]")
plt.show()
I presume this is due to the coordinate system in the Vortex model. r is always positive, whereas cw_x can also be negative.
The issue is resolved when adding
cw_ijlk = np.abs(cw_ijlk)
to the beginning of VortexCylinder.calc_deficit.
A second issues seems to be the definition of:
# radial distance from turbine centre
r_ijlk = np.hypot(dw_ijlk, cw_ijlk)
the wiz package (where this code comes from), uses
r = np.sqrt(Xcp ** 2 + Ycp ** 2)
x and y are defined as the coordinates in the rotor plane. Whereas dw_x is perpendicular to the rotor plane.
r_ijlk is only used here:
ir = (np.abs(r_ijlk / R_il[:, na, :, na] - 1.) < self.limiter) & (np.abs(dw_ijlk / R_il[:, na, :, na]) < self.limiter)
I am not sure if this is an issue. As it is now this would be the point one radius upstream and downstream of the center of the rotor.
/Best Tobias