Velocity profiles in Floris

[fwasilczuk]

Hi,
I stumbled upon an incosistent code behaviour when calculating the velocity profiles.
I'm trying to investigate the impact of rotor averaging on the accuracy of simulations.
When changing the number of rotor points I can see the averaged rotor valocity changing, which results in the change of power and Ct (which is expected).
But when I'm comparing the velocity profiles there is no difference between cases with different number of points (which seems not consistent with the previous point).
I tried to investigate, by dumping the Ct from axial_induction function in /simulations/turbine library. My code (below) does "calculate_wake", where the axial induction function is called once, returning different values for different number of averaging points.
When "calculate_cross_plane" is used later, the axial induction is called three times. First two calls return the same Ct values for cases with different number of rotor averaging points. Only the last call returns the same values as for "calculate_wake", which are different for different number of rotor points. Maybe this has to do something with the fact that the velocity profiles don't differ from each other?
Can someone explain this code behaviour? Or maybe I'm doing something wrong?
Regards
Filip

Minimal working example:

import matplotlib.pyplot as plt
from floris.tools import FlorisInterface

def traverse(input_file="setup_file.yml", grid_points=1):
    fi = FlorisInterface(input_file)

    solver_settings = {
      "type": "turbine_grid",
      "turbine_grid_points": grid_points
    }

    fi.reinitialize(solver_settings=solver_settings)
    fi.calculate_wake()
    cross_plane = fi.calculate_cross_plane(y_resolution=3, z_resolution=41, z_bounds = (0, 400), downstream_dist=400)
    plane_values = cross_plane.df
    traverse = plane_values[(plane_values['x1'] == 0)][["x2", "u", "v", "w"]]
    power_turbine1=fi.get_turbine_powers()[0,0,0]/1000
    ct1=(fi.get_turbine_Cts())[0,0,0]
    WS_eff1=fi.get_turbine_average_velocities()[0,0,0]   

    return WS_eff1, power_turbine1, ct1, traverse 

point1=traverse(grid_points=1)
points49=traverse(grid_points=49)

print("Rotor average velocity: 1 point - %s m/s, 49 points - %s m/s" % (point1[0], points49[0]))
print("Power: 1 point - %s kW, 49 points - %s kW" % (point1[1], points49[1]))
print("Ct: 1 point - %s, 49 points - %s" % (point1[2], points49[2]))

plt.figure()
plt.plot(point1[3]['u'], point1[3]['x2'], label='1 point') 
plt.plot(points49[3]['u'], point1[3]['x2'], label='49 points', linestyle="--") 
plt.legend()
plt.show()

Setup_file:

name: GCH
description: Three turbines using Gauss Curl Hybrid model
floris_version: v3.0.0

logging:
  console:
    enable: true
    level: WARNING
  file:
    enable: false
    level: WARNING

solver:
  type: turbine_grid
  turbine_grid_points: 3

farm:
  layout_x:
  - 0.0
  layout_y:
  - 0.0
  turbine_type:
  - iea_10MW

flow_field:
  air_density: 1.225
  reference_wind_height: -1
  turbulence_intensity: 0.06
  wind_directions:
  - 270.0
  wind_shear: 0.12
  wind_speeds:
  - 10
  wind_veer: 0.0

wake:
  model_strings:
    combination_model: sosfs
    deflection_model: jimenez
    turbulence_model: none
    velocity_model: jensen

  enable_secondary_steering: true
  enable_yaw_added_recovery: true
  enable_transverse_velocities: true

  wake_deflection_parameters:
    gauss:
      ad: 0.0
      alpha: 0.58
      bd: 0.0
      beta: 0.077
      dm: 1.0
      ka: 0.38
      kb: 0.004
    jimenez:
      ad: 0.0
      bd: 0.0
      kd: 0.1

  wake_velocity_parameters:
    cc:
      a_s: 0.179367259
      b_s: 0.0118889215
      c_s1: 0.0563691592
      c_s2: 0.13290157
      a_f: 3.11
      b_f: -0.68
      c_f: 2.41
      alpha_mod: 1.0
    gauss:
      alpha: 0.58
      beta: 0.077
      ka: 0.38
      kb: 0.004
    jensen:
      we: 0.1

  wake_turbulence_parameters:
    crespo_hernandez:
      initial: 0.1
      constant: 0.5
      ai: 0.8
      downstream: -0.32

[rafmudaf]

Hi @fwasilczuk My apologies, I must have missed this message entirely. I just ran your script, and I indeed I see that the rotor average velocity, power, and Ct are different in expected ways. Since the 1-point case doesn't resolve the shear layer, the effective wind speed should be the same value as entered in the input file, 10 m/s. The 49-point case more closely captures the shear layer and it is skewed lower close to the ground, so the average velocity is slightly less than the input. This is reflected in the power and thrust, as well.

Rotor average velocity: 1 point - 10.0 m/s, 49 points - 9.97589992665184 m/s
Power: 1 point - 8114.221106820878 kW, 49 points - 8053.491442550159 kW
Ct: 1 point - 0.7747068376068377, 49 points - 0.7753453865587975

I also see that the two series on your plot are exactly the same. You're plotting a "cross plane" which is a cross section of the wake, so it's in the y-z plane with three points in the y-direction and 41 points in the z-direction. The points in z are distributed between the bounds (0, 400). "x2" in that case is z or vertical, so the plot is the axial component of velocity as a function of height above ground.

Note that your plotting function is using the same number of points in the z-direction for both cases: z_resolution=41, and that's probably intentional so that you can compare the two profiles.

    cross_plane = fi.calculate_cross_plane(y_resolution=3, z_resolution=41, z_bounds = (0, 400), downstream_dist=400)

I think you're expecting to see the difference in Ct above reflected in the calculated wake in the figure. This isn't happening because the FlorisInterface.calculate_NN_plane functions do some things behind the scenes that aren't quite obvious. First, they export the current Floris configuration to a dictionary to use later because the turbine-grid type solver doesn't apply to the visualization. Then, they create a different kind of grid-type, a full flow field grid, based on the resolution and bounds given as arguments. Next, the full flow field solver is run, but it turns out this does need a turbine-grid in order to get the turbine Ct's for the wake models. Instead of using the one saved earlier, it creates a new configuration with 3 grid points (3x3). It completes the calculation, saves the full flow field grid results, and attempts to reset the state of the FlorisInterface.Floris object back to what it was at the beginning of the function with the previously saved dictionary.

All this to say that if you want to do this comparison, it's probably easiest to modify the solver directly and export the data to a Numpy file to use in your plot later. Alternatively, you could connect the previously saved Floris settings through to the full_flow_solver function so that it behaves how you expect it to.

Let me know if you have any more questions around this.

[rafmudaf]

@fwasilcz I should have also mentioned that @vallbog recently added new methods to sample and plot velocity profiles in #699 through FlorisInterface.sample_velocity_deficit_profiles. This method uses the same solver algorithm as FlorisInterface.calculate_cross_plane so the turbine grid is still hard-coded. However, the methods for doing the calculation and plotting the profiles in introduced in #699 are very useful.

[fwasilcz]

Hi @rafmudaf,
Thanks for the explanation, I expected the reason is something like you've explained, but it's nice to know the details :)
Regards
Filip

[rafmudaf]

Great! I'll mark the response above as the answer for the benefit of anyone coming to this discussion in the future.