Experimental Study and Simulation of a Small-Scale Horizontal-Axis Wind Turbine

With global installed wind energy capacity nearing 500 GW and expected to increase to 800 GW by 2020, the need to improve the wind farm efficiency is imperative. Wind farm power output is a function of many variables including atmospheric conditions, geographic terrain, wind turbine design, turbine spacing, and electrical transmission. Understanding the behavior of turbulence generated from wind turbines and wind turbine wake dynamics can lead to more robust wind turbine designs, assist engineers in wind farm layout, and lead to increased wind farm efficiency.

Computational fluid dynamics (CFD) has become a staple in the study of wind turbine wakes, and advances in wind energy technology can be directly attributed to research efforts using CFD in conjunction with wind tunnel experiments. There is an overwhelming amount of literature available regarding computational techniques for modeling wind turbine wakes and applying turbulence closure methods as they relate to horizontal-axis wind turbines (HAWT). Likewise, there are vast publications for wind tunnel experiments to study wake dynamics and wind turbine performance. Research using Reynolds stress model (RSM) for wind turbine simulations is limited to studies where the turbine is modeled as a rotating disk. Studies using RSM to model a fully resolved wind turbine are nonexistent and have likely been avoided due to the computational expense and potential for instabilities. Periodically, researchers have summarized advances and state-of-the-art approaches to wind turbine modeling or wake simulations. Several of these reviews are noteworthy in that they are comprehensive, relevant, and detailed. In 1999, Crespo et al. [1] provided an overview of the different modeling methods used to predict velocity deficit in the wake of single and multiple wind turbines. Their review included discussions on kinematic wake models, field models, terrain effects, and wind farm modeling. In kinematic models, it is assumed that the perturbation profiles of the velocity deficit and turbulence intensity are axisymmetric and follow a self-similar distribution. Kinematic models express the velocity deficit by an analytical expression developed from theoretical work on co-flowing jets and experimental data. Field models are much like today's CFD models in that they calculate the velocity at every point of the flow field and rely on a numerical solution of turbulent transport equations. Early kinematic and field models are still incorporated into software used for wind farm analysis. Vermeer et al. [2] followed in 2003 with an overview of computational methods about horizontal-axis wind turbines and included further discussion of kinematic and field models. The unique aspect of the Vermeer paper was their segregation of experimental and analytical research based on near and far wake studies. They also included a thorough review of experiments that had been performed on a variety of wind turbines. Hansen et al. [3] studied computational methods for wind turbine analysis including blade element momentum (BEM) methods, panel methods, vortex methods, and actuator disk methods. Aeroelastic methods for predicting the dynamic response of the turbine blades from time-dependent aerodynamic loads were also presented. In 2011, Sanderse et al. [4] provided a state-of-the-art review of CFD methods for simulating wind turbines. They classified the different numerical methods used and distinguished between models specific to simulating the rotor versus simulating the wake. For simulating the wake, kinematic and field models, previously discussed, are employed. For simulating the rotor, the body forces in the momentum equations can be represented by an actuator disk (AD), actuator line (AL), or actuator surface (AS).

In the case of the actuator disk with uniform loading, the body force is expressed as a function of the approaching relative velocity and the thrust coefficient. For a nonuniformly loaded disk, the body force is dependent on the radial location but is constant over an annular area as seen in Fig. 1. Sectional lift and drag coefficients are used to determine local forces, and the body force is resolved as a time average of circular line forces. The actuator line method is similar to the computation of local line forces; however, the line forces are time-dependent. Thus, the AL method allows the tip vortex shedding to be calculated. The actuator surface method is also time-dependent, but more complex as the planar surface model forces are determined from pressure and skin friction distribution as well as from lift and drag coefficients. For all nonuniform loading situations, the approaching velocity is computed by interpolating the velocity field at cell locations near the disk surface. Other researchers have provided numerical reviews, but not to the depth as the aforementioned papers.

In more recent studies where the objective is to determine the turbine power output, validate a solver code, evaluate a meshing technique, or validate a turbulence model, the blade is usually modeled directly with a higher cell density near the blades. Ma et al. [5] used a commercial code, ansys fluent, to compare an unsteady delayed-detached eddy simulation (DDES) model to steady Spalart–Allmaras and k–ω turbulence closure models. Their simulations were evaluated against experimental data from wind tunnel tests on a Fortis Montana 5.8 kW HAWT and they concluded the time-averaged Reynolds averaged Navier–Stokes (RANS) methods tended to overestimate the power losses from the flow field around the turbine compared to the DDES. Sagol et al. [6] also utilized ansys fluent to perform a similar study, but focused on applying k–ε and k–ω methods to compare to NREL Phase VI data. Carrion et al. [7] used a unique compressible multiblock solver to simulate wind turbines tested in the large-scale, low-speed, facility of DNW (German–Dutch Wind Tunnels). This testing was commonly referred to as the “model experiments in controlled conditions” (MEXICO). Carrion provided a detailed list of all known CFD studies of the MEXICO experiments along with the solver, turbulence models, and geometry used in each study. AbdelSalam and Ramalingam [8] employed ansys fluent to create a full rotor model of a Danwin 180 kW turbine and used a standard k–ε turbulence model to compare the wake velocity profile at various downstream locations to widely scattered experiment data. The results indicated that the full rotor model predicted the downstream wake more accurately than the actuator disk models. Rocha et al. [9] used an open source CFD package, openfoam, with a k–ω sheer stress transport (SST) turbulence model to study the performance of small (3 m diameter) prototype, rotor and compared the result to experimental values for the power coefficient. Inlet conditions based on the turbulent length scale were shown to provide unrealistic power coefficients for certain values of the length scale. Tran et al. [10] used the commercial software star ccm+ from CD-Adapco to predict the wake for a floating offshore wind turbine using a k–ω SST turbulence model. An overset mesh, sometimes referred to as overlapping or Chimera grid, was incorporated with rigid body motion (RBM) used to simulate the motion of the blade. An alternate approach was also investigated using ansys fluent and a multiple reference frame model and the results of both simulations were compared to established BEM codes, fast, and ubem. Lawson et al. [11] also used star ccm+ to simulate an off-shore HAWT and compared the results to BEM code WT_Perf using a hypothetical wind turbine design. CFD simulations with multiple turbines modeled directly are virtually nonexistent except for Seydel and Aliseda [12]. They modeled two offset NREL Phase VI rotors with ansys fluent using a k–ω SST turbulence model.

Over the past decade, there has been a growth in the number of wind tunnel experiments on small-scale wind turbines. The growth can be attributed to the desire to study wake effects under controlled conditions and the need to provide data for CFD validation. Experimental setup varies depending on the size of the wind tunnel, instrumentation, and test objective. Aubrun et al. [13,14] experimented with a small, three-bladed, 416 mm diameter wind turbine using hot wire anemometry (HWA) to measure the velocity profile at certain distances downstream. The wind tunnel incorporated an upstream mesh to simulate the ABL. For comparison, a nonrotating porous disk was also tested where the disk was designed to have the same velocity deficit at 0.5D_b. The results indicated the porous disk created similar downstream turbulence characteristics, and spectral analysis of the HWA data indicated the tip vortex structure was indistinguishable at x/D_b > 3. Medici [15] performed wind tunnel experiments on several small turbine models with different blade designs. The experiments were aimed at studying the variation in power output as a function of turbine yaw. Power was determined by the electrical output from a small motor acting as a generator. Using the same equipment and test facilities, additional work performed by Medici [16] showed that the presence of a wind turbine affects the upstream flow field more than 3D_b upstream of a wind turbine. Hot wire anemometry has been used to measure velocity profiles and wake turbulence for single small-scale wind turbines and on arrays of wind turbines. Turbulence measurements using particle image velocimetry (PIV) on a single turbine have been captured by Hu et al. [17], Yang et al. [18], and Massouh and Dobrev [19]. PIV was also used to measure turbulence around an array of small turbines by Cal et al. [20] and Lebron et al. [21].

CFD studies implementing the Reynolds stress model (RSM) turbulence closure model for a wind turbine are rare and use an actuator disk model to simulate the wind turbine [22,23]. The majority of CFD studies that model the rotor directly focus on two-equation methods based on the Boussinesq hypothesis. The wide use of two-equation models is mostly because of the high computational expense associated with the RSM or any higher-order models. With faster and more powerful computers, LES methods are becoming more demanded, especially in conjunction with generalized actuator (AD, AL, AS) methods.

Over the past decade, there has been growth in the number of wind tunnel experiments on small-scale wind turbines. The growth can be attributed to the desire to study wake effects under controlled conditions and the need to provide data for CFD validation. Experimental setup varies depending on the size of the wind tunnel, instrumentation, and testing objective. The goal of this study was to examine the wake profile and turbulence downstream and between multiple small-scale turbines where the turbine spacing was considered in the near-wake region. The experimental data were to be used for comparison to CFD simulations.

The wind tunnel in the Wind Energy Laboratory at the University of Wisconsin-Milwaukee (UWM) was used for the experimental study. Welsh [24] provides a detailed description of the wind tunnel facility and the design considerations, but a general synopsis is provided here. The wind tunnel at UWM is an open-circuit, suction-type tunnel that uses an axial fan to draw air through the test section. The inlet settling chamber section was designed with a 7.62 cm long honeycomb and 1 cm hexagonal shaped cells to reduce the large-scale turbulence and to eliminate mean lateral and vertical velocity components. A series of screens with reduced mesh sizes was incorporated to reduce further the turbulence and reduce the variation in the mean longitudinal velocity. The room is temperature controlled which reduces the variation in temperature during testing. Temperature changes with the wind turbine operating continuously for up to 10 h were less than 0.5 °C. The inlet to the contraction section is a little over 9.3 m² and transitions to the 1.4 m² test section in ∼4 m distance. The 6.2 contraction ratio is on the low end of the recommended range, but given the low wind tunnel speed, it was deemed acceptable. The test section is ∼120 cm × 120 cm × 243 cm long and has clear polycarbonate walls to provide a smooth surface. The wind tunnel boundary layer was not measured for this experimentation, but it was estimated to be less than 15 cm in 2.4 m so that the boundary layer does not extend into the turbine wake. The diffuser section was designed with a 2.25 expansion ratio and a 2.5 deg expansion angle. The diffuser section transitions from a square cross section at the exit of the test section to an octagonal shape at the entrance to the 1.83 m fan diameter. The six-bladed fan is attached to a 34 hp motor which is controlled by a variable frequency motor speed drive mounted to the side of the fan enclosure.

Inside the wind tunnel test section, a three-axis traverse system is mounted on the top panel. Stepper motors are attached to each arm to allow the hot wire probe to be positioned accurately upon command. Each stepper motor provides 25.4 mm of movement for each 4000 revolutions giving a 6.35 μm resolution on the position of the hot wire probe. The stepper motors were driven by a Velmex VXM controller and communication between the data acquisition system and the controller was performed via RS232. The baseline (empty tunnel) turbulence intensity was measured at 0.32%.

Three small-scale wind turbines were built for testing in the UWM wind tunnel. Each wind turbine included a custom three-blade rotor made from ABS plastic, a tower made from 12.7 mm diameter steel rod, a 30.5 cm × 30.5 cm × 6.4 mm steel base plate, and a small DC motor (Radio Shack 273-0258) to act as a generator. Figure 2 shows a sketch of the assembly.

The hub height, H, was designed at 30.5 cm and the blade diameter, D_b, was designed at 20.3 cm giving a hub height to blade diameter ratio of 1.5. The overall size of the wind turbines was based on the size of the wind tunnel test section and a desire to fit two turbines side by side, if necessary, with minimal effect from the sides of the wind tunnel. The hub height to blade diameter ratio was selected based on a comparison to several different full-scale wind turbines from various manufacturers. The blade design used for the experiments was arbitrary; however, the blade was designed so that lift and drag coefficients could be determined from published data or airfoil codes. The three-blade model, shown in Fig. 3, was created using pro-engineer software and a.stp file of the model was used to create the blades on a 3D rapid prototyping machine. Each blade had a linear taper with a twist angle ranging from 11 deg to 25 deg along the blade length.

The turbine blade cross section was based on the NACA 4424 profile shown in Fig. 4 with the section thickness scaled from the chord length. Thus, the cross section near the hub was thicker than near the blade tip. The NACA 4424 profile was selected because the coordinates are established and it offered a thicker blade profile relative to the chord length to improve the strength and prevent the blades from breaking during handling.

Hot wire anemometry (HWA) has been used for many years to measure velocity fluctuations in turbulent flows. It is well suited for turbulent flow studies in low Reynolds number air/gas flows with low to moderate turbulent intensities (<25%). It has the advantage of being much lower cost compared to particle image velocimetry (PIV), laser Doppler anemometry (LDA), or planar laser-induced fluorescence (PLIF) systems and does not require particulate in the flow that could contaminate the components of a wind tunnel. HWA is relatively easy to use, provides a continuous analog output, meaning no information is lost, and has an high temporal resolution for spectral measurements; turbulence fluctuations as high as 400 kHz can be measured in principle, HWA is used to measure local velocity by placing a heated, fine wire into a flow stream and controlling the wire current that corresponds to the convective heat loss in the wire. Thus, HWA has two components: a probe containing the heated wire and an anemometer which contains the necessary electronic components to control the wire current (i.e., wire temperature) and produce a conditioned voltage signal that is proportional to the velocity. Hot-wire probes are available, commercially, in a variety of styles and can include multiple wires to measure the velocity in multiple axes. In general, the number of wires in the probe corresponds to the number of velocity components being measured. The hot-wire probe used for this study was a Dantec Model 55P64, dual sensor, cross-wire (X-wire) type probe and is capable of measuring U and V components of the velocity vector. The 55P64 has two platinum-plated tungsten, 5 μm diameter, wires welded to the probe at 45 deg to each other and can measure velocity components within a ±45 deg cone. The probe is capable of measuring velocity between 0.05 m/s to 500 m/s. Two Dantec 54T30 miniature constant temperature anemometers (CTAs) were used to provide a 0–5 V analog output voltage based on the characteristics of the probe. The most important aspect of velocity measurement using hot-wire anemometry is the calibration of the probe and anemometer. Experimental work performed over the years has led to a number of different methods for calibrating hot wire anemometers with inclined wires. The literature on hot wire anemometry is extensive and covers a wide variety of topics. Browne et al. [25] provided a comprehensive list of methods developed to account for longitudinal cooling and yaw angle sensitivity. One method that has not been researched extensively is the use of a biharmonic spline interpolation algorithm, developed by Sandwell [26] for constructing ocean floor topography maps from satellite data. The algorithm has also been used for mapping the total electron content of the ionosphere from global positioning system data [27]. The algorithm is well suited for hot wire anemometry where measurements at different yaw angles and calibration velocities could be irregularly spaced. The algorithm has been implemented in the matlabgriddata.m function with the “-v4” method option. Only a few researchers have implemented the griddata.m function for hot wire anemometry. For this research, the biharmonic interpolation algorithm, sometimes referred to as a surface interpolation, was compared to traditional effective angle and polynomial surface fit techniques for hot wire calibration. The three different methods for converting hot wire anemometer voltages to velocities were examined for accuracy and computation time. Figure 5 shows the surface interpolation method using the griddata.m function was more accurate when applied to the same calibration data. The only drawback was the time to convert voltages to velocity was more than 40 times longer. All hot-wire measurements for this study used the griddata.m function to convert voltages to velocity after compensating for temperature.

Before discussing the experimental results, it is useful to review the aerodynamics of a horizontal-axis wind turbine wake. The wake is typically divided into a near wake and a far wake as shown in Fig. 6.

The near wake is defined as the area directly behind the rotor extending to 1–3 rotor diameters. This region is dominated by pressure and velocity gradients resulting from the extraction of mean flow energy by the rotor. In this part, the blade geometry dictates the shape of the flow field, and the pressure gradient at the rotor is important in developing the wake velocity deficit. The reduction in velocity in the near wake region is directly related to the thrust coefficient since this determines the momentum transferred from the free stream to the turbine. As the lower speed wake convects downstream, the velocity gradient between the wake and the free stream creates shear-generated turbulence, which transfers momentum into the wake and causes mixing. The mixing region spreads to the center of the wake and outward which erodes the velocity deficit and expands the wake. The wake velocity eventually recovers and the profile becomes axisymmetric and Gaussian in shape.

The wind tunnel at UWM was used to measure the velocity profiles and turbulence downstream of a model wind turbine. All hot-wire measurements were taken with the fan motor drive frequency set at 23 Hz which corresponded to 6.65 m/s. At this wind speed, rotor torque and thrust were estimated near 0.005 N·m and 0.3 N, respectively. Experimental uncertainty for velocity measurements was established as 2.83%.

Hot-wire measurements with the single turbine were taken in 6.35 mm increments in both vertical (y) and horizontal (−z) directions at the planes shown in Fig. 7. The measurement planes have been identified as S0–S23 and dimensions for plane locations are listed in Table 1. Except at the planes immediately behind the rotor (S8–S11), vertical measurements were taken at the turbine centerline 20.3 cm below the hub height and extending 20.3 cm above the hub height. The horizontal measurements were taken at the hub height extending laterally 20.3 cm from the centerline. Lateral (horizontal) measurements were restricted to the half-plane to avoid interference between the turbine and the traverse. Directly behind the turbine, vertical measurements were taken from 1.27 cm above the rotor centerline extending to 20.3 cm; horizontal measurements were taken 2.54 cm from the rotor centerline extending 22.9 cm. At each measurement point, the hot-wire voltages were sampled for 0.5 s at 20 kHz and data for each point were saved to a unique text file. Before initiating hot-wire measurements in each plane, the blade speed was recorded with the hot wire probe positioned at the hub height. The arm of the traverse was located approximately 10 cm downstream of the hot-wire probe and approximately 20 cm from the edge of the traverse arm. It was noted that the blade speed was dependent on the position of the traverse. The air temperature was recorded at each measurement plane and values were used to correct for the difference between measurements and hot-wire calibration. Temperature changes were less than 0.4 °C throughout all measurements.

Figure 8 shows velocity deficit profiles in three different regions from vertical measurements on a single turbine. The velocity deficit is plotted as ΔU = 1 − U/U_∞ which indicates a lower measured velocity has a higher deficit. The vertical axis has been normalized by the blade radius with y = 0 equal to the hub height and rotational centerline. For the upstream measurements in Fig. 8(a), the deficit increases as the flow approaches the rotor and the deficit is symmetric about the rotational axis (y/r = 0). As expected, the highest deficit occurs closest to the hub (plane S7) as the flow stagnates near the center of the hub. With the outer edge of the hub at y/r ∼ 0.1, the spread in the upstream deficit between y/r = ±0.5 suggests the blades are creating a substantial pressure gradient in the radial direction. The pressure gradient is specific to the blade design and results from the larger blade width and blending near the root (i.e., higher local solidity). Downstream and within two rotor diameters, Fig. 8(b) shows a nonsymmetric deficit about the rotational axis; the deficit pattern below the hub height is more erratic and the deficit magnitude is higher due to turbulence generated from the tower and interaction with the rotating blades.

Above the hub centerline, there are three distinct regions created by shear layers and pressure gradients in the wake. Above y/r ∼ 1, the deficit is near zero as the velocity is equal to the free stream velocity. A shear layer is present at the outer edge of the blade near y/r = 1 where the wake and free stream interact. Between y/r ∼ 0.3 and y/r ∼ 1, the deficit increases for lower y/r values due to momentum extracted from the flow and turbulence created by the turbine blades. Figure 8 also shows the largest deficit occurs at y/r = 0 and between y/r ± 0.3; the large deficit is likely due to the combination of turbulence created by the blade, pressure gradient from the solid hub, and pressure gradient from the blades. As the turbulence dissipates and the free stream flow mixes with the wake flow, pressure and velocity gradients erode, and the deficit profile becomes Gaussian as shown in Fig. 8(c). The wake does not appear to expand as expected which may be due to lack of turbulent energy in the freestream flow and insufficient mixing in the shear layer. Note that there is still an additional deficit created due to the presence of the tower far downstream. Figure 9 shows there is very little change in the deficit beyond three rotor diameters where the majority of the turbulence and pressure gradients have dissipated. Figure 10 shows the corresponding velocity deficit profiles for lateral (horizontal) measurements at the same plane locations. Comparison to Fig. 8 shows similarity in the shape and magnitude of the deficit for z/r > 0 and indicates the flow would be axisymmetric without the tower. As with measurements in the vertical direction, there are three distinct regions where the shape is different; however, the gradient at the edge of the blade appears to be near z/r ∼ 0.8 which suggests the wake may have shifted toward z/r = 0 due to air flow around the traverse. The total turbulence intensity (%) at each measurement plane is shown in Fig. 11. Within the first two rotor diameters, the turbulence intensity (TI) is highest directly behind the turbine at planes, S8–S18, and at y/r values closer to zero (closer to the hub). For y/r < 0, there is significant turbulence generated by the tower and the interaction with the rotating blade.

As with the velocity deficit, there are three regions of interest for y/r > 0. At y/r = 1, there is a sharp gradient in TI due to the shear layer between the rotating wake and the free stream. While there appears to be a slight spike in the intensity at the S8 plane, a higher turbulence intensity was expected at this location. A helical vortex structure was seen in high-speed images and will be discussed in Sec. 3.4. Hot-wire anemometry is unable to capture the structure of the helical tip vortex, unlike PIV. Between y/r ∼ 0.3 and y/r = 1, the turbulence intensity is higher near the blades in planes S8–S11, but quickly dissipates to a nearly constant 5% value in the radial direction. Figure 12 shows there is little change in the TI beyond one rotor diameter in the downstream direction for this y/r range. Turbulence in this region is likely created by boundary layer separation and shear stresses at the surface of the turbine blades. Between y/r = 0 and y/r ∼ 0.3, the TI is the greatest closer to the turbine blades and dissipates more slowly in the freestream direction compared to y/r > 0.3. Turbulence in this region is dominated not only by pressure gradients, but also from global turbulence from surfaces of the hub, blades, and nacelle (DC motor in this case). From Fig. 12, the turbulence intensity does not dissipate until three rotor diameters downstream. There is a strong correlation in the profiles of the turbulence intensity and the velocity deficit where regions of higher turbulence intensity have a higher deficit; deficit recovery in the freestream x/D_b direction also corresponds to a reduction in turbulence intensity. Turbulence created by the tower and the interaction with the rotating turbine blades is significant and is of the same order or higher than turbulence generated by the blade itself. The turbulence generated from interaction with the tower does not appear to dissipate as quickly as the turbulence generated from the blade. The greater turbulence intensity and velocity deficit in this region (y/r < 0) suggests the tower geometry should not be neglected in wake studies and CFD simulations.

Hot-wire measurements, with two turbines spaced two blade diameters apart, were taken in 3.18 mm increments in both vertical (y) and horizontal (−z) directions at the planes shown in Fig. 13. The measurement planes have been identified as D0–D43 and dimensions for plane locations are listed in Table 1. Additional measurement planes were added after experiments with the single turbine to get a better mapping of the profiles along the length of the tunnel. Except at the planes immediately behind the rotors (D8–D12 and D21–D25), vertical measurements were taken at the turbine centerline 25.4 cm below the hub height and extending 25.4 cm above the hub height. The horizontal measurements were taken at the hub height extending laterally 25.4 cm from the centerline. Figure 14(a) shows the deficit profiles from vertical measurements between the first and second turbines and Fig. 14(b) shows the profiles after the second turbine, but only up to two rotor diameters downstream. Figure 14(c) shows the profile after two rotor diameters; if not for the tower, the profile would be Gaussian-shaped about the rotational centerline. Measurements upstream of the first turbine are not shown, but they were nearly identical to the values from experiments using a single turbine.

The general shape of the vertical profiles between the first and second turbines is very similar to the profile shapes after the second turbine. The profiles after the first turbine are nearly identical to the single turbine measurements which would be expected. As with testing with a single turbine, the velocity deficit profiles do not change appreciably after two rotor diameters downstream of the second turbine. There are distinct gradients in the profiles near y/r = 0.3 and y/r = 1, much like in the single turbine tests. The maximum deficits occurred at planes D13 and D26 which correspond to x/D_b ∼ 0.5. At x/D_b > 2, the deficit profiles show the wake is Gaussian-shaped after the second turbine. Looking at the deficit at several y/r values in the freestream direction, Fig. 15 shows how the second turbine adds to the deficit. There is a sharp gradient in the deficit after each turbine as momentum is extracted from the flow and pressure gradients are created across the turbines. Figure 16 shows the corresponding horizontal measurements with two inline turbines. In comparing to Fig. 14, the profiles indicate the flow is axisymmetric. Similar to experiments with a single turbine, the curves also indicate a shift in the wake with the gradient at zero deficit near z/r = 0.8.

The turbulent intensity for each region is shown in Fig. 17. The trends are the same as in the single turbine testing with the highest turbulence near the blade and very little change beyond three rotor diameters. The effect of the additional turbine can be seen clearly in Fig. 18 with a sharp rise in the turbulence intensity after the second turbine and a sharp decay within one rotor diameter. The sharp decay in intensity within one rotor diameter is a result of dissipation of higher frequency turbulence generated by the turbine blades and tip vortices.

A third turbine was added 40.6 cm downstream of the second turbine. Hot wire measurements, with the three turbines spaced two blade diameters apart, were taken in 3.18 mm increments in both vertical (y) and horizontal (−z) directions at the planes shown in Fig. 19. The measurement planes have been identified as T0–T48 and dimensions for plane locations are listed in Table 1.

Except at the planes immediately behind the rotors (T8–T12, T21–T25, and T34–T38), vertical measurements were taken at the turbine centerline 25.4 cm below the hub height and extending 25.4 cm above the hub height. The horizontal measurements were taken at the hub height extending laterally 25.4 cm from the centerline.

The velocity deficit profiles for vertical measurements are shown in Fig. 20 for three different regions with the probe positioned in line with the center of the turbine. Measurements upstream of the first turbine are not shown, but they were nearly identical to the values from experiments using a single turbine. The profiles between the first and second turbines are similar to the profiles between the second and third turbines. There are distinct gradients in the profiles near y/r = 0.3 and y/r = 1, much like in the single and dual turbine tests. The maximum deficits occurred at planes T13, T26, and T39 which correspond to x/D_b ∼ 0.5. At x/D_b > 2, the deficit profiles show the deficit has eroded and the shape is Gaussian-shaped.

Figure 21 shows the effect of additional turbines on the velocity deficit. After each turbine, the deficit increased sharply for all y/r values less than one due to momentum being extracted from the flow and pressure gradients across the rotor. Figure 22 shows the corresponding horizontal measurements with three inline turbines. In comparing to Fig. 20, the profiles indicate the flow is axisymmetric, but also show a shift in the wake with the gradient at zero deficit near z/r = 0.8. The turbulence intensity profiles shown similar trends as with the single and dual turbine tests, as shown in Fig. 23. The effect of multiple turbines is best seen in Figs. 24 and 25, where the velocity deficit and turbulence intensity have been plotted at normalized distances from the first turbine.

Since experiments with multiple turbines and hot-wire calibrations were performed on different dates, the graphs show excellent repeatability in the measurements as indicated by nearly identical values in the deficit and turbulence intensity values upstream of the first turbine, between the first and second turbine, and between the second and third turbine. Figure 24 also shows the velocity deficit along the shear layer y/r = 1 is unaffected by the number of turbines and the deficit at the centerline appears proportional to the number of turbines. For these experiments, turbulence generated from the tower was substantial and appeared to have a significant effect on the downstream deficit. The additional turbulence from the tower highlights the importance of including the tower in numerical studies.

In addition to hot wire measurements, high-speed filming was performed using a Photron Mini UX50 camera set to record at 2000 frames/sec using 1280 × 1024 pixel resolution. The camera had a 4 GB internal hard drive which limited the recording time to 1.09 s and 2180 total frames. A low-cost fog machine was used to introduce a nontoxic, high density, oil-based mist at the inlet of the wind tunnel. Multiple trials were attempted to direct the mist toward the center of the hub, but that did not always occur. While the high-speed imaging did not yield any quantitative results, it did provide qualitative information about the evolution of vortices in the wake and the turbulent structure in the wake. Figure 26 shows a clear vortex in the wake created by the rotor blades, and the vortex appears to break down approximately 2–3 rotor diameters downstream. The wake is obscured by the frame of the test section access door; however, the vortex dissipation appears to coincide with the x/D_b location where velocity deficit and turbulence intensity become nearly constant as shown in Figs. 9 and 12. There seems to be a trend that the wake shifts downward after the vortices have dissipated. The apparent shift could be due to gravity, phase change in the oil mist, or a general meandering.

For this study, commercial software, star ccm+, was used to simulate the model wind turbine described in Sec. 2. A direct model approach was used where the complete geometry of the rotor blades was discretized within a rotating subdomain. The Reynolds stress transport (RST), also called the Reynolds stress model (RSM), was applied as the turbulence closure method in solving the Reynolds-averaged Navier–Stokes (RANS) equations. The application of the RSM to a fully resolved wind turbine has been avoided due to the computational cost and potential for numerical instabilities. The advantage of using RSM is having access to shear stresses that are not available from two-equation RANS-based closure models like k–ε and k–ω. Access to the shear stresses will aide in the understanding of how the blade design will affect the wake, particularly in the near-wake region. All simulations were performed on the high-performance computing (HPC) cluster at the University of Wisconsin-Milwaukee. The HPC cluster includes 142 computing nodes with eight cores per node for a total of 1136 cores 3616 GB total RAM. For simulations in this study, only 32 cores were used. Batch processing was utilized to take advantage of the parallel computing architecture in star ccm+.

Experimental results, discussed in Sec. 3, indicated that the wind turbine rotational speed was dependent on the location of the vertical arm of a three-axis traverse system used to automate positioning of a hot wire probe. With the traverse arm located near the wind turbine and the wake, flow directed around the traverse arm affected the blade speed and wake velocity. Having the blade speed and wake-dependent on the traverse arm location creates a challenge for CFD simulations in that the location of the traverse has to be accounted for at each measurement position. In essence, the computational domain would have to be unique for each measurement point made in the vertical and lateral directions. As a compromise, multiple CFD models were created with a simulated traverse arm at specific downstream locations and only measurements in the vertical plane were used for comparison. Figures 27–30 show the computational domains for each of the simulations with the traverse arm located at specific positions.

The domains were chosen to simulate the turbine installed in a wind tunnel with dimensions of 1.21 m × 1.21 m × 2.44 m. The turbine blade is located 0.30 m downstream of the inlet. Each domain was separated into two subdomains to facilitate simulating the rotation of the blade using steady-state moving reference frame (MRF) and transient rigid body motion (RBM) approaches. The MRF model is a steady-state approximation in which cell regions move at different rotational or translational speeds. The flow in the moving region is solved using transformed continuity and momentum equations. It is important to note that the MRF does not change the position of the cell vertices and the rotor blades are frozen in space; forces are imposed on the cells that are induced by the rotation. The “frozen rotor” approach results in a solution that represents the time-averaged flow. At the interface between the regions, the local reference frame is transformed to allow flow variables in one region to be used to calculate fluxes at the boundary of adjacent regions. In the RBM method, the mesh vertices of the rotor region move to provide a time-accurate transient solution. The results from the steady-state MRF solution were used as the initial condition for the transient RBM solution. The rotating subdomain shown in Fig. 31 contains the grid elements of the blade and a Ø0.25 m region enclosing the blade. A new rotating reference frame was established for this subdomain with the axis of rotation set through the blade centerline and parallel to the 2.44 m tunnel length. The stationary subdomain included the tunnel, inlet and outlet sections, tunnel walls, tower, and nacelle. A mass flow boundary condition was applied at the inlet with the initial turbulence kinetic energy, k, and turbulence dissipation rate, ε, specified as 6.69 × 10⁻⁰⁴ J/kg and 1.0 × 10⁻⁰⁴ m²/s³, respectively.

Meshing in star ccm+ is mostly automated and based on the meshing model selected and cell size, cell growth, and refinement settings. For this research, the trimmer mesh model was chosen because the domain contains a substantial number of cells where the flow direction is aligned with the Cartesian coordinate system and because it produces a fast, high-quality mesh. A region-based meshing scheme was employed to coincide with the different physics models used for the rotating and stationary regions; however, additional constraints were added to the surfaces for boundary layer resolution and to improve the mesh quality. The rotor region consisted of predominantly hexahedral cells, but the automated meshing routine also created a mix of hexahedral, wedge, pyramid, and polyhedral cells based on localized prism meshing and mesh size definitions.

Six layers of orthogonal prismatic cells, as shown in Fig. 32, were added next to the turbine blade surfaces to ensure the boundary layer was well defined for laminar flow conditions and to ensure the wall y+ < 5 for turbulence modeling. Similar prism meshes were added to the surfaces of the tower, nacelle, tunnel walls, and traverse arm.

The base mesh size was initially set to 0.04 m and refinement through the domain was controlled by growth parameters. A rectangular grid refinement, as seen in Fig. 33, was added to reduce the element volume in the wake. The midplane mesh also shows downstream refinement as a result of the traverse arm. Mesh quality metrics was established for improved accuracy and convergence as follows:

Maximum cell skewness angle < 90 deg for all cells

Face validity > 0.95

Face quality > 0.2

Volume change > 1 × 10⁻¹⁰

For cells that did not meet minimum quality criteria, the cell quality remediation tool in star ccm+ was used. The tool identifies low-quality cells and their neighbors and modifies the gradients in those cells to improve solution robustness. The remediation is confined to the immediate vicinity of bad cells, so the influence on the solution accuracy is minimal.

All CFD studies are based on the fundamental equations of continuity, momentum, and energy. For this study, the fluid domain is considered isothermal and incompressible, so the energy equation is ignored. Defining fluid properties, such as velocity and pressure, as an average value with a perturbation, the Navier–Stokes equations can be transformed into the Reynolds-averaged Navier–Stokes (RANS) equations. The RANS equations include extra terms widely defined as Reynolds stresses.

A full 3D analytical solution to these nonlinear equations does not exist, so various methods of approximating the Reynolds stresses have been developed.

The main goal in RANS-based turbulence modeling is to develop a suitable closure method to predict the Reynolds stresses. Different types of turbulence closure methods have been developed, each with its own set physics, advantages, disadvantages, and applicability for certain types of flows.

Versteeg and Malalasekera [29] provide a good summary of the most popular methods. In this research, the Reynolds stress model (RSM) was selected for turbulence closure because it accounts for the anisotropy in the wake due to strong swirling motion and streamline curvature. Few studies exist using the RSM closure method for wind turbine analysis because it is a more advanced technique than one or two-equation turbulence models, often numerically unstable, and computationally expensive. Studies that do exist use an actuator disk to simulate the rotor. In addition to six equations needed to solve for the Reynolds stress components, an additional model equation is needed for the turbulence dissipation, ε. The RSM modeling strategy originated from Launder et al. [30] and variations of the transport equations have been developed over the years. star ccm+ offers a choice of four different Reynolds stress transport models: a linear pressure strain, quadratic pressure strain, two-layer linear pressure strain, and elliptic blending. For this research, the linear pressure strain was selected for stability and because it allows for a hybrid y+ wall treatment for coarse and fine meshes. The hybrid y+ wall treatment (referred to as an “All y+ Wall Treatment” in star ccm+) uses an exponential weighing function to blend calculation of turbulence quantities, such as dissipation and production, depending on the grid resolution at the wall boundary.

Both steady and unsteady incompressible flow simulations were performed on the direct-modeled wind turbine. The segregated solver of star ccm+ was used which solves the flow equations (one for each velocity component and one for pressure) in a second-order, uncouple fashion. The under-relaxation factors were set to 0.7 for velocity and 0.2 for pressure. For implicit unsteady simulations, the results of the steady simulations were used as initial conditions. A second-order temporal discretization scheme was used with a 2.083 × 10⁻⁰⁴ s time step and 20 inner iterations. The total simulation time was set at 0.375 s which represents ten revolutions of the rotor at the set rotation speed; each time step represented 2 deg of blade rotation. The convective Courant numbers were kept below 1 in the nonrotating region and less than three in the rotating region to maintain numerical stability and prevent the solution from diverging. For all simulations, values of the residuals, rotor thrust, rotor moment, inlet mass flow, outlet mass flow, and average blade surface pressures (direct model only) were monitored for convergence. The MRF and actuator disk simulations were allowed to run for 10,000 iterations. A grid sensitivity analysis was performed to ensure the grid size did not influence the simulation results. Grid size was reduced until the calculated thrust and torque changed by less than 4%. The final grid included 1,885,600 cells in the rotating subdomain and 1,837,728 cells in the stationary domain. For all grid resolutions, y⁺< 5 to ensure the boundary layer could be solved in the viscous sublayer.

From Fig. 7, the measurement locations have been designated as S8, S9…S16 for measurement locations where x/D_b < 1. Measurement locations where x/D_b > 1 have been designated as S17, S18,…S23. As mentioned, the location of the traverse affected the velocity measurement at each location. Initial simulations neglecting the effect of the traverse indicated errors of 50% compared to experimental data. To validate the application of the Reynolds stress model using experimental data, simulations were performed to compensate for the traverse location.

Three additional models were created with the traverse positioned to correspond to measurements at x/D_b ∼ 1, x/D_b ∼ 3, x/D_b ∼ 5. The locations were selected based on three regions having different turbulence characteristics. At x/D_b ∼ 1, the vortex structure created at the blade tips is strong, while, in the x/D_b ∼ 3 region, the vortices are mixing with the free stream the process of dissipating. At x/D_b ∼ 5, the tip vortices have completely dissipated. Comparative measurement planes are S17, S19, and S21. In the RBM simulations, vortices in the wake are captured due to the time-dependent nature of the analysis. Vorticity is a vector field and a measure of the local rotation of the fluid. Figures 34–36 show contour plots of the z-direction vorticity component with the traverse arm located at x/D_b ∼ 1, x/D_b ∼ 3, and x/D_b ∼ 5, respectively. In the contour plots, positive values indicate rotation in the same direction as the turbine blades and negative values indicate rotation in the opposite direction of the turbine blades. The contour plots show the helical vortex structure in the wake that results from vortices shed at the blade tips. The helical vortex rotates in the opposite direction as the turbine blades. The contour plots also show vortices shed from the traverse arm. As the traverse is moved farther away from the turbine, the turbine wake straightens and the wake behind the hub is extended as there is less mixing with the free stream flow. The wake behind the traverse appears to shorten, but this is likely a numerical effect due to the proximity to the outlet boundary. As a wind turbine extracts energy from the wind, there is a resulting wind speed decrease in the wake. The wake deficit is defined as the change in wind speed divided by the free stream wind speed, $ΔU/U∞$⁠. For this study, the different CFD methods of simulating the wind turbine are compared using velocity deficit profiles and experimental data. For the MRF simulations, the residuals and values of the thrust and torque tended to oscillate around a slightly meandering steady-state value which indicates that the influence of flow transients (i.e., separation and vortex shedding at the blade) may not allow the solution to converge completely in the steady-state solution. The velocity profiles also tended to change slightly even after 10,000 iterations. Figure 37 shows the effect of including the traverse arm in the simulations. With the traverse included in the simulation, flow diverted toward the wake created more mixing between free stream flow and wake flow which eroded the velocity deficit. The RSM turbulence model combined with the rigid body motion method shows excellent agreement with experimental data. Figure 38 shows the percent difference between the experimental data and the simulation results with the traverse location corresponding to measurements. Clearly, the RBM method is better at predicting the deficit with simulation results within 12% of the experimental results when accounting for the effects of the traverse. In addition to the strong agreement in the velocity deficit profiles using the Reynolds stress model with the transient RBM method, there is strong visual correlation between high-speed images and the numerical results. Figure 39 shows the isosurface of the Q-criterion scalar variable contoured by vorticity magnitude. The Q-criterion [31] represents regions where rotation dominates strain in the flow. The vortex compare well with the high-speed image shown in Fig. 26. The simulation showed the tip vortices dissipated approximately three rotor diameters downstream, which coincides with the high-speed images. Table 2 is a summary of the simulations performed for this study where the results of the MRF simulations were used as initial conditions for the RBM simulations.

This research has shown that the Reynolds stress model (RSM) for turbulence closure can be successfully applied to a fully discretized wind turbine and produce an accurate, numerically stable solution. In providing a solution to the RANS equations, RSM provides direct computation of the Reynolds shear stress components, as opposed to the eddy viscosity approach used in k–ε or k–ω turbulence models. Knowledge of the Reynolds stresses can help wind turbine designers and wind farm designers understand where regions of high turbulence kinetic energy occur in the wake by correlating to the degree of anisotropy of the Reynolds stresses. Transient effects, such as vortex shedding at the blade tips or from the tower, can be simulated and help mitigate noise or structural vibration concerns in the blades. Wind tunnel experiments have shown the interaction between the tower and the rotor can create significant turbulence that can be present in the far wake. Thus, the effects of the tower should not be ignored in CFD simulations.

This research was funded, in part, by grants from the National Science Foundation-Chemical, Bioengineering, Environmental, and Transport Systems (CBET) Division, the Department of Energy (DOE), and We Energies.

c = chord length (m)

D_b = turbine blade diameter (m)

H = hub height (m)

k = turbulence kinetic energy (J/kg)

r = turbine blade radius (m)

Re = Reynolds number, based on D_b

U = velocity component parallel to free stream (m/s)

U_∞ = free stream velocity (m/s)

x = downstream distance from the turbine blade (m)

y = distance to the nearest wall, vertical direction (m)

y⁺ = dimensionless wall distance

z = distance in horizontal direction (m)

Greek Symbols

ε = turbulent kinetic energy dissipation rate (m²/s³)

λ = tip speed ratio

ν = kinematic viscosity (m²/s)

ρ = density (kg/m³)

ω = specific turbulence dissipation rate (1/s)

Ω = angular rotational rate (rad/s)