From 0d875ad2cd17ed40680defade7845ed49b70a9a8 Mon Sep 17 00:00:00 2001
From: Alex <alrf@dtu.dk>
Date: Thu, 8 Jul 2021 11:22:13 +0200
Subject: [PATCH] added RathmannScaled from paper and TurboPark as well as
tests
codestyle and test crash fix
---
py_wake/deficit_models/__init__.py | 2 +-
py_wake/deficit_models/noj.py | 114 ++++++++--
py_wake/deficit_models/rathmann.py | 201 +++++++++++++++++-
.../test_blockage_models.py | 11 +-
.../test_deficit_models.py | 10 +-
py_wake/tests/test_utils/test_model_utils.py | 5 +-
py_wake/turbulence_models/stf.py | 5 +-
7 files changed, 317 insertions(+), 31 deletions(-)
diff --git a/py_wake/deficit_models/__init__.py b/py_wake/deficit_models/__init__.py
index 68ae1d2cc..7d062657c 100644
--- a/py_wake/deficit_models/__init__.py
+++ b/py_wake/deficit_models/__init__.py
@@ -5,7 +5,7 @@ from .rankinehalfbody import RankineHalfBody
from .vortexcylinder import VortexCylinder
from .vortexdipole import VortexDipole
from .rathmann import Rathmann
-from .noj import NOJDeficit, NOJLocalDeficit
+from .noj import NOJDeficit, NOJLocalDeficit, TurboNOJDeficit
from .gaussian import BastankhahGaussianDeficit, IEA37SimpleBastankhahGaussianDeficit, \
NiayifarGaussianDeficit, ZongGaussianDeficit
from .fuga import FugaDeficit, FugaYawDeficit
diff --git a/py_wake/deficit_models/noj.py b/py_wake/deficit_models/noj.py
index 92516950c..0987b3902 100644
--- a/py_wake/deficit_models/noj.py
+++ b/py_wake/deficit_models/noj.py
@@ -1,7 +1,5 @@
from numpy import newaxis as na
-
import numpy as np
-from py_wake.deficit_models import DeficitModel
from py_wake.superposition_models import SquaredSum, LinearSum
from py_wake.wind_farm_models.engineering_models import PropagateDownwind
from py_wake.utils.area_overlapping_factor import AreaOverlappingFactor
@@ -21,10 +19,10 @@ class NOJDeficit(NiayifarGaussianDeficit, AreaOverlappingFactor):
def _calc_layout_terms(self, WS_ilk, WS_eff_ilk, D_src_il, D_dst_ijl, dw_ijlk, cw_ijlk, **kwargs):
WS_ref_ilk = (WS_ilk, WS_eff_ilk)[self.use_effective_ws]
R_src_il = D_src_il / 2
- wake_radius, k_ijlk = self._wake_radius(D_src_il, dw_ijlk, **kwargs)
- term_denominator_ijlk = (1 + k_ijlk * dw_ijlk / R_src_il[:, na, :, na])**2
+ wake_radius_ijlk = self.wake_radius(D_src_il, dw_ijlk, **kwargs)
+ term_denominator_ijlk = (wake_radius_ijlk / R_src_il[:, na, :, na])**2
term_denominator_ijlk += (term_denominator_ijlk == 0)
- A_ol_factor_ijlk = self.overlapping_area_factor(wake_radius, dw_ijlk, cw_ijlk, D_src_il, D_dst_ijl)
+ A_ol_factor_ijlk = self.overlapping_area_factor(wake_radius_ijlk, dw_ijlk, cw_ijlk, D_src_il, D_dst_ijl)
with np.warnings.catch_warnings():
np.warnings.filterwarnings('ignore', r'invalid value encountered in true_divide')
@@ -32,18 +30,16 @@ class NOJDeficit(NiayifarGaussianDeficit, AreaOverlappingFactor):
def calc_deficit(self, WS_ilk, WS_eff_ilk, D_src_il, D_dst_ijl, dw_ijlk, cw_ijlk, ct_ilk, **kwargs):
if not self.deficit_initalized:
+ kwargs['ct_ilk'] = ct_ilk
self._calc_layout_terms(WS_ilk, WS_eff_ilk, D_src_il, D_dst_ijl, dw_ijlk, cw_ijlk, **kwargs)
ct_ilk = np.minimum(ct_ilk, 1) # treat ct_ilk for np.sqrt()
term_numerator_ilk = (1 - np.sqrt(1 - ct_ilk))
return term_numerator_ilk[:, na] * self.layout_factor_ijlk
- def _wake_radius(self, D_src_il, dw_ijlk, **kwargs):
+ def wake_radius(self, D_src_il, dw_ijlk, **kwargs):
k_ijlk = np.atleast_3d(self.k_ilk(kwargs.get('TI_eff_ilk', 0)))[:, na]
wake_radius_ijlk = (k_ijlk * dw_ijlk + D_src_il[:, na, :, na] / 2)
- return wake_radius_ijlk, k_ijlk
-
- def wake_radius(self, D_src_il, dw_ijlk, **kwargs):
- return self._wake_radius(D_src_il, dw_ijlk, **kwargs)[0]
+ return wake_radius_ijlk
def calc_deficit_convection(self, WS_ilk, WS_eff_ilk, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, **kwargs):
raise NotImplementedError("calc_deficit_convection not implemented for NOJ")
@@ -129,6 +125,57 @@ class NOJLocal(PropagateDownwind):
turbulenceModel=turbulenceModel)
+class TurboNOJDeficit(NOJDeficit, AreaOverlappingFactor):
+ """
+ Modified definition of the wake expansion given by Nygaard [1], which
+ assumes the wake expansion rate to be proportional to the local turbulence
+ intensity in the wake. Here the local turbulence intensity is defined as
+ the combination of ambient and wake added turbulence. Using the added
+ wake turbulence model by Frandsen and integrating, an analytical expression
+ for the wake radius can be obtained.
+ The definition in [1] of ambient turbulence is the free-stream TI and for
+ this the model constant A has been tuned, however a fully consistent
+ formulation of the model should probably use the locally effective TI,
+ which includes the added turbulence from upstream turbines.
+ [1] Nygaard 2020 J. Phys.: Conf. Ser. 1618 062072
+ https://doi.org/10.1088/1742-6596/1618/6/062072
+ """
+ args4deficit = ['WS_ilk', 'WS_eff_ilk', 'D_src_il', 'D_dst_ijl',
+ 'dw_ijlk', 'cw_ijlk', 'ct_ilk', 'TI_ilk', 'TI_eff_ilk']
+
+ def __init__(self, A=.6, cTI=[1.5, 0.8], use_effective_ws=False, use_effective_ti=False):
+ self.A = A
+ self.use_effective_ws = use_effective_ws
+ self.use_effective_ti = use_effective_ti
+ self.cTI = cTI
+
+ def wake_radius(self, D_src_il, dw_ijlk, **kwargs):
+ TI_ref_ilk = (kwargs['TI_ilk'], kwargs['TI_eff_ilk'])[
+ self.use_effective_ti]
+ ct_ilk = kwargs['ct_ilk']
+ # constants from Frandsen
+ c1, c2 = self.cTI
+
+ # constants related to ambient turbulence
+ alpha_ilk = c1 * TI_ref_ilk
+ # avoid zero division
+ ct_ilk = np.maximum(ct_ilk, 1e-20)
+ beta_ilk = c2 * TI_ref_ilk / np.sqrt(ct_ilk)
+
+ fac_ilk = self.A * TI_ref_ilk * D_src_il[..., na] / beta_ilk
+ term1_ijlk = np.sqrt(
+ (alpha_ilk[:, na] + beta_ilk[:, na] * dw_ijlk / D_src_il[:, na, :, na])**2 + 1)
+ term2_ilk = np.sqrt(1 + alpha_ilk**2)
+ term3_ijlk = (term1_ijlk + 1) * alpha_ilk[:, na]
+ term4_ijlk = (term2_ilk[:, na] + 1) * (alpha_ilk[:, na] +
+ beta_ilk[:, na] * np.abs(dw_ijlk) / D_src_il[:, na, :, na])
+
+ wake_radius_ijlk = 0.5 * (D_src_il[:, na, :, na] + fac_ilk[:, na] * (
+ term1_ijlk - term2_ilk[:, na] - np.log(term3_ijlk / term4_ijlk)))
+
+ return wake_radius_ijlk
+
+
def main():
if __name__ == '__main__':
from py_wake.examples.data.iea37._iea37 import IEA37Site
@@ -143,16 +190,25 @@ def main():
wf_model = NOJ(site, windTurbines)
wf_model_local = NOJLocal(site, windTurbines, turbulenceModel=STF2017TurbulenceModel())
-
+ wf_model_turbo = PropagateDownwind(site, windTurbines, rotorAvgModel=RotorCenter(),
+ wake_deficitModel=TurboNOJDeficit(
+ use_effective_ws=True, use_effective_ti=False),
+ superpositionModel=LinearSum(),
+ turbulenceModel=STF2017TurbulenceModel())
+ # wf_model_turbo = NOJLocal(
+ # site, windTurbines, turbulenceModel=STF2017TurbulenceModel())
# run wind farm simulation
sim_res = wf_model(x, y)
sim_res_local = wf_model_local(x, y)
+ sim_res_turbo = wf_model_turbo(x, y)
# calculate AEP
aep = sim_res.aep().sum()
aep_local = sim_res_local.aep().sum()
+ aep_turbo = sim_res_turbo.aep().sum()
# plot wake map
- fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4.5), tight_layout=True)
+ fig, (ax1, ax2, ax3) = plt.subplots(
+ 1, 3, figsize=(11, 4.5), tight_layout=True)
levels = np.arange(0, 10.5, 0.5)
print(wf_model)
flow_map = sim_res.flow_map(wd=30, ws=9.8)
@@ -167,10 +223,44 @@ def main():
flow_map.plot_windturbines(ax=ax2)
ax2.set_title('Local Jensen, AEP: %.2f GWh' % aep_local)
+ # plot wake map
+ print(wf_model_turbo)
+ flow_map = sim_res_turbo.flow_map(wd=30, ws=9.8)
+ flow_map.plot_wake_map(levels=levels, ax=ax3, plot_colorbar=False)
+ flow_map.plot_windturbines(ax=ax3)
+ ax3.set_title('Turbo Jensen, AEP: %.2f GWh' % aep_turbo)
+
plt.figure()
flow_map.plot_ti_map()
plt.title('TI map for NOJLocal with STF2017 turbulence model')
plt.show()
+ # plot wake width as in Nygaard 2020
+ D = 1
+ D_src_il = np.array([[D]])
+ x = np.linspace(0, 60, 100)
+ dw_ijlk = x[na, :, na, na]
+
+ noj = NOJDeficit(k=0.04)
+ noj_wr = noj.wake_radius(D_src_il, dw_ijlk)
+
+ ct_ilk = np.array([[[8 / 9]]]) # thrust coefficient
+ TI_ilk = np.array([[[0.06]]])
+ TI_eff_ilk = np.array([[[0.06]]])
+ tj = TurboNOJDeficit()
+ tj_wr = tj.wake_radius(
+ D_src_il, dw_ijlk, ct_ilk=ct_ilk, TI_ilk=TI_ilk, TI_eff_ilk=TI_eff_ilk)
+
+ plt.figure()
+ plt.title(
+ 'Wake width comparison, NOJ orig and TurboNOJ (Nygaard2020) TI=6%')
+ plt.plot(x, noj_wr[0, :, 0, 0], label='NOJ, k=0.04')
+ plt.plot(x, tj_wr[0, :, 0, 0], label='TurboNOJ')
+ plt.xlabel('x/D')
+ plt.ylabel('y/D')
+ plt.grid()
+ plt.legend()
+ plt.show()
+
main()
diff --git a/py_wake/deficit_models/rathmann.py b/py_wake/deficit_models/rathmann.py
index 603ce499c..9a6c651d4 100644
--- a/py_wake/deficit_models/rathmann.py
+++ b/py_wake/deficit_models/rathmann.py
@@ -5,7 +5,7 @@ from py_wake.deficit_models import BlockageDeficitModel
class Rathmann(BlockageDeficitModel):
"""
- Ole Sten Rathmann (DTU) developed in 2020 an approximation to the vortex
+ Ole Sten Rathmann (DTU) developed in 2020 an approximation [1] to the vortex
cylinder solution (E. Branlard and M. Gaunaa, 2014). In speed it is
comparable to the vortex dipole method, whilst giving a flow-field
nearly identical to the vortex cylinder model for x/R < -1. Its centreline
@@ -14,6 +14,10 @@ class Rathmann(BlockageDeficitModel):
a point upstream. To simulate the speed-up downstream the deficit is mirrored in the
rotor plane with a sign change.
+ References:
+ [1] A Meyer Forsting, OS Rathmann, MP van der Laan, N Troldborg, B Gribben, G Hawkes, E Branlard -
+ Verification of induction zone models for wind farm annual energy production estimation -
+ Journal of Physics: Conference Series 1934 (2021) 012023
"""
args4deficit = ['WS_ilk', 'D_src_il', 'dw_ijlk', 'cw_ijlk', 'ct_ilk']
@@ -96,18 +100,72 @@ class Rathmann(BlockageDeficitModel):
return deficit_ijlk
+class RathmannScaled(Rathmann):
+ """
+ Vortex cylinder based models consistently underestimate the induction, due to missing
+ wake expansion [1]. A simple fix alleviating this issue is to simply scale the results
+ with respect to the thrust coefficient as demonstrated in [1]. There is also a small
+ depenancy on the distance fron the rotor.
+
+ References:
+ [1] A Meyer Forsting, OS Rathmann, MP van der Laan, N Troldborg, B Gribben, G Hawkes, E Branlard -
+ Verification of induction zone models for wind farm annual energy production estimation -
+ Journal of Physics: Conference Series 1934 (2021) 012023
+ """
+
+ args4deficit = ['WS_ilk', 'D_src_il', 'dw_ijlk', 'cw_ijlk', 'ct_ilk']
+
+ def __init__(self, sct=1.0, limiter=1e-10, exclude_wake=True):
+ BlockageDeficitModel.__init__(self)
+ # coefficients for BEM approximation by Madsen (1997)
+ self.a0p = np.array([0.2460, 0.0586, 0.0883])
+ # limiter to avoid singularities
+ self.limiter = limiter
+ # coefficient for scaling the effective forcing
+ self.sct = sct
+ # if used in a wind farm simulation, set deficit in wake region to
+ # zero, as here the wake model is active
+ self.exclude_wake = exclude_wake
+ # scaling coefficients for Eq.11-13 in [1]
+ self.sd = np.array([1.02, 0.1554, 0.0005012, 8.45, 0.025])
+
+ def deficit_scaling(self, D_src_il, dw_ijlk, cw_ijlk, ct_ilk):
+ """
+ Scaling function defined in [1], Eq. 11-13 forcing the output closer to the
+ CFD results.
+ """
+ r = np.hypot(cw_ijlk, dw_ijlk) / D_src_il[:, na, :, na]
+ fac = 1. + self.sd[4] * (r - 5.)
+ fac[fac > 1.] = 1.
+ mval = 1. - 4. * self.sd[4]
+ fac[fac < mval] = mval
+ boost_ijlk = self.sd[0] * \
+ np.exp(self.sd[1] * fac * ct_ilk[:, na]) + \
+ self.sd[2] * np.exp(self.sd[3] * fac * ct_ilk[:, na])
+
+ return boost_ijlk
+
+ def calc_deficit(self, WS_ilk, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, **_):
+ boost_ijlk = self.deficit_scaling(D_src_il, dw_ijlk, cw_ijlk, ct_ilk)
+ return boost_ijlk * Rathmann.calc_deficit(self, WS_ilk, D_src_il, dw_ijlk, cw_ijlk, ct_ilk, **_)
+
+
def main():
if __name__ == '__main__':
from py_wake.examples.data.iea37._iea37 import IEA37Site
from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines
+ from py_wake.site._site import UniformSite
+ from py_wake.wind_turbines import OneTypeWindTurbines
from py_wake.superposition_models import LinearSum
from py_wake.wind_farm_models import All2AllIterative
from py_wake.deficit_models.no_wake import NoWakeDeficit
+ from py_wake.rotor_avg_models import RotorCenter
from py_wake.deficit_models.vortexcylinder import VortexCylinder
from py_wake.deficit_models.vortexdipole import VortexDipole
import matplotlib.pyplot as plt
from py_wake import HorizontalGrid
from timeit import default_timer as timer
+ import time
# setup site, turbines and wind farm model
site = IEA37Site(16)
@@ -115,10 +173,13 @@ def main():
windTurbines = IEA37_WindTurbines()
d = windTurbines.diameter()
ra = Rathmann()
+ ras = RathmannScaled()
grid = HorizontalGrid(x=np.linspace(-6, 6, 100) * d, y=np.linspace(0, 4, 100) * d)
noj_ra = All2AllIterative(site, windTurbines, wake_deficitModel=NoWakeDeficit(),
superpositionModel=LinearSum(), blockage_deficitModel=ra)
+ noj_ras = All2AllIterative(site, windTurbines, wake_deficitModel=NoWakeDeficit(),
+ superpositionModel=LinearSum(), blockage_deficitModel=ras)
noj_vc = All2AllIterative(site, windTurbines, wake_deficitModel=NoWakeDeficit(),
superpositionModel=LinearSum(), blockage_deficitModel=VortexCylinder())
noj_vd = All2AllIterative(site, windTurbines, wake_deficitModel=NoWakeDeficit(),
@@ -126,29 +187,40 @@ def main():
t1 = timer()
flow_map = noj_ra(x=[0], y=[0], wd=[270], ws=[10]).flow_map(grid=grid)
t2 = timer()
- flow_map_vc = noj_vc(x=[0], y=[0], wd=[270], ws=[10]).flow_map(grid=grid)
+ flow_map_ras = noj_ras(x=[0], y=[0], wd=[270], ws=[10]).flow_map(grid=grid)
t3 = timer()
- flow_map_vd = noj_vd(x=[0], y=[0], wd=[270], ws=[10]).flow_map(grid=grid)
+ flow_map_vc = noj_vc(x=[0], y=[0], wd=[270], ws=[10]).flow_map(grid=grid)
t4 = timer()
- print(t2 - t1, t3 - t2, t4 - t3)
+ flow_map_vd = noj_vd(x=[0], y=[0], wd=[270], ws=[10]).flow_map(grid=grid)
+ t5 = timer()
+ print(t2 - t1, t3 - t2, t4 - t3, t5 - t4)
plt.figure()
- clevels = np.array([.6, .7, .8, .9, .95, .98, .99, .995, .998, .999, 1., 1.005, 1.01, 1.02, 1.05]) * 10.
+ clevels = np.array(
+ [.6, .7, .8, .9, .95, .98, .99, .995, .998, .999, 1., 1.005, 1.01, 1.02, 1.05]) * 10.
flow_map.plot_wake_map(levels=clevels)
- plt.contour(flow_map.x, flow_map.y, flow_map.WS_eff[:, :, 0, -1, 0], levels=clevels, colors='k', linewidths=1)
- plt.contour(flow_map.x, flow_map.y, flow_map_vc.WS_eff[:, :, 0, -1, 0], levels=clevels, colors='r', linewidths=1, linestyles='dashed')
- plt.contour(flow_map.x, flow_map.y, flow_map_vd.WS_eff[:, :, 0, -1, 0], levels=clevels, colors='b', linewidths=1, linestyles='dotted')
+ plt.contour(flow_map.x, flow_map.y,
+ flow_map.WS_eff[:, :, 0, -1, 0], levels=clevels, colors='k', linewidths=1)
+ plt.contour(flow_map.x, flow_map.y, flow_map_ras.WS_eff[:, :, 0, -1, 0],
+ levels=clevels, colors='g', linewidths=1, linestyles='dashed')
+ plt.contour(flow_map.x, flow_map.y, flow_map_vc.WS_eff[:, :, 0, -1, 0],
+ levels=clevels, colors='r', linewidths=1, linestyles='dashed')
+ plt.contour(flow_map.x, flow_map.y, flow_map_vd.WS_eff[:, :, 0, -1, 0],
+ levels=clevels, colors='b', linewidths=1, linestyles='dotted')
+ plt.plot([0, 0], [0, 0], 'k-', label='Rathmann')
+ plt.plot([0, 0], [0, 0], 'g--', label='scaled Rathmann')
+ plt.plot([0, 0], [0, 0], 'r--', label='vortex cylinder')
+ plt.plot([0, 0], [0, 0], 'b:', label='vortex dipole')
plt.title('Rathmann')
plt.ylabel("Crosswind distance [y/R]")
plt.xlabel("Downwind distance [x/R]")
+ plt.legend()
plt.show()
# run wind farm simulation
sim_res = noj_ra(x, y, wd=[0, 30, 45, 60, 90], ws=[5, 10, 15])
-
# calculate AEP
aep = sim_res.aep().sum()
-
# plot wake map
plt.figure()
print(noj_ra)
@@ -157,5 +229,114 @@ def main():
plt.title('Rathmann model, AEP: %.3f GWh' % aep)
plt.show()
+ # run wind farm simulation
+ sim_res = noj_ras(x, y, wd=[0, 30, 45, 60, 90], ws=[5, 10, 15])
+ # calculate AEP
+ aep = sim_res.aep().sum()
+ # plot wake map
+ plt.figure()
+ print(noj_ras)
+ flow_map = sim_res.flow_map(wd=0, ws=10)
+ flow_map.plot_wake_map(levels=clevels, plot_colorbar=False)
+ plt.title('Rathmann model, AEP: %.3f GWh' % aep)
+ plt.show()
+
+ class epfl_model_wt(OneTypeWindTurbines):
+ def __init__(self):
+ OneTypeWindTurbines.__init__(self, 'NREL 5MW', diameter=2,
+ hub_height=1, ct_func=self._ct,
+ power_func=self._power, power_unit='W')
+
+ def _ct(self, u):
+ ct = 0.798
+ return ct * u / u
+
+ def _power(self, u):
+ cp = 0.5
+ A = np.pi
+ rho = 1.225
+ return 0.5 * rho * u**3 * A * cp
+
+ wt = epfl_model_wt()
+ d = wt.diameter()
+ h = wt.hub_height()
+ wt_x = np.array([-6. * d, -3. * d, 0. * d, 3. * d, 6. * d])
+ wt_y = np.array([0., 0., 0., 0., 0.])
+
+ class epfl_wt(UniformSite):
+ def __init__(self):
+ p_wd = [1]
+ ws = [1]
+ ti = [0.06]
+ UniformSite.__init__(self, p_wd=p_wd, ti=ti, ws=ws)
+ self.initial_position = np.array([wt_x, wt_y]).T
+
+ site = epfl_wt()
+ blockage_models = [
+ ('Rathmann', Rathmann()),
+ ('RathmannScaled', RathmannScaled())]
+
+ def pywake_run(blockage_models):
+
+ grid = HorizontalGrid(x=np.linspace(-10 * d, 10 * d, 100),
+ y=np.linspace(-15, 2 * d, 100))
+
+ res = {}
+ out = {}
+ for nam, blockage_model in blockage_models:
+ print(nam, blockage_model)
+ # for dum, rotor_avg_model in rotor_avg_models:
+ wm = All2AllIterative(site, wt, wake_deficitModel=NoWakeDeficit(),
+ superpositionModel=LinearSum(), blockage_deficitModel=blockage_model,
+ rotorAvgModel=RotorCenter())
+ res[nam] = wm(wt_x, wt_y, wd=[180., 195., 210., 225.], ws=[1.])
+ tic = time.perf_counter()
+ flow_map = res[nam].flow_map(grid=grid, ws=[1.], wd=225.)
+ toc = time.perf_counter()
+ elapsed_time = toc - tic
+ out[nam] = flow_map['WS_eff'] / flow_map['WS']
+ out[nam]['time'] = elapsed_time
+ return out, res
+
+ out, res = pywake_run(blockage_models)
+ # CFD data from Meyer Forsting et al. 2017
+ p_cfd = np.array([[-1.12511646713429, 0.268977884040651, 0.712062872514373, 1.08033923355738, 1.97378837188847],
+ [-0.610410399845213, 0.355339771667814,
+ 0.670255435929930, 0.915154608331424, 1.52808830519513],
+ [-0.0988002822865217, 0.451698279664756,
+ 0.636987630794206, 0.751760283763044, 1.03629687168984],
+ [0.477918399858401, 0.628396438496795, 0.663799054750132, 0.628396438496795, 0.479688468006036]])
+
+ plt.figure()
+ lines = ['.-', 'o--', '^', '-.', '.-', '.-']
+ theta = [0, 15, 30, 45]
+ viridis = plt.cm.viridis(np.linspace(0, 0.9, 5))
+ jj = 0
+ tno = np.arange(1, 6, 1)
+ for i in range(4):
+ ii = len(theta) - i - 1
+ plt.plot(tno, ((p_cfd[ii, :] / 100. + 1.) / (p_cfd[ii, 0] / 100. + 1.) - 1.) * 100, 's:', color=viridis[i], label='CFD: ' + str(theta[i]) + 'deg',
+ lw=2)
+
+ for nam, blockage_model in blockage_models:
+ for i in range(4):
+ if i == 3:
+ plt.plot(tno, (res[nam].Power[:, i, 0] - res[nam].Power[0, i, 0]) / res[nam].Power[0, i, 0] * 100, lines[jj],
+ label=nam, color=viridis[i],
+ lw=1, alpha=0.8)
+ else:
+ plt.plot(tno, (res[nam].Power[:, i, 0] - res[nam].Power[0, i, 0]) / res[nam].Power[0, i, 0] * 100, lines[jj],
+ color=viridis[i],
+ lw=1, alpha=0.8)
+ jj += 1
+
+ plt.grid(alpha=0.2)
+ plt.xlabel('Turbine no.')
+ plt.ylabel('Power change, $(P-P_1)/P_1$ [%]')
+ plt.xticks([0, 1, 2, 3, 4, 5])
+ plt.xlim([.9, 5.1])
+ plt.legend(fontsize=11)
+ plt.show()
+
main()
diff --git a/py_wake/tests/test_deficit_models/test_blockage_models.py b/py_wake/tests/test_deficit_models/test_blockage_models.py
index 1ee3e4986..fb4b9a5ac 100644
--- a/py_wake/tests/test_deficit_models/test_blockage_models.py
+++ b/py_wake/tests/test_deficit_models/test_blockage_models.py
@@ -13,7 +13,7 @@ from py_wake.deficit_models.selfsimilarity import SelfSimilarityDeficit, SelfSim
from py_wake.deficit_models.vortexdipole import VortexDipole
from py_wake.deficit_models.rankinehalfbody import RankineHalfBody
from py_wake.deficit_models.hybridinduction import HybridInduction
-from py_wake.deficit_models.rathmann import Rathmann
+from py_wake.deficit_models.rathmann import Rathmann, RathmannScaled
from py_wake.flow_map import XYGrid
debug = False
@@ -49,6 +49,9 @@ def setup():
(Rathmann,
[9.950432, 9.926598, 9.882538, 9.792603, 9.597037, 10.041663, 10.397465, 10.204838, 10.116258, 10.072778],
[9.94641, 9.917526, 9.858037, 9.707178, 9.190401, 10.0, 10.0, 10.0, 10.0, 10.0]),
+ (RathmannScaled,
+ [9.929791, 9.897708, 9.838744, 9.719092, 9.456743, 10.056169, 10.535845, 10.277511, 10.159648, 10.101453],
+ [9.924278, 9.88539, 9.805788, 9.60523, 8.908535, 10.0, 10.0, 10.0, 10.0, 10.0]),
][::-1])
def test_blockage_map(setup, blockage_model, center_ref, side_ref):
site, windTurbines = setup
@@ -98,6 +101,9 @@ def test_blockage_map(setup, blockage_model, center_ref, side_ref):
(Rathmann,
[9.950432, 9.926598, 9.882538, 9.792603, 9.597037, 10.041663, 10.397465, 6.428758, 6.899183, 7.299176],
[9.94641, 9.917526, 9.858037, 9.707178, 9.190401, 4.560631, 5.505472, 6.223921, 6.782925, 7.226399]),
+ (RathmannScaled,
+ [9.929791, 9.897708, 9.838744, 9.719092, 9.456743, 10.056169, 10.535845, 6.501432, 6.942573, 7.327852],
+ [9.924278, 9.88539, 9.805788, 9.60523, 8.908535, 4.560631, 5.505472, 6.223921, 6.782925, 7.226399]),
][::-1])
def test_wake_and_blockage(setup, blockage_model, center_ref, side_ref):
site, windTurbines = setup
@@ -132,7 +138,8 @@ def test_wake_and_blockage(setup, blockage_model, center_ref, side_ref):
(VortexDipole, 0.4321021),
(RankineHalfBody, 0.4321021),
(HybridInduction, 0.5372064),
- (Rathmann, 0.4233034298949511)
+ (Rathmann, 0.4233034298949511),
+ (RathmannScaled, 0.5928161590163314)
][::-1])
def test_aep_two_turbines(setup, blockage_model, blockage_loss):
site, windTurbines = setup
diff --git a/py_wake/tests/test_deficit_models/test_deficit_models.py b/py_wake/tests/test_deficit_models/test_deficit_models.py
index ec925d089..7cffc3816 100644
--- a/py_wake/tests/test_deficit_models/test_deficit_models.py
+++ b/py_wake/tests/test_deficit_models/test_deficit_models.py
@@ -8,7 +8,7 @@ from py_wake.deficit_models.gaussian import BastankhahGaussianDeficit, IEA37Simp
ZongGaussianDeficit, NiayifarGaussianDeficit, BastankhahGaussian, IEA37SimpleBastankhahGaussian, ZongGaussian,\
NiayifarGaussian
from py_wake.deficit_models.gcl import GCLDeficit, GCL, GCLLocal
-from py_wake.deficit_models.noj import NOJDeficit, NOJ, NOJLocalDeficit, NOJLocal
+from py_wake.deficit_models.noj import NOJDeficit, NOJ, NOJLocalDeficit, NOJLocal, TurboNOJDeficit
from py_wake.deficit_models import VortexDipole
from py_wake.examples.data.hornsrev1 import Hornsrev1Site
from py_wake.examples.data.iea37 import iea37_path
@@ -41,6 +41,10 @@ class GCLLocalDeficit(GCLDeficit):
23558.34198, 36738.52415, 38663.44595, 21056.39764, 12042.79324,
13813.46269, 30999.42279, 63947.61202, 17180.40299, 11334.12323,
6972.14345])),
+ (TurboNOJDeficit(), (354154.2962989713, [9320.85263, 8138.29496, 10753.75662, 13398.00865, 21189.29438,
+ 24190.84895, 37081.91938, 41369.66605, 23488.54863, 12938.48451,
+ 14065.00719, 30469.75602, 71831.78532, 16886.85274, 11540.51872,
+ 7490.70156])),
(BastankhahGaussianDeficit(), (355971.9717035484,
[9143.74048, 8156.71681, 11311.92915, 13955.06316, 19807.65346,
@@ -102,6 +106,8 @@ def test_IEA37_ex16(deficitModel, aep_ref):
[3.27, 3.27, 9.0, 7.46, 7.46, 7.46, 7.46, 7.31, 7.31, 7.31, 7.31, 8.3, 8.3, 8.3, 8.3, 8.3, 8.3]),
(NOJLocalDeficit(),
[3.09, 3.09, 9., 9., 5.54, 5.54, 5.54, 5.54, 5.54, 5.54, 6.73, 6.73, 6.73, 6.73, 6.73, 6.73, 6.73]),
+ (TurboNOJDeficit(),
+ [3.09, 3.09, 9., 9., 5.54, 5.54, 5.54, 5.54, 5.54, 5.54, 6.73, 6.73, 6.73, 6.73, 6.73, 6.73, 6.73]),
(BastankhahGaussianDeficit(),
[0.18, 3.6, 7.27, 8.32, 7.61, 6.64, 5.96, 6.04, 6.8, 7.69, 8.08, 7.87, 7.59, 7.46, 7.55, 7.84, 8.19]),
(IEA37SimpleBastankhahGaussianDeficit(),
@@ -156,6 +162,8 @@ def test_deficitModel_wake_map(deficitModel, ref):
# test that the result is equal to last run (no evidens that these number are correct)
[(NOJDeficit(), [100., 75., 150., 100., 100.]),
(NOJLocalDeficit(), [71., 46., 92., 71., 61.5]),
+ (TurboNOJDeficit(), [99.024477, 61.553917,
+ 123.107833, 92.439673, 97.034049]),
(BastankhahGaussianDeficit(),
[83.336286, 57.895893, 115.791786, 75.266662, 83.336286]),
(IEA37SimpleBastankhahGaussianDeficit(),
diff --git a/py_wake/tests/test_utils/test_model_utils.py b/py_wake/tests/test_utils/test_model_utils.py
index 79c34566c..d0778e257 100644
--- a/py_wake/tests/test_utils/test_model_utils.py
+++ b/py_wake/tests/test_utils/test_model_utils.py
@@ -23,14 +23,13 @@ def test_get_models():
def test_get_signature():
assert get_signature(NOJDeficit) == "NOJDeficit(k=0.1, use_effective_ws=False)"
- print()
assert get_signature(NOJDeficit, indent_level=1) == """NOJDeficit(
k=0.1,
use_effective_ws=False)"""
assert (get_signature(STF2017TurbulenceModel) ==
- "STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=LinearSum(), weight_function=FrandsenWeight())")
+ "STF2017TurbulenceModel(c=[1.5, 0.8], addedTurbulenceSuperpositionModel=LinearSum(), weight_function=FrandsenWeight())")
assert (get_signature(STF2017TurbulenceModel, {'addedTurbulenceSuperpositionModel': SqrMaxSum}) ==
- "STF2017TurbulenceModel(addedTurbulenceSuperpositionModel=SqrMaxSum(), weight_function=FrandsenWeight())")
+ "STF2017TurbulenceModel(c=[1.5, 0.8], addedTurbulenceSuperpositionModel=SqrMaxSum(), weight_function=FrandsenWeight())")
assert(get_signature(XRSite) ==
"XRSite(ds, initial_position=None, interp_method='linear', shear=None, distance=StraightDistance(), default_ws=[3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25], bounds='check')")
diff --git a/py_wake/turbulence_models/stf.py b/py_wake/turbulence_models/stf.py
index 15b99fc06..ef31fe09a 100644
--- a/py_wake/turbulence_models/stf.py
+++ b/py_wake/turbulence_models/stf.py
@@ -59,16 +59,17 @@ class STF2017TurbulenceModel(TurbulenceModel):
args4addturb = ['dw_ijlk', 'cw_ijlk', 'D_src_il', 'ct_ilk', 'TI_ilk']
- def __init__(self, addedTurbulenceSuperpositionModel=LinearSum(),
+ def __init__(self, c=[1.5, 0.8], addedTurbulenceSuperpositionModel=LinearSum(),
weight_function=FrandsenWeight(), **kwargs):
TurbulenceModel.__init__(self, addedTurbulenceSuperpositionModel, **kwargs)
+ self.c = c
self.apply_weight = weight_function.apply_weight
def max_centre_wake_turbulence(self, dw_ijlk, cw_ijlk, D_src_il, ct_ilk, **_):
dist_ijlk = np.sqrt(dw_ijlk**2 + cw_ijlk**2) / D_src_il[:, na, :, na]
# In the standard (see page 103), the maximal added TI is calculated as
# TI_add = 1/(1.5 + 0.8*d/sqrt(Ct))
- return 1 / (1.5 + 0.8 * dist_ijlk / np.sqrt(ct_ilk)[:, na])
+ return 1 / (self.c[0] + self.c[1] * dist_ijlk / np.sqrt(ct_ilk)[:, na])
def calc_added_turbulence(self, dw_ijlk, cw_ijlk, D_src_il, TI_ilk, **kwargs):
""" Calculate the added turbulence intensity at locations specified by
--
GitLab