Theory error in Gaussian wake models

[RyanDunn729]

Gaussian wake models have 2 errors that do not align with theory

This doesn't fit the typical issue template since it is just comparing code with theory. I will do my best to justify everything with references.

The first reference is the research paper that the Gaussian wake model is based off of:

Bastankhah, Majid, and Fernando Porté-Agel. 
"Experimental and theoretical study of wind turbine wakes in yawed conditions." 
Journal of Fluid Mechanics 806 (2016): 506-541.

which may be accessed publicly here. The main equations I am looking at are Eq. (6.4)
$$\frac{u_R}{u_\infty} = \frac{C_T \cos \gamma}{2 ( 1 - \sqrt{1-C_T \cos \gamma} )}$$
and Eq. (6.16)
$$\frac{x_0}{D} = \frac{\cos \gamma (1+\sqrt{1-C_T})}{\sqrt{2} [ 4 \alpha I + 2 \beta (1-\sqrt{1-C_T})]}$$

I am specifically referencing the code within the following files:

floris/simulation/wake_velocity/gauss.py

(herein referred to as wake_velocity) and

floris/simulation/wake_deflection/gauss.py

(herein referred to as wake_deflection)

Error 1 in wake_velocity

The first error occurs in wake_velocity, with respect to Eq. (6.4) where $u_R$ is calculated on line 87.

        uR = u_initial * ct_i / ( 2.0 * (1 - np.sqrt(1 - ct_i) ) )

This equation does not match up with Eq. (6.4), and it should be

        uR = u_initial * ct_i * cosd(yaw_angle) / ( 2.0 * (1 - np.sqrt(1 - ct_i) ) )

(note that cosine is multiplied in the numerator)

For additional verification, I also noticed that in wake_deflection lines 150-156, $u_R$ is calculated by

        uR = (
            freestream_velocity
          * ct_i
          * cosd(tilt)
          * cosd(yaw_i)
          / (2.0 * (1 - np.sqrt(1 - (ct_i * cosd(tilt) * cosd(yaw_i)))))
        )

further confirming the inconsistency with theory.

Error 2 in wake_deflection

The first error occurs in wake_deflection, with respect to Eq. (6.16) where $x_0$ is calculated on lines 160-166.

        x0 = (
            rotor_diameter_i
            * (cosd(yaw_i) * (1 + np.sqrt(1 - ct_i * cosd(yaw_i))))
            / (np.sqrt(2) * (
                4 * self.alpha * turbulence_intensity_i + 2 * self.beta * (1 - np.sqrt(1 - ct_i))
            )) + x_i
        )

This equation does not match up with Eq. (6.16), and it should be

        x0 = (
            rotor_diameter_i
            * (cosd(yaw_i) * (1 + np.sqrt(1 - ct_i)))
            / (np.sqrt(2) * (
                4 * self.alpha * turbulence_intensity_i + 2 * self.beta * (1 - np.sqrt(1 - ct_i))
            )) + x_i
        )

(note that cosine is removed within the square root)

For additional verification, I also noticed that in wake_velocity line 102, $x_0$ is calculated by

        x0 *= rotor_diameter_i * cosd(yaw_angle) * (1 + np.sqrt(1 - ct_i) )

further confirming the inconsistency with theory.

Floris version

v3.4 (latest)
I installed by

git clone git@github.com:NREL/floris.git

Before submitting I made sure to git pull and it appears this has not yet been spotted.

[xxxiiinn]

There is another difference that I found between the code and the theory.

The reference paper is also the :

Bastankhah, Majid, and Fernando Porté-Agel. 
"Experimental and theoretical study of wind turbine wakes in yawed conditions." 
Journal of Fluid Mechanics 806 (2016): 506-541.

In order to determine the potential core area at x=0, we should refer to the content in page 530.

And according to the formulas (6.4) and (6.7) given in reference paper:

$$\frac{u_R}{u_{\infty}} = \frac{C_T\cos{\gamma}}{2(1-\sqrt{1-C_T\cos{\gamma}})}$$

$$\frac{u_0}{u_{\infty}} = \sqrt{1-C_T}$$

The expression of $d\sqrt{u_R/u_0}$ should be:

$$d\sqrt{\frac{u_R}{u_0}} = d\sqrt{\frac{C_T\cos{\gamma}}{2(1-\sqrt{1-C_T\cos{\gamma}})(\sqrt{1-C_T})}}$$

But in current file, from line 131 to line 145

floris/simulation/wake_velocity/gauss.py

While determining the wake expansion in near wake zone, the linear ramp from 0 to 1, from the start of the near wake to the start of the far wake is calculated as:

near_wake_ramp_up = (x - xR) / (x0 - xR)
near_wake_ramp_down = (x0 - x) / (x0 - xR)

sigma_y = near_wake_ramp_down * 0.501 * rotor_diameter_i * np.sqrt(ct_i / 2.0)
sigma_y += near_wake_ramp_up * sigma_y0
sigma_y *= np.array(x >= xR)
sigma_y += np.ones_like(sigma_y) * np.array(x < xR) * 0.5 * rotor_diameter_i

sigma_z = near_wake_ramp_down * 0.501 * rotor_diameter_i * np.sqrt(ct_i / 2.0)
sigma_z += near_wake_ramp_up * sigma_z0
sigma_z *= np.array(x >= xR)
sigma_z += np.ones_like(sigma_z) * np.array(x < xR) * 0.5 * rotor_diameter_i

I place here a figure (0 yaw angle case) to make the above formulas easier to understand.

After doing a little geometry, the $\sigma_y$ in near wake zone can be calculated as:

$$\sigma_y = \text{near wake ramp up}*(\sigma_{y_0}) + \text{near wake ramp down} * (\frac{d}{2})\sqrt{\frac{u_R}{u_0}}$$

So the code should be like:

near_wake_ramp_up = (x - xR) / (x0 - xR)
near_wake_ramp_down = (x0 - x) / (x0 - xR)

sigma_y = near_wake_ramp_down * 0.501 * rotor_diameter_i * cosd(yaw_angle_i) * np.sqrt(ct_i * cosd(yaw_angle_i)/2.0/(1 - np.sqrt(1 - ct_i * cosd(yaw_angle_i)))/(np.sqrt(1 - ct_i)))
sigma_y += near_wake_ramp_up * sigma_y0
sigma_y *= np.array(x >= xR)
sigma_y += np.ones_like(sigma_y) * np.array(x < xR) * 0.5 * rotor_diameter_i

sigma_z = near_wake_ramp_down * 0.501 * rotor_diameter_i * np.sqrt(ct_i * cosd(yaw_angle_i)/2.0/(1 - np.sqrt(1 - ct_i * cosd(yaw_angle_i)))/(np.sqrt(1 - ct_i)))
sigma_z += near_wake_ramp_up * sigma_z0
sigma_z *= np.array(x >= xR)
sigma_z += np.ones_like(sigma_z) * np.array(x < xR) * 0.5 * rotor_diameter_i

Furthermore, @RyanDunn729, in wake_velocity line 102, I think the formula for $x_0$ is correct, because it is not finished, $x_0$ is calculated via line 102 and line 103.

x0 *= rotor_diameter_i * cosd(yaw_angle) * (1 + np.sqrt(1 - ct_i) )
x0 /= np.sqrt(2) * (
            4 * self.alpha * turbulence_intensity_i + 2 * self.beta * (1 - np.sqrt(1 - ct_i) )
        )