Question regarding turbine power calculations

[MarcusBinder]

Hello

I am trying to implement axial induction control into floris, but I have some questions for the power() function in the turbine.py file.

Are there any references, for why the power is calculated the way it is?
I can see that the comments mentions the paper from Ervin Bossanyi “Optimising yaw control at wind farm level”, but it looks like that has more to do with the todo comment, then the rest of the formulas.
The formula might be common knowledge for people more experienced with wind farm calculations, but I have only used the P = 1/2 rho A Cp V^3.

Also this might be more python technical, but how does this line function: power_interp[turb_type](yaw_effective_velocity)
I am not the most well versed in python, so I was curious how this ends up calculating the power output.

[Bartdoekemeijer]

Hi Marcus,

It is true that the standard equation from actuator disk theory for the turbine power is

P = 0.50 * rho * A *  V^3 * Cp(V)

with P the turbine power production, rho the air density, A the rotor swept area, V being the inflow wind speed, and Cp the power coefficient including any losses being a function of V. Actuator disk theory is only valid when the rotor is aligned with the inflow wind direction. Now, when the rotor is misaligned with the inflow and we have a nonzero yaw angle, we follow the cosine rule as widely described in the literature:

P = 0.50 * rho * A * V^3 * Cp(V) * cos(yaw)^Pp

where Pp typically has a value between 1.4 and 3.0.

Now this equation works mostly fine in regular region 2 operation. However, consider the situation where we have a NREL 5MW turbine operating at a very high wind speed, like 20 m/s. It will produce 5MW of power in that situation. If we then yaw by 20 degrees, assuming Pp=2.0, the equation above would predict a turbine power production of 4.415 MW. This is wrong, since the effective wind speed is still high enough for the turbine to be producing 5MW of power, even with the misalignment. Hence, instead, what we do is introduce an "effective wind speed". Thus, we include the cos(yaw)^Pp rule into the wind speed V and then use the nominal equation as:

pW = pP / 3.0
V_effective = V * cos(yaw)^pW
P = 0.50 * rho * A *  V_effective^3 * Cp(V_effective)

This ensures proper behavior in region 2.5 and region 3 operation of the turbine. It also is theoretically more correct since the turbine always operates on its power curve.

Finally, for Python, square brackets are to access a specific entry in a list, array or dictionary. In this case, power_interp is a Python dictionary containing a power interpolant function for every turbine type defined in the farm. Then, turb_types is simply a list (or numpy array) containing all the turbine types that we have in the wind farm. In most situations where we just have a single turbine type, we have:

power_interp = {'nrel_5MW': <scipy.interpolate.i...7dcb08f40>}
turb_types = ['nrel_5MW']

power_interp[turb_type] means that we are accessing the entry turb_type from the power_interp dictionary. Thus power_interp[turb_type] accesses power_interp['NREL_5MW'] and will return us the <scipy.interpolate> power interpolation function for the NREL 5MW wind turbine. Then, the round brackets indicate that power_interp[turb_type] is a function and we are providing it a single input, yaw_effective_velocity, which is what I called V_effective above. So power_interp[turb_type] contains the function y(input) = 0.50 * rho * A * input^3 * Cp(input), and we provide it input=yaw_effective_velocity.

For more general usage, I would suggest reading up on how to debug Python code in your editor. This is a hugely useful tool to understand the variables and the code structure. I hope this helps!

Best,
Bart

[MarcusBinder]

Hello Bart

Thanks for the fast reply!
That makes a ton of sense!
Can I ask why you divide the pP by 3 to get the pW, and what is going on with the (rho/1.225)**(1/3)?
Also I can see that the final power output is being multiplied with 1.225 before returning. Where does this factor come from?

Best regards
Marcus

[Bartdoekemeijer]

Hi Marcus,

Happy to help! The pW = pP/3 comes from pulling the cos(yaw) term into V in the second equation in the above post:

0.50 * rho * A * V^3 * Cp(V) * cos(yaw)^pP = 0.50 * rho * A * (V * (cos(yaw)^pP)^(1/3)) ^3  * Cp(V) 

where

V * (cos(yaw)^pP)^(1/3) = V * cos(yaw)^(pP/3) = V * cos(yaw)^(pW) = V_effective

The 1.225 term is best explained here: Issue #211. Let me know if you have additional questions after reading that post!

Cheers,
Bart

[MarcusBinder]

Hey Bart

Thanks for the explanation. Reading the other post makes perfect sense.
Again, thanks for taking the time to help me out!

Best
Marcus