Discrepancies: Verification of FLORIS results with WAsP and manual calculation

[niajhsaka]

I am trying to verify the AEP calculation of FLORIS with WAsP and through manual calculation in a spreadsheet, however I am getting a difference of about 4%-5% between FLORIS and WAsP. WAsP and the manual calculation match closely.

These are the conditions that I am simulating in WAsP:

1. Elevation – uniform elevation of 0 m
2. Roughness – Uniform roughness of 0.1 m. The roughness both inside and outside the site map boundary is set to 0.1 m.
3. Resource grid calculated at 100m
4. Single turbine - Turbine hub height used is 100 m, rotor diameter is 156 m and rated power is 3.3 MW.
5. Reference Air density at which power production curve is defined – 1.225 kg/m3
6. Air density at the site – 1.225 kg/m3. This is ensured by turning off the setting – “Do interpolate/ extrapolate power curves to match air density”
7. Generalized wind climate (GWC) file for the site is obtained from Global Wind Atlas

The above setup in WAsP are matched in FLORIS by:

1. Elevation -I don’t think we can explicitly specify orographic and landcover information in FLORIS. To ensure the same conditions as simulated in WAsP, a wind rose is derived from sectoral A and K values (Weibull distribution), from the GWC file, at 100 m elevation and 0.1 m roughness. This wind rose is saved in a csv file that matches the FLORIS template requirement (columns – ws, wd and freq_val) and is input using the parameter wind_rose. The wind rose is specified with wind speed bins of 0.25 and starting at 0.125 m/s and extending to 30.125 m/s.
2. Roughness - Same as 1
3. Resource grid – Same as 1
4. Single turbine – Turbine model specifications are defined in a yaml file. WTG specifications used are the same as WAsP. In addition, wind speed array and power array are also specified in the yaml file. Since the wind rose is specified with wind speed bins of 0.25 m/s, the power production curve is also specified using the same wind speed intervals. Since it’s a single turbine simulation and there is no wake loss consideration, the thrust_coeffcient array values are set to 0. The turbine model is called by specifying the turbine name in the gch yaml file using the parameter turbine_type.
5. Reference air density is set using the ref_air_density parameter in the turbine model yaml file and is set to 1.225 kg/m3.
6. Air density at the site – This is specified using the parameter air_density in the gch yaml file. Value used in 1.225 kg/m3.
7. The wind climate again is specified using the wind rose file derived from the A and K values specified in the GWC, at elevation 100 m and roughness 0.1 m. The wind rose.
   WAsP yields an AEP of 8.169 GWh.

I ran the calculation in a spreadsheet by matching the wind speed bins with power production curve and multiplying with frequency. I used 0.25 m/s granularity. This yielded 8.171 GWh, closely matching WAsP. The power density in the spreadsheet (153.75 W/m2) also closely matches the value in WAsP (154 W/m2). I have attached this spreadsheet with the message.

When I run the same calculation in FLORIS, I get 8.57 GWh as the AEP. I have attached the FLORIS files as well
Could you do a quick review and advice as to why there is such a big discrepancy between FLORIS and WAsP.
Compute AEP_Test.xlsx
FLORIS Files.zip

[misi9170]

Hi @niajhsaka ,

Thanks for clearly articulating your question and providing the source code needed to generate the issue---that's much appreciated!

Sorry that it has taken me some time to get back to you. I seem to have found why there is a discrepancy. The issue is that the WindRose.read_csv_long() method doesn't exactly respect the inputs provided---rather, the wind direction and wind speed bins are resampled. They use the same spacing (0.25m/s for wind speed), but instead of the bins being specified at 0.125, 0.375, 0.625, ... they are now at 0.0, 0.25, 0.50, ... . The change is not major, but does result in a change to the AEP as you noticed.

You can get around this by specifying the wind directions and wind speeds for the wind rose manually, i.e. using code such as

    wind_rose = WindRose(
        wind_directions=wds,
        wind_speeds=wss,
        freq_table=frequencies,
        ti_table=0.06
    )

Here is a modified version of the code you sent.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from floris import (
    FlorisModel,
    WindRose,
)

rebuild_wind_rose = True


fmodel = FlorisModel("inputs/gch_env_test.yaml")

# Load the wind rose from csv as in example 003
wind_rose = WindRose.read_csv_long(
    "inputs/wind_rose_V7.csv", wd_col="wd", ws_col="ws", freq_col="freq_val", ti_col_or_value=0.0
)
print(wind_rose.wind_speeds) # Don't exactly match the wind speeds in the csv file

df = pd.read_csv("inputs/wind_rose_V7.csv")
wds = df["wd"].unique()
wss = df["ws"].unique()


if rebuild_wind_rose: # Manually create wind rose
    wind_rose = WindRose(
        wind_directions=wds,
        wind_speeds=wss,
        freq_table=df["freq_val"].values.reshape(len(df["wd"].unique()), len(df["ws"].unique())),
        ti_table=0.06
    )
    print(wind_rose.wind_speeds) # Exactly match the wind speeds in the csv file

# Store some values
n_wd = len(wind_rose.wind_directions)
n_ws = len(wind_rose.wind_speeds)

# Store the number of elements of the freq_table which are 0
n_zeros = np.sum(wind_rose.freq_table == 0)

# Set the wind rose
fmodel.set(wind_data=wind_rose, wind_shear=0.0)

# Run the model
fmodel.run()

# Note that the frequency table contains 0 frequency for some wind directions and wind speeds
# and we've not selected to compute 0 frequency bins, therefore the n_findex will be less than
# the total number of wind directions and wind speed combinations
print(f"Total number of wind direction and wind speed combination: {n_wd * n_ws}")
print(f"Number of 0 frequency bins: {n_zeros}")
print(f"n_findex: {fmodel.n_findex}")

# Get the AEP
aep = fmodel.get_farm_AEP()

# Print the AEP
print(f"AEP from wind rose: {aep/1E9:.3f} (GWh)")

# Run the model again, without wakes, and use the result to compute the wake losses
fmodel.run_no_wake()

# Get the AEP without wake
aep_no_wake = fmodel.get_farm_AEP()

# Compute the wake losses
wake_losses = 100 * (aep_no_wake - aep) / aep_no_wake

# Print the wake losses
print(f"Wake losses: {wake_losses:.2f}%")

tpowers = fmodel.get_turbine_powers() / 1000.0  # Convert to kW

fig, ax = plt.subplots(1,1)
ax.scatter(wind_rose.wind_speeds, tpowers[0,:], marker="o", color="C0", label="simulated")
turb_dict_powercurve = fmodel.core.farm.turbine_definitions[0]["power_thrust_table"]
ax.scatter(turb_dict_powercurve["wind_speed"], turb_dict_powercurve["power"], marker="x",
           color="red", label="power curve")
ax.legend()
ax.set_xlabel("Wind Speed (m/s)")
ax.set_ylabel("Power (kW)")

plt.show()

If you set rebuild_wind_rose = False at the top of the script, the resulting AEP produced is AEP from wind rose: 8.573 (GWh) as you noted, and the following figure is now produced:

In particular, notice that there is a wind speed just above cut-out where there is an interpolated power value. This will probably contribute significantly to the AEP, even though it is not explicitly a point on the power curve.

On the other hand, you can try setting rebuild_wind_rose = True. Now, the wind speeds used in the wind rose exactly match those in the power curve, and the resulting AEP produced is AEP from wind rose: 8.167 (GWh) and the figure shows that the sampled wind speeds match the specified points:

Best,
Misha

[niajhsaka]

Thanks for the useful reply, this works!

[misi9170]

Great! I do have one follow-on suggestion---when using turbine power curves in FLORIS, you may want to specify the cut-out speed with very close resolution; that is, specify the (rated) power value at the cut out wind speed, and then specify 0 power at the rated wind speed plus a very small delta. This will help to avoid the possibility shown in the first figure above that the power value is interpolated to something between the rated power and zero if the wind speed is a bit above rated.

You'll find that we do this for the prepackaged turbines, e.g. here for the NREL 5MW.

You'll also need to add the extra term to your thrust coefficient table to match.

We also do the same thing for the cut-in wind speed, but there it is less critical because the power contribute in that band is very low anyway.