Load Models

Stokes’ theory of ocean waves assumes irrotational flow of an incompressible fluid, described by a scalar velocity potential \(\mathrm {\Phi(\vec{x},t)}\) which satisfies

with a velocity field \(\mathrm {\vec{u}}\) given by

Assuming an impermeable sea-floor of constant mean depth \(\mathrm {h}\), we have \(\mathrm {u_z(z=-h)=0}\), which we write in terms of the wave potential as

The free surface of the wave is given by \(\mathrm {z=\eta (x,y,t)}\). For sea states with long-crested waves the \(\mathrm {y}\) coordinate is irrelevant, while for short-crested waves it is sometimes useful to transform between a cartesian and a polar coordinate system, such that \(\mathrm {z=\eta (r,\theta ,t)}\). Now, \(\mathrm {r}\) plays the role which \(\mathrm {x}\) plays in long-crested waves, and \(\mathrm {\eta }\) results from the superposition of several \(\eta _{\mathrm i}(r,t)\equiv\eta (r,\theta _{\mathrm i},t)\).

Two boundary conditions are imposed at the wave surface. The kinematic boundary condition

which states that a particle at the surface moves together with the surface, and the dynamic boundary condition

which follows from assuming constant pressure along the surface. Here, \(\mathrm {g}\) is the acceleration of gravity.

Together with Equation (26), which applies within the bulk of the ocean, Equation (28), specify the set of differential equations from which \(\mathrm {\Phi}\) is obtained.

In order to obtain analytic results for \(\mathrm {\Phi}\), a perturbative approach is taken.

The velocity potential at the wave surface is extrapolated onto the mean surface level \(\mathrm {z=0}\) using a Taylor expansion, with terms of higher order in the expansion expected to contribute at higher order in the ordering parameter \(\mathrm {\epsilon}\).

The expansion is, to order \(\mathrm {\mathcal{O}(\epsilon^2)}\),

Likewise, the surface elevation is written as a series expansion,

where it is usually assumed that \(\mathrm {\eta ^{(1)}}\) is a linear wave. This is a construction which allows us to study the non-linear effects which appear at second order and beyond.

The surface elevation is a function of \(\mathrm {x}\), \(\mathrm {y}\) and \(\mathrm {t}\), or equivalently and for many purposes more conveniently, a function of \(\mathrm {r}\), \(\mathrm {\theta }\) and \(\mathrm {t}\). Staying with the former set of coordinates, the surface elevation in Fourier space is a function of wavenumber \(\mathrm {\vec{k}=(k_x,k_y)}\) and angular frequency \(\mathrm {\omega }\). They are related through the dispersion relation for surface waves at finite water-depth,

which is complete to second order. At third order, the dispersion relation includes an additional term.

The first order contribution \(\mathrm {\eta ^{(1)}}\) is provided to WaveKin by its Fourier components \(a_{\mathrm i}\),

where the sums run over the number of positive frequencies \(\mathrm {N}\) and the number of directions \(\mathrm {M}\) in the discretized directional spectrum, \(\psi_{\mathrm {ik}}=\omega _{\mathrm i}t-\vec{k_{\mathrm {ik}}}\cdot \vec{x}-\varepsilon _{\mathrm {ik}}\) is the phase function, with angular frequency $_{i}={i}$. The phase is \(\varepsilon _{\mathrm {ik}}\) and the wavenumber \(\vec{k_{\mathrm {ik}}}\) is obtained from \(\omega _{\mathrm i}\) through use of the dispersion relation and wave direction \(\theta _{\mathrm {ik}}\). The two subscripts denote a discretization of the Fourier components of a directional spectrum. In the following, the second subscript \(\mathrm {k}\) is dropped in favor of a compact notation. It should be implicitly understood that for directional spectra, all sums over Fourier components also imply a sum over directions.

For a simulated wave from a given spectrum, the random phase will typically be drawn from a uniform probability distribution. For a directional spectrum this involves \(\mathrm {MN}\) independent and identically distributed random numbers. See Section 3.4.2 for a further discussion on this.

The sum of Equation (3) can be re-written to run over positive frequencies only. Because \(\mathrm {\eta (x,y,t)}\) is real-valued, \(a_{-\mathrm i}=\bar{a}_{\mathrm i}\). This gives

where the sum over directions is now implicit.

Velocities and accelerations are also real-valued; WaveKin therefore works with one-sided spectra only. Time-series may have a non-zero mean value, therefore the one-sided spectra include the zero-frequency components. The zero-frequency component for surface elevation \(\mathrm {a_0}\) is also included in the one-sided spectrum, though its value should be zero or near zero. See Section 3.4.1.6 for how symmetries allow summations of first and second order terms to run over positive frequencies only and Section 3.4.1.7 for a further discussion on the zero-frequency components.

The expansion is in terms of an ordering parameter \(\mathrm {\epsilon}\). It can be seen that higher orders in the expansions of the boundary conditions at the wave surface results in terms with higher powers of spatial derivatives, surface elevation, and the velocity potential. In Fourier space, spatial derivatives result in additional factors of \(\mathrm {k}\), and we have constructed the first order surface elevation to be a sum of harmonics with amplitude \(\mathrm {|a|}\).

The combination of factors of \(\mathrm {k}\) and \(\mathrm {a}\) provides a physical meaning to the ordering parameter. The wave steepness is proportional to \(\mathrm {ak}\), for a harmonic wave elevation \(\mathrm {\eta _k=a\sin kx}\). In order for the series expansion to converge, the linear wave should have a spectral content which satisfies \(|a_{\mathrm i}|k_{\mathrm i}<\epsilon\) for all \(\mathrm {i}\).

In addition to the Taylor series expansion of the wave surface boundary conditions and the expansion of \(\mathrm {\eta }\), the velocity potential is written as a series expansion

A solution is obtained for each term of the expansion by ordering the terms in the wave surface boundary conditions, setting \(\mathrm {z=\eta }\) in Equation (3). Terms at equal order in \(\mathrm {\epsilon}\) are solved for independently.

Solving for the first order velocity potential gives

Following Marthinsen and Winterstein (1992) for long-crested waves and Marthinsen and Winterstein (1992) for short-crested waves, the velocity potential and second order contribution to surface elevation is given to second order by a set of rather complex expressions, organized in terms of sum-frequency and difference-frequency contributions. The expressions are analytic, and therefore support differentiation, allowing the analytic calculation of velocities and accelerations.

Solving for the second order surface elevation gives

Solving for the second order velocity potential gives

where the latter term is independent of spatial coordinates and therefore irrelevant for kinematics.

The second order kernels \(E_{\mathrm {i,j}}\) and \(P_{\mathrm {i,j}}\) are given by

and

The solutions for both \(\mathrm {\eta ^{(2)}}\) and \(\mathrm {\Phi^{(2)}}\) depend on products, sums and differences of Fourier components. Because the two-sided spectrum can be expressed as a one-sided spectrum, the second order contributions are naturally separated into sum-frequency and difference-frequency contributions.

Because the sum and difference-frequency components only depend on the sum and difference of frequencies (wavenumbers) and the Fourier components satisfy \(a_{-\mathrm i}=\bar{a}_{\mathrm i}\), the calculations of second order contributions to surface elevation, velocities and accelerations can be considerably simplified.

The kernels from Equation (41) and Equation (42) can be written in terms of sums and differences of frequencies and corresponding wavenumbers,

and

where \(\mathrm {i,j}>0\).

The symmetry properties of the second order kernels are

where \(\mathrm {i,j}>0\).

For the second order wave elevation of Equation (3) we now have \(\mathrm {\eta ^{(2)}=\eta ^{}\eta ^{-}}\) with

The second order velocity potential of Equation (40) can now be written as \(\mathrm {\Phi^{(2)}=\Phi^{}\Phi^{-}}\) with

where space-independent terms have been omitted.

The difference-frequency kernels \(E_{\mathrm {ij}}^-\) and \(P_{\mathrm {ij}}^-\) are singular for \(\mathrm {i}=\mathrm {j}\). If one chooses to take the analytic limit, one gets a non-zero contribution to the zero-frequency contribution of the difference-frequency transfer functions. For \(E_{\mathrm {ii}}^-\) this gives

where

is the group velocity.

Equation (51) follows Marthinsen and Winterstein (1992).

With non-zero contributions from the difference-frequency kernels when \(\mathrm {i}=\mathrm {j}\), the resulting transfer functions give finite contributions to the zero-frequency component of the second order elevation and velocity potential. This corresponds to finite mean-values for \(\mathrm {\eta }\), \(\mathrm {u_x}\) and \(\mathrm {u_y}\). This can optionally be avoided by setting all zero-frequency components to zero, (equivalently, setting the diagonal entries of the difference-frequency kernels to zero) however this procedure lacks physical or theoretical justification.

In particular, the so-called mean sea state set-down leads to an inconsistency with the boundary condition applied at the sea floor, i.e. at \(\mathrm {z=-h}\), as the sea floor depth \(\mathrm {h}\) should properly be measured from the mean sea level, \(\mathrm {<\eta >}\), not from the first order contribution to the same, \(\mathrm {<\eta ^{(1)}>}\). The inconsistency might be small for most practical uses, but deserves mentioning. This is discussed occasionally in the literature, e.g., Marthinsen and Winterstein (1992) and The default behaviour of WaveKin in RIFLEX is to set the diagonal terms equal to zero.

Depending on how the target spectrum is filtered to obtain the first order wave elevation defined by \(\{a_{\mathrm {i}}\}\), the sum-frequency transfer functions may result in significant high frequency contributions to the second order wave elevation and kinematics. The acceleration terms are particularly sensitive to the high frequency part of the spectrum. Care must be taken when defining the first order wave elevation to avoid un-physical high frequency terms arising at second order. Otherwise, un-physically large accelerations could result.

The issue of large accelerations, i.e., rapidly varying velocities, is most acute near the mean sea surface and above (\(\mathrm {z\ge0}\)).

By target spectrum is meant the empirically established power spectrum \(\mathrm {S(\omega )}\) for long-crested waves or the directional spectrum \(\mathrm {S(\omega ,\theta )}\) for directional sea states.

In order to generate an irregular sea state, phases are typically drawn from a uniform distribution, with the amplitude of Fourier components fixed by a filtered power spectrum. The resulting wave elevation is then referred to as a linear wave and identified as the first order contribution in the perturbative expansion of \(\mathrm {\eta }\).

Various heuristics are presented in the literature for how the filtering should be done, see e.g. Marthinsen and Winterstein (1992) and Marthinsen and Winterstein (1992), but no solid theoretical foundation underlies these methods. For many purposes, results may be insensitive to the details of the filtering, however results from WaveKin indicate that filtering has a significant effect on wave kinematics near the sea surface at water depths of \(20-40{\mathrm {~m}}\) and significant wave heights in the range of \(5-6{\mathrm {~m}}\).

Equating the linear wave with the first order wave elevation is a matter of convenience, but could under-estimate the non-linear nature of Stokes’ waves. The assumption of un-correlated phases is justified by very weak wave-wave interactions, which in turn applies to low wave steepness. While the first order wave elevation will always have a more narrow banded spectrum than the full wave, it will nevertheless contain steep waves.

A directional spreading function \(\mathrm {D(\theta ,\omega )}\) is defined by

together with the normalization condition

The power spectrum \(\mathrm {S(\omega )}\) for directional sea states is obtained by integrating over the directional coordinate

The simplest approach to describing a directional sea state is to assume a separable directional spectrum,

where the power spectrum \(\mathrm {S(\omega )}\) (e.g., a JONSWAP spectrum) and the directional spreading function \(\mathrm {D(\theta )}\) are given separately and together define the sea state.

A common choice for \(\mathrm {D}\) is a frequency-independent cosine-2s distribution

for \(\mathrm {|\theta -\theta _0|\leq\frac{\pi }{2}}\) and zero otherwise. \(\mathrm {C_s}\) is a normalization factor. The width of the distribution is parametrized by \(\mathrm {s}\) and the mean direction is given by \(\mathrm {\theta _0}\).

Frequently, convective terms are neglected in the calculation of wave kinematics. In the current implementation of WaveKin the inclusion of these terms is optional. If they are omitted, only the contribution from partial time-derivatives are included

Convection terms can optionally be included in the accelerations returned by WaveKin. These are calculated by taking spatial derivatives of the velocity field (in the frequency domain) and multiplying with the velocity field (in the time domain).

Convection terms contribute to the accelerations of momentum-carrying fluid particles, and thus contribute to the inertia term of Morison’s load model. They emerge when we take the total derivative of velocity with respect to time,

where \(\mathrm {\vec{v}}\) is the velocity vector and \(\mathrm {\vec{u}=\vec{u}(\vec{x},t)}\) is the velocity field.

The inclusion of convection terms is turned off by default in WaveKin. This follows an assumption that they are small and can be neglected for practical purposes. This assumption should be questioned for severe sea states. The last term in Equation (61) constitutes the convection terms. It should be seen that its magnitude depends on the gradient of the velocity field, meaning that rapidly varying velocities, as found near the sea surface in severe sea states, may lead to significant contributions from convection terms.

The velocity potential \(\mathrm {\Phi(x,y,z,t)}\) is extrapolated above \(\mathrm {z=0}\) according to the Taylor expansion of Equation (3). Kinematics are obtained from the extrapolated potential by spatial and temporal differentiation in the same manner as for the velocity potential at and below \(\mathrm {z=0}\). However, it should be noted that the second order contribution at \(\mathrm {z>0}\) only involves the first order velocity potential \(\mathrm {\Phi^{(1)}}\). There are no quadratic transfer functions involved above mean surface level until third order contributions.

High frequencies resulting from providing an improperly filtered spectrum as the input first order wave may cause unphysical accelerations near to and above the mean surface level. To compensate for this, an option is provided to filter out large frequencies prior to extrapolating the kinematics above the mean surface level.

Together with, or instead of this option, there is an option to set a threshold surface level \(z_{\mathrm {lim}}<0\), above which kinematics are extrapolated from the potential at \(z=z_{\mathrm {lim}}\). As the large frequencies are rapidly dampened even at very small values of \(\mathrm {z<0}\), setting \(z_{\mathrm {lim}}\) to only a few centimeters below \(\mathrm {z=0}\) may be sufficient to avoid the excessive spikes in acceleration.

The default behaviour of WaveKin in RIFLEX is to use the entire provided spectrum when extrapolating kinematics above the mean surface level, and to extrapolate from \(\mathrm {z=0}\).

In order to extrapolate kinematics between two nodes, one of which is dry, it can be impractical to set kinematics to zero above the sea surface \(\mathrm {z>\eta }\) (sometimes referred to as a dry node). Instead, the option exists to set kinematics at \(\mathrm {z>\eta }\) equal to the the kinematics at the surface: \(\mathrm {z=\eta }\). This is the default behaviour of WaveKin in RIFLEX.

For kinematics on the surface and at or above the mean surface level, \(\mathrm {z=\eta \ge0}\), the kinematics are obtained at second order by extrapolation of the first order kinematics at \(\mathrm {z=0}\). For kinematics on the surface and below the mean surface level, \(\mathrm {z=\eta <0}\), there is a choice between using the second order potential, which exists at all \(\mathrm {z\le0}\) or using the first order potential at \(\mathrm {z=\eta <0}\) together with the first order potential extrapolated from \(\mathrm {z=0}\). The latter method provides a consistent calculation of kinematics on the surface regardless of whether it is above or below \(\mathrm {z=0}\). However, it causes a discontinuity in the kinematics on and directly beneath the surface.

The magnitude of the discrepancy between the two methods should be of third order, if the series expansion behaves consistently. Thus, a comparison of the two methods is a test of the series expansion.

BEM theory is not valid for induction factors greater than 0.4. The Glauert correction is used for large induction factors: for \(\mathrm {a>0.4}\), the empirical thrust curve recommended by Burton et al (2011, p67) is used:

where \(\mathrm {a_2=1.0}\), \(\mathrm {C_{T2}=1.82}\), \(\mathrm {a_1=1.0-0.5\sqrt{C_{T2}}}\), \(\mathrm {C_{T1}=4a_1(1-a_1)}\), and \(\mathrm {F}\) is the Prandtl factor Equation (141). For \(\mathrm {F=1}\), the thrust curve used in the calculations is compared to the momentum theory result shown in the figure Figure 17.

To account for the tip and hub loss on a blade due to a finite number of blades, the Prandtl factor is implemented. The Prandtl factor \(\mathrm {F}\) as a function of the radial location \(\mathrm {r}\) is computed as:

where

\(\mathrm {B}\) is the number of blades, \(\mathrm {R}\) is the outer radius of the rotor, and \(\mathrm {\phi }\) is the angle that the trailing vorticies make with the rotorplane. \(\mathrm {\phi }\) is assumed to be the same as the angle of incoming flow, neglecting tangential induced velocity. The default program behavior is to correct \(\mathrm {\phi }\) for yawed inflow.

The dynamic wake effect is the time lag in induced velocities due to the shedding and downstream convection of vorticity. Dynamic wake effects are most pronounced for heavily loaded rotors, corresponding to high induction factors (low wind speeds). This effect can be modeled by the Stig Øye dynamic inflow model, which acts as a filter for induced velocities, as shown below (Snel and Schepers, 1995). \(\mathrm {\boldsymbol{W_{qs}}}\) is the quasi static induced velocity vector, and \(\mathrm {\boldsymbol{W}}\) is the new induced velocity vector. \(\mathrm {\tau_1}\) and \(\mathrm {\tau_2}\) are time constants.

Dynamic stall refers to time-dependent variations in the lift and drag coefficients of an airfoil due to changes in the angle of attack. The Stig Øye model is implemented and gives unsteady lift by filtering the trailing edge separation point with an empirical time constant (see also Hansen, 2008).

The degree of stall (\(\mathrm {f_s}\)) affects the lift as

where \(C_{L,\mathrm {inv}}\) is the lift coefficient without separation, and \(\mathrm {C_{L,fs}}\) is the fully separated lift coefficient. A static value \(\mathrm {f_{s}^{st}}\) can be found to reproduce the static airfoil data. The degree of stall is assumed to chase the static value:

where \(\mathrm {\tau}\) is a time constant. By integration:

For flow that is not fully stalled, the lift coefficient is computed as:

where \(\mathrm {\frac{dC_L}{d\alpha }}\) is computed at the full-stall limit (depending on the sign of \(\mathrm {\alpha }\)), and \(\mathrm {\alpha _0}\) is the angle of attack where the lift is zero.

For fully stalled flow:

where \(\mathrm {C_{L,qs}}\) is the quasi-static lift.

In order to initialize the dynamic stall calculation, the program identifies several important points in the lift curve: the angle of attack where there is zero lift (\(\mathrm {\alpha _0}\)), the maximum slope of the lift curve in the linear region for both positive and negative lift values (\([dC_L/d\alpha ]_{\mathrm {max}}(1)\) and \([dC_L/d\alpha ]_{\mathrm {max}}(2)\)), and the angle of attack corresponding to full separation for both positive and negative lift values (\(\mathrm {\alpha _{fs}(1)}\) and (\(\mathrm {\alpha _{fs}(2)}\))) For typical airfoils, these points are relatively easy to identify. For more complex airfoils, the user may input these values directly.

Skewed inflow may occur due to rotor tilt or a yaw angle between the rotor and the oncoming wind. Glauert developed the basic formulation for correcting the induction factor due to skewed inflow: \(\mathrm {a}\) is the induction factor, \(\mathrm {\chi}\) is the wake skew angle, and \(\mathrm {\Psi}\) is the azimuth angle that is zero at the most downwind position of the rotor.

The tower influence (shadow) effect is implemented using potential flow theory. Since the incoming wind has to travel around the tower, the tower has an effect on the local inflow, even for an upwind turbine. The velocity deficit due to the tower is based on the 2D potential solution for constant flow around a circle, with the possibility to include corrections for the tower drag and the Bak coefficient. The coordinate system for the tower influence calculation is shown in figure Figure 19.

For an upwind wind turbine, the non-dimensional influence at a location \(\mathrm {(x,y,z)}\) is evaluated as:

where \(D_{\mathrm {tow}}\) is the local tower diameter for the given \(\mathrm {z}\) level, \(\mathrm {b}\) is the Bak coefficient, and \(\mathrm {C_D}\) is the local tower drag coefficient for the given \(\mathrm {z}\) level. If the effect of tower drag is not considered, the Bak factor is set to zero and potential flow solution is applied.

In practice, the horizontal wind direction and speed are computed at the location \(\mathrm {(x,y,z)}\). This speed is multiplied by the velocity factors \(x_{\mathrm {infl}}\) and \(y_{\mathrm {infl}}\), and the modified speed is then transposed back to the initial wind direction (Moriarty and Hansen, 2005).

Note that the \(\mathrm {x}\)-coordinate here is upwind and the velocity is set to zero for any point inside the tower. The tower is assumed to be approximately vertical and no effect on the vertical wind speed is included.

When the blade segment is above the tower height, the velocity is smoothened to match the free wind. The resulting velocity deficit factors is shown in the figures below. The blade is pointing upwards when the azimuth angle is 0 deg, and is closest to the tower at 180 deg. In these figures the distance from the hub to the center of the segment is 90 m, the tower top is located 6 m below the hub and the tower diameter is ranging from 9 m (azimuth=180 deg) to 6.54 m (tower top). There is no cone on the blade line and no tilt of the shaft.

In order to model downwind wind turbines, a different method for computing the velocity deficit due to the presence of the tower is needed. The most common models, which are implemented in AeroDyn and Bladed, for example, are based on work by Powles. Powles’ model, as implemented in RIFLEX, assumes a cosine-squared shape for the tower shadow.

The cosine-squared model for the wake factor \(u_{\mathrm {wake}}\) is: $ u_\{}=

$ where \(\mathrm {d=\sqrt{x^2+y^2}}\). Together with the potential flow model, the local wind velocity is found as: $ U_\{}=(u-u_\{})U_$ $ V_\{}=(v-u_\{})U_$

where \(\mathrm {U_\infty}\) is the total instantaneous horizontal velocity, and \(\mathrm {u=1}\) when the cosine-squared model is in effect. In Powles’ model, the flow is assumed to not reverse: for \(u-u_{\mathrm {wake}}<0\), we assume \(U_{\mathrm {local}}=0\). The computed \(U_{\mathrm {local}}\) and \(V_{\mathrm {local}}\) are then transformed back to the local global reference frame. Thus, the wake is assumed to align with the instantaneous horizontal wind vector.

The resulting velocity deficit for a downwind turbine is shown in figure Figure 22.