2.1 One-dimensional mean flow turbine analysis

The coordinate directions in the streamwise and (any normal pair of) transverse directions are denoted by $x_{1}$ , $x_{2}$ , $x_{3}$ respectively. The corresponding mean flow velocity components are $(\overline{U}_{1},\overline{U}_{2},\overline{U}_{3})$ . The turbulent velocity components are $(u_{1},u_{2},u_{3})$ , with $(\unicode[STIX]{x1D705}_{1},\unicode[STIX]{x1D705}_{2},\unicode[STIX]{x1D705}_{3})$ the three components of the wavenumber vector $\unicode[STIX]{x1D73F}$ , in the corresponding directions. The suffix $\infty$ is used to indicate conditions sufficiently far upstream to be independent of the rotor and the suffix $0$ to indicate values at the rotor plane, usually at the origin which is at the centre of the rotor disc (figure 1).

The simplest model for the mean flow field of a HAWT (or HATT) rotor is the actuator disc model proposed by Froude (Reference Froude1911), Betz (Reference Betz1920) and Joukowsky (Reference Joukowsky1920). This replaces the disc swept by the rotor blades by a thin sheet of resistance to the flow. The model is the same as that proposed by Taylor (Reference Taylor1944) for a porous plate and indeed porous discs are often used as experimental analogues of rotors for wind tunnel testing of arrays. Within limits the mean flow field is represented quite well by this model but it omits the mean swirl in the wake.

From an axial momentum balance Betz (Reference Betz1920) showed that the rate of conversion of flow energy into mechanical energy is directly related to the induction of backflow by the rotor wake, which slows the mean flow into the rotor so that $\overline{U}_{1}(x_{1}=0)=(1-a)\overline{U}_{1\infty }$ where $a$ is the wake induction factor. In this simple model the rotor wake lies within a cylindrical sheet of vorticity and the induction factor and hence the mean streamwise velocity at the rotor plane are both constant over the disc. The strength of the vortex sheet is $2a\overline{U}_{1\infty }$ representing the helical vorticity shed from the rotor blade tips but omitting the small axial component which generates swirl. Far downstream the axial velocity within the wake is therefore $\overline{U}_{1w}=(1-2a)\overline{U}_{1\infty }$ and outside the wake it is $\overline{U}_{1\infty }$ .

Axial momentum analysis shows that this one-dimensional flow model predicts the power coefficient to be

where $P$ is the power extracted from the flow and $A(=\unicode[STIX]{x03C0}R^{2})$ is the swept area of the rotor disc.

This result leads to the well-known Betz limit that the maximum value possible for the power coefficient of an isolated open rotor is $16/27$ attained when $a=1/3$ . The result is sufficiently accurate for most real HAWTs and HATTs to be operated as close to this condition as possible to maximise efficiency. The model was later developed for propellers by Glauert (Reference Glauert1935) into the BEM method which forms the basis of industrial HAWT and HATT rotor design to the present day.

The actuator disc model of the mean flow field will be used in the present paper for the purpose of assessing distortion of turbulence incident on a rotor as it is the simplest model which is sufficiently realistic. This is so for high tip-speed-ratio rotors even if they have only two or three blades as is now the usual configuration of HAWTs and HATTs.

Similarly it will be assumed for the purposes of assessing the effect of distortion of incident turbulence that any turbulent flow of interest in the atmosphere or in a tidal channel can be represented by homogeneous isotropic turbulence convected by a uniform stream, conditions being based on those at the hub height of the rotor in the turbulent shear layer of the ABL or tidal flow. In what follows the mean free stream velocity will be taken as $(\overline{U}_{1\infty },\overline{U}_{2\infty },\overline{U}_{3\infty })=(1,0,0)$ so that all velocities are in effect non-dimensionalised by the free stream mean velocity and the radius of the actuator disc will similarly be taken to be $R=1$ so that all length scales are non-dimensionalised by the rotor radius.

The effect of the turbine on the turbulence when the turbulence length scale is much smaller than the rotor diameter: $L_{1\infty }\ll 1$ is examined first.

Figure 1.

Mean flow streamlines through an actuator disc, induction factor $a=1/3$ .

Following the RDT analysis of BP for locally homogeneous strain the only information required is of the mean flow conditions at the position at which the distorted turbulent velocities are to be calculated relative to the upstream conditions. In the present work the distorted velocity field is to be calculated at the rotor plane in order to assess the effect of distortion on what would normally be assumed to be frozen incident turbulence impacting the rotor. Therefore the blocking effect of the blade surfaces should not be included in the distortion calculation. The issue of rotor mean wake distortion by the turbulent velocities is considered later but has no effect at the rotor plane in the limit of very small turbulence length scale. Since the buffet forces on rotor blades operating at high tip speed ratio (as is usual) are mainly sensitive to the axial $u_{1}$ component of the turbulence, therefore in this paper this velocity component will be calculated. If the disc centre ((0, 0, 0) in figure 1) is taken as a representative point the simplifications of axisymmetry and plane strain can be used which reduce the computational effort very considerably. But off the axis elements of the turbulence vorticity are subject to rotation (tilting) in addition to the stretching or compression by the pure strain $(\unicode[STIX]{x2202}\overline{U_{i}}/\unicode[STIX]{x2202}x_{i})~i=1,2$ or 3. Durbin (Reference Durbin1981) has derived this contribution for the small length scale case of turbulence incident on an axisymmetric body. His approach will be followed here for points off the axis.

2.2 Distortion of small length scale turbulence

The distorted streamwise turbulent velocity $u_{1}$ at the centre of the disc (0,0,0) is first considered. The mean flow axial velocity at (0,0,0) is

which corresponds to a mean flow expansion.

The actuator disc flow field close to the axis is the opposite of that for a wind tunnel contraction and is similarly a pure strain. The BP analysis as given in Batchelor (Reference Batchelor1953) develops the rapid distortion analysis for this case in terms of the mean flow contraction coefficient $c$ $(=\overline{U}_{1}(final)/\overline{U}_{1\infty }(initial))>1$ . For the actuator disc case the equivalent ‘contraction’ coefficient $c<1$ . Writing $\unicode[STIX]{x1D73B}(\boldsymbol{x})=\unicode[STIX]{x1D735}\times \boldsymbol{u}(\boldsymbol{x})$ as the vorticity of the turbulent fluctuations at a general point $\boldsymbol{x}$ with corresponding value far upstream $\unicode[STIX]{x1D73B}_{\infty }(\unicode[STIX]{x1D743}_{\infty })$ then the Lagrangian specification of the motion used by Batchelor shows that in the absence of viscous forces

where $\boldsymbol{x}$ is the coordinate of a fluid particle originally upstream with Lagrangian coordinate $\unicode[STIX]{x1D743}_{\infty }$ on the same streamline. Thus the vorticity is intensified or diminished and distributed between components according to the strain tensor $(\unicode[STIX]{x2202}x_{i}/\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{j\infty })$ .

The turbulent velocity field expressed in wavenumber form

where $Gk\unicode[STIX]{x1D714}$ is angular frequency, can be calculated from the vorticity field through the relationship

If the strain tensor may be considered to vary very little over the length scale of the turbulence, which is the strain homogeneity assumption for small length scale turbulence, then (2.6) may be solved locally with the strain tensor terms as constants rather than solving over the whole strain field. Therefore in wavenumber form the velocity at the disc can be expressed in terms of a distorted wavenumber $\unicode[STIX]{x1D74C}$ by

where the distorted wavenumber $\unicode[STIX]{x1D74C}$ is related to the undistorted wavenumbers $\unicode[STIX]{x1D73F}$ by

and $\unicode[STIX]{x1D750}_{ijm}$ is the unit alternating tensor. In the case of pure strain as for a wind tunnel contraction or for the rotor flow near the axis only the diagonal terms of the strain tensor, $(\unicode[STIX]{x2202}x_{1}/\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{1\infty })=c,(\unicode[STIX]{x2202}x_{2}/\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{2\infty })=(\unicode[STIX]{x2202}x_{3}/\unicode[STIX]{x2202}\unicode[STIX]{x1D709}_{3\infty })=1/\sqrt{c}$ , are non-zero. The relationship between the spectral tensors $\unicode[STIX]{x1D731}_{11}$ of the distorted and undistorted turbulence, time mean products of the wavenumber amplitudes $\hat{u} _{1}$ with their complex conjugates, then simplifies greatly to

The spectrum tensor for the undistorted upstream streamwise turbulence velocity, assumed isotropic and homogeneous, has the form (Batchelor Reference Batchelor1953)

where $\unicode[STIX]{x1D70F}^{2}=\unicode[STIX]{x1D705}_{2}^{2}+\unicode[STIX]{x1D705}_{3}^{2}$ and $E(\unicode[STIX]{x1D705})$ is the turbulence energy spectrum function. Therefore, its $(\unicode[STIX]{x1D705}_{1})$ frequency spectrum is

and the frequency spectrum for the distorted streamwise turbulence velocity at the disc in the region of the axis becomes

Integrating these two spectra with respect to $\unicode[STIX]{x1D705}_{1}$ , substituting $\unicode[STIX]{x1D705}_{1}=\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D705}\sin \unicode[STIX]{x1D703}$ , the ratio of the mean square of the streamwise turbulence velocity at the disc centre to its upstream value is

where $s=c^{-3}-1$ . This result, which may be obtained directly from the BP result setting $c<1$ and taking account of the imaginary square root, was given by Farr & Hancock (Reference Farr and Hancock2012). It is close to but just less than the upper limit of the amplification ratio of intensity due to distortion over all length scales. It is independent of the form of the turbulence energy spectrum function $E(\unicode[STIX]{x1D705})$ which divides out in taking the ratio. Figure 2 shows a plot of this ratio versus the induction factor $a=(1-c)$ . In this case, $L_{1\infty }\rightarrow 0$ , the distortion of the streamwise turbulence intensity $\overline{u_{1}^{2}}/\overline{u^{2}}_{1\infty }=1.341$ when $a=1/3$ . For $a>1/2$ actuator disc theory predicts mean flow reversal in the wake and although wake reversal does occur in reality, the model breaks down.

Figure 2.

Mean square turbulence intensity ratio $\overline{u_{1}^{2}}/\overline{u_{1\infty }^{2}}$ versus induction factor $a$ for $L_{1\infty }\rightarrow 0$ .

In order to calculate the spectrum $\unicode[STIX]{x1D719}_{11}$ of the distorted turbulence it is necessary to specify the form of the energy spectrum function $E(\unicode[STIX]{x1D705})$ . The turbulence of the wind in a neutral ABL is usually represented by an empirical spectrum such as the von Kármán spectrum (Deaves & Harris Reference Deaves and Harris1978) that incorporates a $-5/3$ exponent at high frequencies. This distribution is also considered to be representative of moderately homogeneous turbulence in a wide range of flows, including grid-generated turbulence relevant to some wind tunnel tests of turbulence effects on model rotors, turbulence in tidal streams and possibly wake-induced turbulence within large wind farm arrays. The von Kármán spectrum will be assumed in the present paper to represent the turbulence adequately. Mann (Reference Mann1994) discusses other forms of spectra taking into account the anisotropy of the ABL and the ground plane blockage. It is unlikely that the results obtained here using the von Kármán distribution would change greatly if alternatives such as these, for example, the frequently used Kaimal spectrum (Kaimal et al. Reference Kaimal, Izumi, Wyngaard and Cote1972) were used instead.

The von Kármán spectrum representing homogeneous, isotropic turbulence with a substantial inertial subrange, is described by the triple wavenumber energy spectrum function:

where

$g_{0}$ is a ratio of Gamma functions

and the constant $C=(220/9)L_{1\infty }$ .

Hence the triple wavenumber spectrum function of the turbulence is

For spectra involving cross-stream components a simple method to accommodate some effects of anisotropy in ABL turbulence is to replace $\overline{u^{2}}_{1\infty }$ in (2.17) by the relevant turbulence intensity. But for the streamwise component $i=j=1$ :

This triple wavenumber spectrum tensor is integrated with respect to the radial wavenumber $\unicode[STIX]{x1D70F}$ to give the one-dimensional frequency spectrum:

The above integral has been evaluated numerically and the resulting power spectral density (PSD) of the $\unicode[STIX]{x1D705}_{1}$ spectrum for the distorted turbulence in the optimum operating case, $a=1/3$ , i.e. $c=2/3$ is compared with the undistorted upstream turbulence spectrum $\unicode[STIX]{x1D719}_{11\infty }$ in figure 3.

The amplification ratio for large values of $\unicode[STIX]{x1D705}_{1}^{\ast }\rightarrow \infty$ in the small length scale limit is

where from Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik1980), integral result 3.197, $_{2}F_{1}$ is the Hypergeometric function.

Therefore $(\unicode[STIX]{x1D719}_{11}/\unicode[STIX]{x1D719}_{11\infty })(\unicode[STIX]{x1D705}_{1}\rightarrow \infty )=1.033$ when $a=1/3$ . The ratio of the zero frequency intercept value of the spectrum can be obtained by using (2.19) and substituting $\unicode[STIX]{x1D705}_{1}^{\ast }=0$ . The expression for this ratio then simplifies to $c^{-2}$ with value $=2.25$ when $a=1/3$ . Both these asymptotic ratios are consistent with the spectra shown in figure 3.

Figure 3.

Frequency spectrum $\unicode[STIX]{x1D719}_{11}$ for induction factor $a=1/3$ ; solid line distorted, dashed line undistorted.

2.3 Transverse length scale of the distorted turbulence

In this case the magnitudes of the length scales of the turbulence at the rotor change due to distortion by the mean strain field. In a propeller or fan flow field $c>1$ causes extension of $L_{1}$ and contraction of $L_{2}$ (Majumdar & Peake Reference Majumdar and Peake1998) whereas in a turbine flow, $c<1$ , the opposite is the case and the transverse integral length scale $L_{2}$ is increased by the distortion. This potentially increases the correlation of the unsteady forces on the rotor blades and hence the r.m.s. of the integrated unsteady axial force on the rotor as a whole. The transverse length scale $L_{2}$ of the turbulence after distortion is

since the distortion factor in the integral is homogeneous, independent of $y$ in the limit $L_{1\infty }\rightarrow 0$ where $y^{\prime }=c^{-1/2}y$ due to the transverse displacement of the streamlines of the mean flow.

Simplifying by carrying out the integrations with respect to $y$ and $\unicode[STIX]{x1D705}_{2}$ and substituting the integral expression for $\overline{u_{1}^{2}}(0)$ leads to

In this case the result depends on the form of $E(\unicode[STIX]{x1D705})$ since the integrals involving the turbulence energy can no longer be divided out.

Assuming the von Kármán form for $E(\unicode[STIX]{x1D705})$ as before, the integrals for the transverse length scale ratio can be evaluated to give

with $s=c^{-3}-1$ . The transverse length scale amplification ratio $L_{2}/L_{2\infty }$ is plotted against the induction factor $a$ in figure 4. This ratio is smaller than $c^{-1/2}$ the transverse divergence ratio of the mean streamlines. When $a=1/3$ the amplification factor for the transverse length scale ratio is 1.118.

For small length scales the integrated fluctuating force induced on the rotor will be approximately proportional to the fluctuating axial velocity multiplied by the correlation area. The effect of distortion of the incident turbulence on this integrated buffeting force induced by the turbulence on the rotor as a whole is significant, amplifying the r.m.s. axial rotor force amplitude by a factor of approximately

Figure 4.

Transverse length scale amplification ratio $(L_{2}/L_{2\infty })$ versus induction factor $a$ , $L_{1\infty }\rightarrow 0$ .

2.4 Distorted turbulence intensity off the axis

Durbin (Reference Durbin1981) has developed the small length scale analysis for distortion at general points in the flow field of an axisymmetric obstacle. As for the BP analysis, the strain field is assumed to be slowly varying so that the strain terms can be treated locally as constants simplifying the calculation of the distorted turbulent velocity from the distorted vorticity. The strain tensor $\unicode[STIX]{x2202}x_{i}/\unicode[STIX]{x2202}x_{j\infty }$ in this general case contains non-zero off-diagonal terms $(i\neq j)$ which represent tilting and bending of vortex lines in the $(r,x_{1})$ plane. Since the mean axial velocity $\overline{U}_{1}=(1-a)$ is constant over the actuator disc the pure strain terms on the diagonal of the strain tensor are unchanged. The result of this is that the distorted turbulent velocity vector is equal to the value resulting from pure strain alone multiplied by a rotation matrix leading to (Durbin Reference Durbin1981, equation (23))

where $\unicode[STIX]{x1D707}(c)$ and $\unicode[STIX]{x1D708}(c)$ are the distortion factors given in Batchelor (Reference Batchelor1953) for in-line and transverse mean square turbulence velocities. $\unicode[STIX]{x1D707}(c)$ for $c<1$ is given in (2.13). Similarly $\unicode[STIX]{x1D708}(c)$ taking account of the imaginary square root when $c<1$ , is

where $s=c^{-3}-1$ . In (2.25) $\unicode[STIX]{x1D6FD}$ is the rotation angle due to the gradient across the mean flow streamlines of the drift function $\unicode[STIX]{x0394}x_{1}(x_{1},r)$ . This is the $x_{1}$ -displacement of a fluid particle convected by the non-uniform mean flow field relative to the position it would have reached at the same time in uniform flow $\overline{\boldsymbol{U}_{\infty }}(=(1,0,0))$ and is given explicitly in (3.23). $\unicode[STIX]{x1D6FF}=\tan ^{-1}((\sqrt{\overline{U}_{2}^{2}+\overline{U}_{3}^{2}})/\overline{U}_{1})$ is the rotation angle due to the inclination of the mean flow streamlines. For the actuator disc mean flow field, except close to the tip region $(r=1,x_{1}=0)$ , these two rotation angles are in opposite directions, are of comparable size and are small so that the combined rotation angle $(\unicode[STIX]{x1D6FD}+\unicode[STIX]{x1D6FF})$ is small. Evaluating the right-hand side of (2.25) the ratio $\overline{u_{1}^{2}}/\overline{u_{1\infty }^{2}}$ falls gradually with radial distance from the axis and rapidly at the tip as shown in figure 5 for the case when $a=1/3$ .

The radial component of the mean velocity given by actuator disc theory is logarithmically singular at the rotor tip. Therefore the rotation angle $\unicode[STIX]{x1D6FF}$ based on this mean flow field takes the value $\unicode[STIX]{x03C0}/2$ there. The angle $\unicode[STIX]{x1D6FD}$ is finite at the tip but discontinuous between $r<1$ and $r>1$ . The actuator disc flow field has been retained for the present calculations but being unrepresentative in the tip region a tip correction factor is normally applied in BEM theory to cause the blade loading to approach zero appropriately at the tip.

Figure 5.

Radial variation of distorted streamwise intensity for zero length scale and $a=1/3$ .

3.1 Vorticity transport by the mean flow field of an actuator disc

The general case of arbitrary turbulence length scale to rotor radius ratio will now be considered. The vorticity transport equation for unsteady incompressible flow normalised on the mean free stream velocity as before, is

where $\boldsymbol{Z}=\unicode[STIX]{x1D735}\times \boldsymbol{U}$ is the vorticity of the flow and $Re=UR/\unicode[STIX]{x1D708}$ is the mean flow Reynolds number.

Splitting the velocity and other flow variables into mean and fluctuating parts:

and

The flow is assumed to satisfy the requirements for RDT that in terms of the non-dimensionalised quantities

Under these conditions the equation for the fluctuating component of the vorticity, $\unicode[STIX]{x1D73B}$ , can be linearised and outside wall boundary-layer regions viscous effects may be considered negligible. Hence

This equation describes the distortion of the turbulent vorticity field $\unicode[STIX]{x1D73B}$ by the mean velocity field $\overline{\boldsymbol{U}}$ of the rotor, [ $A$ ] by drift and transverse displacement, [ $C$ ] by stretching and tilting, together with two terms involving the mean vorticity field $\overline{\unicode[STIX]{x1D73B}}$ of the rotor being perturbed by the turbulence $\boldsymbol{u}$ , [ $B$ ] by displacement and [ $D$ ] by stretching.

3.2 Mean velocity field of an actuator disc

In order to calculate distortion in the general length scale case spatial integrals of the mean strain field are required and hence the entire mean velocity field $\overline{U}_{i}(x_{j})$ of the actuator disc. This may be obtained as the flow field of a cylindrical vortex sheet extending from the periphery of the disc, to far downstream, and representing the helical vorticity shed from the tips of the rotor blades. The primary component of this vorticity is directed azimuthally, inducing streamwise backflow and radial outflow, and may be represented by a one-dimensional array of vortex ring elements along the rotor axis from $x_{1}=0$ at the rotor to far downstream $x_{1}=\infty$ . The other component of the sheet vorticity lies in radial planes from the axis in the flow direction and induces the swirl in the wake. It does not directly affect the streamwise velocity at the rotor. Expressions for the velocity components induced by an actuator disc derived from this vortex sheet model of the flow have been given by Hough & Ordway (Reference Hough and Ordway1965) in terms of Legendre functions $Q_{-1/2}$ and Heuman’s lambda function or alternatively in terms of integrals of $Q_{-1/2}$ . More recently Conway (Reference Conway1995) derived alternative forms by using the axisymmetric Green’s functions of the flow which are Bessel functions of the first kind. In the simplest case in which the wake vortex sheet is represented as a cylinder of constant radius and constant strength in steady flow, neglecting the downstream expansion, Conway’s formulation gives the streamwise component of the velocity field, including the free stream, as

and

where $\text{J}_{0}$ and $\text{J}_{1}$ are Bessel functions of the first kind.

This cylindrical array of vortex rings or source-dipole discs is equivalent to a source sheet of uniform density covering the actuator disc at the upstream end of the array. The source density $\unicode[STIX]{x1D70E}$ is equal to twice the reduction in axial velocity at the disc, thus

A uniform axial velocity $-\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FF}_{i1}$ within the wake region $r<r_{w}$ , $x_{1}>0$ provides continuity of axial velocity through the disc. This is the model proposed by Taylor (Reference Taylor1944) for a porous plate.

Since the flow field is axisymmetric about the streamwise axis $x_{1}$ it is conveniently expressed in cylindrical coordinates $(r,\unicode[STIX]{x1D703},x_{1})$ and can be obtained by carrying out the integration of the velocity induced by the source elements $\unicode[STIX]{x1D70E}r\unicode[STIX]{x1D6FF}r\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D703}$ with respect to $\unicode[STIX]{x1D703}$ . The following integral expression is then obtained for the axial velocity $\overline{U}_{1}$ :

for $x_{1}<0$ all r, or for $x_{1}>0$ with $r>r_{w}$ , and

for $x_{1}>0$ if $r<r_{w}$ .

and $r_{w}=1$ here is the radius of a wake cross-section. Here $E$ is a complete elliptic integral of the second kind. This alternative formulation to that of Conway (Reference Conway1995) requires numerical integration only over a finite interval (0,1). The streamlines of the mean flow when $a=1/3$ are shown in figure 1.

The mean velocity on the axis $r=0$ can easily be integrated in either formulation to give

3.3 Distortion of the incident turbulence by the mean flow

The effect of the two terms [ $A$ ] and [ $C$ ] representing stretching and displacement are now considered following the Lagrangian approach as taken by Hunt (Reference Hunt1973), to evaluate the distorted turbulent vorticity field. This is then substituted into the Biot–Savart integral relationship between velocity and vorticity to compute the distorted velocity field. As stated previously the turbulent axial velocity will be evaluated at the disc centre (0, 0, 0) only to give a measure of the effect of the distortion on the streamwise velocity component which is the one that most affects the unsteady loading of the rotor.

Computation of the complete unsteady loading of a rotor would require evaluation of the distorted streamwise velocity component of the turbulence (and if appropriate also the azimuthal component) at a range of points along a radius of the swept actuator disc. Away from the centre of the disc most of the simplifications due to axisymmetry can no longer be applied hugely increasing the computational effort required. Therefore in the present paper which focuses on assessing the relative magnitude of the distortion effect the full computation has only been carried out for the velocity at the centre of the disc (0,0,0) and estimates have been made of the variation of distortion over the rest of the disc as a function of radial distance from the centre.

Since the RDT analysis is linear it is possible to consider individually the distortion of each Fourier component of the incident turbulence

in which

and then integrate the result over all wavenumbers. In (3.18) $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D705}_{1}$ since the turbulence far upstream is assumed to be convected frozen by the free stream $\overline{\boldsymbol{U}}=(1,0,0)$ .

In order to calculate the distorted turbulent velocity component $u_{1}$ at the disc it is convenient to represent the turbulent vorticity field as rings of vorticity concentric about an axis parallel to the $x_{1}$ axis and passing through the point at which the velocity is to be calculated (figure 5). Using a cylindrical coordinate system $(r,\unicode[STIX]{x1D703},x_{1})$ , with this axis and $x_{1}=0$ at the disc, the vorticity field is expressed at each point of each ring element as radial, azimuthal and axial components $(\unicode[STIX]{x1D701}_{r},\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}},\unicode[STIX]{x1D701}_{x1})$ . Using this system the axial velocity $u_{1}$ (0,0,0) depends only on the azimuthal vorticity component $\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}}$ . However any given ring in the flow field will have arrived at its position from its initial state far upstream after convection and distortion by the mean flow field. If its end state is circular its initial state in general cannot have been so because of the non-uniform distortion in transit. The exception to this is when the mean flow field is also coaxial with the $(r,\unicode[STIX]{x1D703},x_{1})$ axis which is the case when the velocity is to be calculated at the centre of the disc. The rings of vorticity are then only subject to radial expansion and axial displacement as a whole as they are convected. The initial and final states are both circular and are simply related by a radius ratio $r/r_{\infty }$ and their axial $x_{1}$ coordinates. The $\unicode[STIX]{x1D701}_{r}$ and $\unicode[STIX]{x1D701}_{x1}$ components will both be stretched and subject to displacement from their undistorted positions, but since, although distorted, they remain in the $[r,x_{1}]$ plane, neither of them contribute to the axial velocity induced at (0,0,0). The $\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}}$ component is subject to stretching because the mean streamlines diverge, and to radial displacement for the same reason. It is also subject to streamwise displacement due to the slowdown of the streamwise velocity. Although $\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}}$ varies with $\unicode[STIX]{x1D703}$ around each ring the axisymmetry of the geometry and distortion greatly simplifies the computation and this therefore is the main case which has been computed.

When the velocity is calculated at other points of the disc, radially removed from the centre, the mean flow field and hence the distortion are no longer concentric with this axis system. Since the initial state of a ring whose end state is circular will now not be circular it will include $\unicode[STIX]{x1D701}_{r}$ and $\unicode[STIX]{x1D701}_{x1}$ contributions making for much greater computational work. For these cases therefore approximate estimates of $u_{1}$ will be calculated.

The actuator disc model omits the mean swirl component in the wake which is small except at low tip speed ratios. The helical twisting of the $\unicode[STIX]{x1D701}_{x1}$ component by this swirl is therefore a relatively small effect. It only occurs in the wake, being zero everywhere outside it, and in particular upstream of the rotor. It generates a very small contribution to $\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}}$ and hence to the axial velocity component of distorted turbulent velocity at the disc centre and this effect is omitted throughout. The whole of the distortion effect computed for the centreline streamwise fluctuating velocity and the majority of it for the velocity off the centreline, is due to the expansion of the rings of vorticity as they are convected towards and past the actuator disc together with a lag effect due to drift. The former intensifies the azimuthal vorticity by stretching it and the latter leads to phase changes.

3.4 Distortion of the turbulence impacting the centre of the actuator disc

The components $\unicode[STIX]{x1D701}_{2\infty }$ and $\unicode[STIX]{x1D701}_{3\infty }$ in the cross-flow plane $x_{1}=-\infty$ of the incident vorticity can be expressed as radial $\unicode[STIX]{x1D701}_{r\infty }$ and azimuthal $\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}\infty }$ components at any point $(r_{\infty },\unicode[STIX]{x1D703})$ where $\unicode[STIX]{x1D703}=0$ , arbitrary when axisymmetric, is taken as aligned with the direction of the radial wavenumber $\unicode[STIX]{x1D70F}_{\infty }$ . Substituting from (3.18) into

where $\unicode[STIX]{x1D6FD}$ is the angle between the radial wavenumber $\unicode[STIX]{x1D70F}_{\infty }$ and $\unicode[STIX]{x1D705}_{2\infty }$ i.e. the $x_{2}$ axis, gives

$\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}}$ at a general point $(r,\unicode[STIX]{x1D703},x_{1})$ is therefore related as a result of stretching and radial displacement by the mean flow field of the disc to its far upstream value $\unicode[STIX]{x1D701}_{\unicode[STIX]{x1D703}\infty }$ by

where $r_{\infty }$ is the radial coordinate far upstream of the streamline passing through $(r,\unicode[STIX]{x1D703},x_{1})$ . Applying continuity of volume flow from far upstream $(-\infty )$ within the streamtube of radius $r$ at $x_{1}$ ,

Combined with this radial displacement from $r_{\infty }$ to $r(x_{1})$ there is an axial displacement which changes the phase of the Fourier component arriving at $x_{1}$ by $-\text{i}\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x0394}x_{1}$ , where

where $\unicode[STIX]{x1D6F9}(r)$ is the mean flow streamline passing through $r(x_{1})$ .

Figure 6.

Distortion of components of the turbulent vorticity by the mean flow: 1  $=$  streamwise, compressed and rotated in radial plane; 2  $=$  radial, stretched and rotated in radial plane; 3  $=$  azimuthal, stretched; $A$ is the mean flow axis for $u_{1}$ computation at disc centre; $B$ is the shifted mean flow axis for $u_{1}$ at $r^{\prime }\neq 0$ (mean streamlines not indicated for this case).

This is Lighthill’s drift function discussed by Hunt (Reference Hunt1973). In the case of a planar two-dimensional body which has mean drag the far upstream velocity field is like that of a line source in a uniform flow for which the upstream perturbation only decays as $(-x_{1})^{-1}$ . As a result the integral forming the drift function cannot be evaluated simply and the associated phase change requires special treatment. However for an axisymmetric source-like flow the perturbation decays as $(-x_{1})^{-2}$ , the integral is convergent for all finite $x_{1}$ and a finite phase change is given directly by (3.23). (Equation (3.23) is not convergent for $x_{1}\rightarrow \infty$ within the far wake because $\overline{U}_{1}\rightarrow (1-2a)$ , but this is a limited region $r<1$ and has negligible effect at the rotor plane.)

Applying the Biot–Savart integral, the axial turbulent velocity induced at the disc centre $0$ by the distorted vorticity field is

Converting this into $(r,\unicode[STIX]{x1D703})$ coordinates:

where $p=\text{i}\unicode[STIX]{x1D70F}\hat{u} _{1\infty }-\text{i}(\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D705}_{2}/\unicode[STIX]{x1D70F})\hat{u} _{2\infty }-\text{i}(\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D705}_{3}/\unicode[STIX]{x1D70F})\hat{u} _{3\infty }$ and $q=\text{i}(\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D705}_{3}/\unicode[STIX]{x1D70F})\hat{u} _{2\infty }-\text{i}(\unicode[STIX]{x1D705}_{1}\unicode[STIX]{x1D705}_{2}/\unicode[STIX]{x1D70F})\hat{u} _{3\infty }$ .

Carrying out the $\unicode[STIX]{x1D703}$ integration using the results that

and

from the condition for divergence-free flow, therefore at the centre of the disc

In (3.30) $f(r,x_{1})$ is the stretching factor $r/r_{\infty }$ and $\text{J}_{1}$ is a Bessel function of the first kind.

In the event that the rotor is absent or equivalently the induction factor $a=0$ there will be no distortion of the incident turbulence. In that case $f=1$ , $\unicode[STIX]{x0394}x_{1}=0$ and making use of integral results (Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik1980), integrals 3.771(2), 6.521(2))

from which

which is consistent with there being no distortion.

The main difficulty with the Biot–Savart integration method is the slow convergence of the radial integral when the wavenumber $\unicode[STIX]{x1D70F}$ is small. Computation of this can be speeded up by subtracting (3.32) from (3.30) and therefore only numerically integrating the difference between the distorted and undistorted vorticity. This integrand $\rightarrow 0$ rapidly as $r\rightarrow \infty$ so that the radial integral converges much more rapidly.

3.5 Computation of spectra and intensities of the distorted turbulence

Having evaluated $\hat{u} _{1,0}$ from (3.30) above, velocity spectra and intensities of the distorted turbulence at the disc centre can be computed by integrating over the wavenumbers $\unicode[STIX]{x1D705}_{2}$ and $\unicode[STIX]{x1D705}_{3}$ , as for the small length scale case

where $\tilde{()}$ indicates a complex conjugate. Converting to radial $(\unicode[STIX]{x1D70F})$ integration, taking account of the axisymmetry of $I_{D}$

and integrating with respect to $\unicode[STIX]{x1D705}_{1}$

Evaluation requires substitution of an expression for the energy function $E(\unicode[STIX]{x1D705})$ and the von Kármán distribution will be used for this as before. Therefore,

and

Figure 7 shows examples of the spectrum $\unicode[STIX]{x1D719}_{11}(\unicode[STIX]{x1D705}_{1})$ computed for the optimum value of the induction factor $a=1/3$ and a range of values of the length scale ratio from 0.1 to 10.0, the curves for 0.1 and 0.3 being hardly distinguishable in the log plot.

Figure 7.

PSD of $u_{1}(0)$ versus frequency $k_{1}^{\ast }$ non-dimensionalised by $L_{1\infty }$ of turbulent axial velocity at the centre of the actuator disc. Distorted (solid lines) $L_{1\infty }=0.1,0.3,1.0,3.0,10.0$ as indicated top to bottom at left (0.1 and 0.3 curves are indistinguishable), undistorted (dashed line).

The PSD of both the undistorted and distorted spectra are approximately flat for low frequencies, $\unicode[STIX]{x1D705}_{1}<O(1)$ , followed by a $\unicode[STIX]{x1D705}_{1}^{-5/3}$ asymptotic behaviour as $\unicode[STIX]{x1D705}_{1}^{\ast }\rightarrow \infty$ . Therefore, the effect of distortion on the lower frequency range of the spectrum is characterised approximately by the ratio of the distorted to undistorted zero frequency intercepts:

The ratio $\unicode[STIX]{x1D6EC}$ can be computed relatively easily by setting $\unicode[STIX]{x1D705}_{1}=0$ in (3.37) and is plotted in figure 8 against length scale $L_{1\infty }$ . Its value when $L_{1\infty }\rightarrow 0$ is consistent with $\unicode[STIX]{x1D6EC}=(1-a)^{-2}$ given by the BP homogeneous strain analysis for $L_{1\infty }=0$ . It is notable in this figure that $\unicode[STIX]{x1D6EC}$ increases initially with increasing $L_{1\infty }$ from its value at $L_{1\infty }=0$ reaching a maximum at $L_{1\infty }\simeq 0.3$ and then decreasing thereafter to slowly approach the large length scale asymptotic value of 1.0.

Figure 8.

Variation of zero frequency intercept ratio, distorted to undistorted, against $L_{1}$ .

In the same way as was done in the small length scale case the intensity of the distorted to undistorted turbulence can be evaluated by integrating the spectra over $\unicode[STIX]{x1D705}_{1}$ . The results of this are shown in figure 9. Also the ratio of the transverse length scales, $L_{2}/L_{2\infty }$ , of distorted to undistorted streamwise velocity of the turbulence can similarly be calculated, requiring integration over $\unicode[STIX]{x1D705}_{1}$ and one transverse wavenumber while setting the other to 0, or equivalently formulating the integral in terms of the transverse radial wavenumber $\unicode[STIX]{x1D70F}$ . However evaluation of $L_{2}$ requires spatial integration of a correlation function involving the distortion integral $I_{D}$ , equation (3.30), over the radius $r$ in the disc plane. Evaluation of $I_{D}(r)$ cannot then make use of the axisymmetric simplification and as a result is much more expensive computationally. In order to retain simplicity in assessing the magnitude of $L_{2}$ here, $I_{D}$ and $\unicode[STIX]{x0394}x_{1}$ have been approximated by their values at the disc centre, but the factor $\text{e}^{\text{i}\unicode[STIX]{x1D70F}r/f}$ which is the strongest influence on the correlation at given $y$ (or $r$ ), is not approximated in this way. As before for general length scale cases both the distorted intensity and transverse length scale depend on the form of the undistorted spectrum but it is unlikely that the results would change significantly if the assumed spectral distribution were changed.

Evaluating the transverse length scale at the disc centre:

The $y$ integration is infinite when $\unicode[STIX]{x1D705}_{2}=0$ , giving a Dirac delta function of $\unicode[STIX]{x1D705}_{2}$ which is integrated separately. The finite part of the integral is then integrated by converting it to axisymmetric form. Hence, using the von Kármán form for the energy distribution:

with

Here

and $\text{J}_{0}$ is a Bessel function of the first kind.

In the limit $L_{1\infty }\rightarrow 0$ the integrals for the spectra, intensity and transverse length scales, equations (3.37), (3.38) and (3.41) take their dominant contribution from the region $\unicode[STIX]{x1D70F}^{\ast },\unicode[STIX]{x1D705}_{1}^{\ast }\sim 1$ and therefore $\unicode[STIX]{x1D70F},\unicode[STIX]{x1D705}_{1}\sim \infty$ . In that case the double integral $I_{D}(\unicode[STIX]{x1D70F})$ given in (3.31) can be shown to take the value

and as a result, with $s=c^{-3}-1$ ,

as $L_{1\infty }\rightarrow 0$ which agree with the values obtained in the $L_{1\infty }=0$ BP analysis.

Figure 9 shows a plot of $L_{2}/L_{2\infty }$ for several values of the integral length scale $L_{1\infty }$ when the induction factor $a=1/3$ . The increase in $L_{2}$ due to this distortion is therefore small ( ${<}5\,\%$ when $L_{1\infty }>1.0$ ) with $a=1/3$ . This would seem to agree with measurements of the transverse correlation coefficient $R_{11}(y)$ by Hancock & Pascheke (Reference Hancock and Pascheke2014) who found little measurable change for turbulence with $L_{1\infty }>1.0$ .

Figure 9.

Variation of $(\overline{u_{1}^{2}}/\overline{u_{1\infty }^{2}})$ (o) and $(L_{2}/L_{2\infty })$ (+) with integral length scale $L_{1\infty },a=1/3,r=0$ .

3.6 Distorted turbulence intensity off the axis

Evaluating the distorted turbulent velocity at a general point $(r^{\prime },\unicode[STIX]{x1D703}^{\prime },0)$ away from the axis the simplifications resulting from axisymmetry can no longer be applied. As a result the distortion of all components of the incident vorticity must be considered in the full three-dimensional space and the computations are correspondingly orders of magnitude more intensive. As an alternative here an estimate has been made of the distortion of the streamwise component of the turbulence away from the axis. In the analysis for the small length scale limit in § 2.4 it was noted that the combined rotation angle $\unicode[STIX]{x1D6FD}+\unicode[STIX]{x1D6FF}$ for the mean flow through the actuator disc was small except close to the tip region and therefore in making this estimate the effect of this angle on tilting the vorticity components is ignored. It is then convenient if the $(r,\unicode[STIX]{x1D703},x_{1})$ axis system used in § 3.3 for axisymmetric distortion is retained passing through the selected point $(r^{\prime },\unicode[STIX]{x1D703}^{\prime },0)$ of the actuator disc where the distorted streamwise turbulent velocity is to be estimated. In this axis system the actuator disc itself and its mean flow field which are concentric with axis B in figure 6, are therefore shifted radially by an amount $r^{\prime }$ off the $r=0$ axis in a direction which can be taken to be $\unicode[STIX]{x1D703}=0$ since the distortion is independent of $\unicode[STIX]{x1D703}^{\prime }$ and the velocity point is retained at (0,0,0). As before, only components of vorticity which are directed azimuthally around rings concentric with the $r=0$ axis contribute to the velocity at (0,0,0). But each such ring in the vicinity of the disc will now have been subject to a mean strain field no longer concentric with the ring’s axis which therefore varies with $\unicode[STIX]{x1D703}$ around the ring. This strain field introduces contributions from $\unicode[STIX]{x1D701}_{r}$ and $\unicode[STIX]{x1D701}_{x1}$ and phase differences. These will be ignored here for the purposes of making an estimate, since the effect of the combination of tilting and drift gradients was found to be small in the small length scale case except in the tip region, § 2.4. Also it will be assumed that the velocity induced by the azimuthal vorticity can be satisfactorily estimated from the average around each ring of the pure strain component of the strain tensor to which it has been subjected. Thus, in (3.30) $f(r,x_{1})$ is replaced $\overline{f}(r,x_{1})=(1/2\unicode[STIX]{x03C0})\int _{0}^{2\unicode[STIX]{x03C0}}f(r,\unicode[STIX]{x1D703},x_{1})\,\text{d}\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x0394}x_{1}(r,x_{1})$ is replaced by $\overline{\unicode[STIX]{x0394}x_{1}}(r,x_{1})=(1/2\unicode[STIX]{x03C0})\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x0394}x_{1}(r,\unicode[STIX]{x1D703},x_{1})\,\text{d}\unicode[STIX]{x1D703}$ the average values of the stretching function $f$ and the drift $\unicode[STIX]{x0394}x_{1}$ around each ring. Therefore,

Equation (3.47) is essentially the same as (3.30) except that the functions $f$ and $\unicode[STIX]{x0394}x_{1}$ are replaced by precalculated ring averages over $\unicode[STIX]{x1D703}$ . The computations may then be repeated to provide estimates of the distortion of the streamwise turbulent intensity at points along a radius of the actuator disc converging to the previously computed value at the centre. Estimates at four radial distances from the centre have been calculated by this procedure when $a=1/3$ and are shown in figure 10. There is a gradual reduction in intensity with increasing radius, with an increased rate of fall off as the tip region is approached.

Figure 10.

Radial variation of distorted streamwise intensity for $L_{1\infty }/R=1.0$ and $a=1/3$ .