Abstract
This paper details a computational fluid dynamic (CFD) study of a constantly loaded actuator disc model featuring different boundary conditions; these boundary conditions were defined to represent a channel and a duct flow. The simulations were carried out using the commercially available CFD software ANSYS-CFX. The data produced were compared to the one-dimensional (1D) momentum equation as well as previous numerical and experimental studies featuring porous discs in a channel flow. The actuator disc was modelled as a momentum loss using a resistance coefficient related to the thrust coefficient (). The model showed good agreement with the 1D momentum theory in terms of the velocity and pressure profiles. Less agreement was demonstrated when compared to previous numerical and empirical data in terms of velocity and turbulence characteristics in the far field. These models predicted a far larger velocity deficit and a turbulence peak further downstream. This study therefore demonstrates the usefulness of the duct boundary condition (for computational ease) for representing open channel flow when simulating far field effects as well as the importance of turbulence definition at the inlet.
1. Introduction
The actuator disc method has been used together, with the Reynolds averaged Navier-Stokes (RANS) equations, for many years and for many applications including helicopter rotors [1], horizontal axis wind turbines [2], and horizontal axis tidal turbines [3, 4] alike. The actuator disc method represents a turbine as a simple disc of similar dimensions to the rotor and is used to approximate the forces applied to the flow. The forces are implemented as body loads or as negative momentum source terms on the flow as it passes through the disc.
The actuator disc approximation has a number of benefits over modelling the full rotor geometry. The most significant benefit amongst these is the reduction in computational expense especially for multiple rotor simulations. Full rotor simulations require a fine mesh to capture the boundary layer and separation along the blade surface, as well as the solution of the unsteady compressible Navier-Stokes equations. A full transient rotor simulation is needed, allowing the rotor blades to rotate in order to capture the wake. The actuator disc method allows for coarser meshes to be used and the incompressible Navier-Stokes equations to be solved as long as the Mach number is below 0.3, for this study the mach number is below 0.00021. Additionally, steady-state solutions can be obtained, vastly reducing the computational expense.
The actuator disc method has been a key tool of the renewable energy industry and has been used in a large number of studies [2, 5]. Even though there are more complex models such as the actuator line and full rotor models the low computational expense of the actuator disc method means it is still widely used [6, 7] and can be used to model multiple turbine interactions and wind farm simulations [8]. Although the actuator disc method has been used for many years, the majority of studies have used in-house code, as opposed to commercially available software, to conduct their studies.
The work described in this paper is essentially a benchmarking study, that is, a comparison of modelling data of previous theoretical, numerical, and experimental studies. Here the actuator disc method was implemented using the commercially available computational simulation suite ANSYS-Workbench with ANSYS-CFX v13 [9] and compared to the one-dimensional (1D) momentum theory as described in [10], as well as previous studies featuring both numerical [3, 4] and experimental data [11, 12].
2. Benchmarked Studies
Three previous studies were chosen for benchmarking featuring one theoretical [10], one numerical [4], and one experimental study [11]. The 1D momentum theory [10] also known as the simple actuator theory is an application of the 1D momentum equation applied to an idealized turbine. It uses control volume analysis to consider an infinitely thin frictionless disc with a constant momentum sink within an inviscid and incompressible fluid. The experimental data to be compared with is detailed in [11] and features three different porous discs to simulate different turbines. The experiment was conducted in a water channel measuring 21 m by 1.37 m with a depth of 0.3 m. Three 0.1 m diameter discs of various porosities were placed into the channel. The various porosities were used to represent different thrust coefficients (), which were measured using a pivot arm mounted onto a load cell. The water velocities were measured at various locations using an acoustic Doppler velocimeter (ADV) at a sample rate of 50 Hz and the data was averaged over 3 minutes. The previous numerical studies chosen for comparison are described in [3, 4] and used ANSYS-CFX to reproduce analogous experimental data [11]. These studies were chosen, because they were conducted using the same software as in the work presented in this paper, hence providing a benchmark to verify modelling methods.
3. Numerical Method
The numerical simulations used ANSYS-Workbench specifically ANSYS-CFX [9] and the steady-state solution of the Reynolds averaged Navier-Stokes (RANS) equations [13] together with the - SST turbulence model [14]. This model was chosen over the - model based on the literature and on some preliminary simulations which showed that the - SST model performs better in flows featuring adverse pressure gradients [15] in terms of accuracy to predict the flow properties. The - SST model was also used in the benchmark studies [3, 4]. The model domain was defined to the dimensions of the experimental channel setup [11]. It featured a 2 m long inlet, a 3 m outlet, and a 0.3 m deep-water column along with a 0.1 m diameter disc with a thickness of 0.001 m at the centre. The flow was assumed to be symmetrical, allowing a symmetry plane to be setup through the centre of the disc, dividing the domain in half creating a width of 0.685 m as opposed to the 1.37 m width of the experimental channel; this therefore reduces the computational expense. All simulations were carried out using water at 25 degrees centigrade corresponding to a density of kg/ and dynamic viscosity of kg/ms. As part of this study three discs were simulated with two different sets boundary conditions to represent a channel and a duct each with two different inlets totaling 12 simulations.
3.1. Boundary Conditions
The inlet velocity was defined in the same manner as the numerical study and based on the empirical data [11]. The equation used to define the inlet velocity was given in [4] as where is the inlet velocity across the width of the domain, is the friction velocity, is the depth of the water, is the kinematic viscosity, and is a constant. Curve fitting methods were used to define and . The numerical paper [4] used values of m/s and m/s were also used in this study. Figure 1 shows the inlet velocity used in this work and the experimental study [11], normalized with a free stream velocity of m/s for the experimental study [4] and m/s in this work. The vertical height was also normalized with the diameter of the disc .
The turbulence intensity, which is defined by (2), was described in two different ways to define two simulations referred to in this paper as inlet 1 and inlet 2. Both inlets were set with a turbulence intensity of at the inlet to produce agreement with the experimental data [11] for . The difference between the inlets is that inlet 2 was also defined with a length scale of 0.3 (height of the domain). Both these approaches are different to [4] which defined the turbulent kinetic energy and eddy dissipation
In (2)–(4) is the turbulence intensity, is the root-mean-square of the turbulent velocity fluctuations, is the mean velocity, and is the turbulent kinetic energy and is the velocity in the directions.
The outlet was defined as a static pressure outlet with a relative pressure of zero. The floor and far side of the domain were defined as a nonslip wall. In this study two separate models were produced featuring different boundary conditions at the top or roof of the domain; the first featuring an opening creating a channel and the second featuring a nonslip wall creating a duct. These boundary condition sets were analogues of those in [3, 4]. Although a free-surface approach may be considered more suitable, as the experiment was carried out in a channel featuring water and air interactions, it was shown to only produce a 0.2% depth change at the disc [4].
3.2. Disc Definition
The disc was defined with a diameter of 0.1 m and a thickness of 0.001 m as a subdomain with a uniform momentum loss across the disc in the longitudinal (-) direction. The momentum loss was defined using a directional loss model, which added a momentum source term () to the flow, which was defined as where is the resistance coefficient, is the density, and is the velocity. The resistance is applied as the loss across the disc thickness and so was specified by the user as , where is the thickness of the disc. The work described in this paper and previous numerical studies [3, 4] used identical resistance coefficient of 1, 2, and 2.5 in separate simulations to represent the three different porous discs used in the experimental study [11]. The resistance coefficient was derived in [3, 4] based on the thrust coefficient observed in the experimental data and was estimated using (6) which is a theoretical relationship between and [3, 4]. Here the value of is obtained from to render the required momentum deficit of 80%, 66%, and 61%, in this study:
3.2.1. Thrust Coefficient
The thrust coefficient () is a nondimensional variable used to describe rotor’s characteristics. The greater the value the greater the wake expansion and turbulence levels within the wake. The value of the porous discs in the experimental study [11] was measured using a pivot arm attached to a load cell. The thrust coefficient can be described numerically using (7). It requires the thrust () to be estimated which can be achieved in a number of ways. In [4] the thrust coefficient was estimated from the results using (8) to define the thrust. However in this work, as in [3], the thrust was calculated using (9):
In (7)–(9) is the thrust coefficient, is the thrust, is the density, the free stream velocity, is the disc area, is the resistance coefficient, is the velocity at the rotor, and is the change in pressure over the disc. The free stream velocity was calculated as the average velocity between at the inlet and was m/s in the experimental study [11], m/s in [4], and m/s in this work.
3.2.2. Mesh
In this work an unstructured hybrid mesh was constructed consisting of various mesh densities ranging from to cells. The majority of the domain was constructed out of tetrahedral cells with an inflated zone of wedge cells at the boundary of the floor and symmetry plane. The data presented in this work corresponds to a mesh density of approximately cells unless otherwise stated. Figure 2 shows the velocity profiles of various mesh densities along the centre line behind the disc and at 14 radii () downstream of the disc. There was very little difference between the predictions of the four different mesh densities showing little advantage in refining the mesh. Figure 2(a) demonstrates a realistic velocity recovery beyond the peak velocity drop just before . Figure 2(b) shows that the main differences between the different mesh densities are within the floor boundary layer and at the peak velocity deficit.
(a)
(b)
4. Results
All the data presented in this paper were produced using ANSYS-CFX and calculated with a root-mean-square residual of which was in line with [4]. The velocities which are compared to both previous numerical and experimental results were normalized using the free stream velocity of the flow described between at the inlet which was m/s in the experimental study [11], m/s in [4], and m/s in the work described in this paper. This range is being consistent with the empirical data. The inlet velocity profile, as shown in Figure 1, shows good agreement with the experimental data [11].
4.1. Influence of the Boundary Conditions
Before detailing the results a comparison of the effects of the boundary types is needed. To do this the same domain and mesh were set up excluding the momentum loss to observe how the velocity profiles develop without the influence of the discs. Figure 3 shows the differences between the channel (solid lines) and duct (dashed lines) velocity profiles as they develop through the domain. The figure shows how the channel flow is almost unchanged as the inlet was defined with a channel velocity profile. The duct profile changes significantly as expected with the additional wall boundary causing a sharp decrease in velocity at the top of the domain which forces the central velocity to increase to maintain the same mass flow rate.
4.2. Wake Predictions
The model in this work was compared to the 1D momentum theory as described in [10], specifically the pressure and velocity profiles along the centre line of the domain. Figure 4(a) shows the pressure and velocity profiles given by the 1D momentum equation. Figure 4(b) shows the pressure and velocity profiles produced by the model in this work. The overall profiles are in good agreement with only the magnitudes of the graphs changing depending on the value and the characteristics of the disc.
(a)
(b)
The experimental study [11] conducted experiments using discs with different porosity measurements to represent different values of . was measured in the experimental study [11] to be , and , respectively for each experiment. The simulations in the numerical paper [4] produced values of , and , respectively, for each disc using (8). The work presented in this paper calculated values of , , and for the channel simulation and , , and for the duct simulation, respectively, for each disc using (9).
Figure 5 shows the velocity profiles along the centre line of the domain and shows good agreement for all simulations in terms of the velocity characteristics, although the velocity magnitude is underpredicted compared to the numerical [4] and experimental data [11] for inlet 1 which has a delayed velocity recovery and appears to be offset from the other data sets. Inlet 2 shows a much better prediction of the experimental data [11] and both inlets show the duct has a quicker velocity recovery.
(a)
(b)
(c)
Figure 6 shows the turbulence intensity along the centre line of the domain and the difference between inlets 1 and 2 (Section 3.1) with inlet 1 having a lower starting turbulence intensity and subsequent peak. Both inlets show very little change in turbulent intensity just behind the rotor and then an almost linear increase up to the maximum intensity. This increase is most likely due to the presence of the wake edge shear layer with the maximum turbulence intensity indicating the merger of the layers and subsequent end of the near wake region. The near wake region of the flow which is defined behind the disc up until the wake edge shear layers meet at the centre line of the wake. The near wake region varies in distance generally from about to around downstream depending on the disc geometry and flow conditions. Figure 6 shows that the model is able to predict the intensity accurately far down stream of the disc, although it is unable to predict the peak in the turbulence intensity behind the rotor both in terms of magnitude and location.
(a)
(b)
(c)
Figures 7 and 8 show the velocity and turbulence intensity profiles at various distances downstream of the disc and for the three different discs. The distances downstream correspond to , , , , and downstream of the disc. These locations where chosen as they were the locations where the experimental data [11] was measured.
(a)
(b)
(c)
(a)
(b)
(c)
For all profiles in Figure 7 agreement was achieved (at least from a qualitative viewpoint) for the majority of the profile characteristics, such as the locations of highest and lowest velocities with the main numerical discrepancy at the maximum velocity deficit for all simulations. While the initial velocity drop is overpredicted at the centre, the free stream and floor boundary layer features are predicted well. Figure 7 shows that the duct simulations, displayed as the dashed line, predicted a smaller velocity deficit than the channel simulations at the centre. However, the duct model predicts higher velocity values towards the boundaries. Figure 7 shows quite well how the velocity deficit of the experimental data recovers quicker than the numerical data with inlet 1 simulations recovering the slowest. The experimental data seems to have almost recovered by and completely recovered by downstream whereas all numerical simulations still show some velocity deficit at .
Figure 8 shows the turbulence intensity of the models in comparison with the experimental data. There is little difference between the solid and dashed lines representing the channel and duct flows, respectively, and all models predicted intensities below that of the experiment data. The figure shows how the experimental data peaks earlier and higher than the modelled numerical data. Beyond approximately downstream of the rotor the modelled and experimental data are very close.
5. Discussion
Two types of simulations were carried out in the work described in this paper to represent a channel and duct flow. It was observed that the duct had a higher central velocity magnitude and marginal lower turbulence intensity than the channel flow with two different turbulent inlets. This is to be expected due to the presence of the additional wall at the top of the domain. The wall creates an additional boundary layer which restricts and slows the flow near the wall. This deceleration along the wall focuses the flow increasing the central velocity magnitude in order to maintain the same mass flow rate which is clearly visible in Figure 3. However, the influence in the region is minimal, meaning that representing the open channel flow as a duct incurs minor error, whilst reducing computational expense.
The difference between inlet 1 and inlet 2, through defining the turbulence length scales, had a significant effect on the simulation results. Inlet 2 predicted a more realistic velocity and turbulence intensity profile when compared to the experimental data [11]. This is due to a reduction of the turbulent dissipation throughout the domain prolonging the turbulence generated at the inlet and disc which was overly dissipated using inlet 1.
The models detailed in this work seem to have an inherent weakness in the definition of the momentum source as a predefined constant unidirectional loss. Allowing this, they performed well and predicted some characteristics of the velocity profile and turbulence levels. The model generally predicted both the velocity and turbulence intensity magnitudes lower than the experimental data [11] in the near wake region, with the discrepancy reducing as the flow moves downstream.
All the simulations carried out in this work predicted the turbulence intensity peak at far lower magnitude and further downstream of the disc than the experimental data [11]. This is due to the definition of the discs within the model as opposed to the physical discs. The discs within the experimental study [11] were porous discs with different values created through different porosities. These porous discs extracted momentum from the flow by converting the velocity into small scale turbulence and, thus, creating a high level of turbulence behind the disc. However the discs within the model extract momentum from the flow explicitly, reducing the velocity with no added turbulence. This explains the very high levels of turbulence behind the rotor for the experimental data and the lack of this peak in the modelled results. The turbulence intensity of the model peaked further downstream than observed experimentally (Figure 6); this is due to the merger of the boundary layers created by the velocity deficit. The porous discs within the experimental study [11] produced a variable 3D momentum loss. Although the discs in this work were defined with an isotropic 1D momentum loss, this reduced the amount of mixing and therefore produced a longer wake, implying the presence of anisotropic momentum losses.
The differences between the values calculated from the channel and duct flow can be attributed to the added boundary layer and subsequent small velocity increase. The values were calculated using a free stream velocity based on the velocity inlet which is perfectly reasonable for the channel flow as this velocity profile is fully developed. However, this is not the case for the duct flow. Taking this into account the values were recalculated, using a new free stream velocity of m/s obtained at the domain origin in the absence of the disc. Application of (9) produced new values of , , and for the three discs, respectively.
Use of the - SST model was based on the literature and Figure 4 implies that the adverse pressure gradient before the disc is well predicted as well as the floor effects shown in Figures 3 and 7. The - SST model is most appropriate in situations with adverse pressure gradient and 3D flow phenomena featuring strong swirl but the work described here only considered a 1D momentum source. The ability of the - SST model to prevent the overprediction of eddy viscosity may have inadvertently reduced turbulence in the wake and led to the longer wake seen in Figure 7 when compared with the experimental data [11] that had higher turbulence levels and 3D effect from the porous discs.
6. Conclusion
The work described in this paper has used the steady-state RANS solution method resident within the commercially available ANSYS-CFX [9] to benchmark an actuator disc model without rotation with the 1D momentum theory [10] and previous numerical [3, 4] and experimental studies of porous discs [11, 12]. This study has compared four different boundary condition sets, a channel and a duct, each with two different turbulent inlets containing different actuator discs. The discs were simulated featuring different resistance coefficients to represent different porous discs used within the experimental study [11]. These simulations were found to be in good agreement with the 1D momentum equation [10] in terms of the velocity and pressure profiles. When compared to experimental data [11], model predictions deteriorated with respect of velocity and turbulence intensity magnitude, just behind the disc, although the agreement improved further downstream. This discrepancy can be attributed to small scale turbulence present in the experiments and the momentum extraction method employed by the models.
This paper shows that the model method was sufficient to predict the far field velocity characteristics of a porous disc. Our future studies will model three-dimensional anisotropic effects at the disc by using variable momentum sources and include an additional turbulent source terms to account for the discrepancies found. Equally, more sophisticated modelling techniques such as adding rotation, the actuator line, surface model, or using a more sophisticated solver such as large-eddy simulation (LES) may produce closer agreement with field data. Moreover, such techniques might go some way to explaining the inherent poor prediction of turbulent intensity.
The main achievement of this study was demonstrating the usefulness of the duct boundary conditions (for computational ease) for representing an actuator disc in open channel flow when simulating far field effects, given the particular velocity profile, which is (1), applied at the inlet. It was found that the channel and duct simulations predicted very similar results with the duct predicting a slightly higher velocity magnitude for the majority of the domain. The main discrepancies observed when compared to the experimental study [11] can be attributed to the definition of the momentum source, which explicitly extracts momentum from the flow rather than converting it into small scale turbulence. Moreover, only a unidirectional momentum source was used, which did not account for the three-dimensional effects of the real discs used in the experimental study [11].
Nomenclature
| : | Boundary layer model constant |
| : | Area of the rotor |
| : | Thrust coefficient |
| : | Rotor diameter |
| : | Turbulence intensity |
| : | Resistance coefficient |
| : | Turbulence kinetic energy |
| : | Pressure |
| : | Thrust |
| : | Velocity |
| : | Free stream velocity |
| : | Velocity at the rotor |
| : | Velocity in th (=, , or ) direction |
| : | Friction velocity |
| : | Mean velocity |
| : | Inlet velocity |
| : | Root-mean-square of the turbulent velocity fluctuations |
| : | Distance in the widthwise direction |
| : | Distance in the depth-wise direction |
| : | Distance in the stream-wise direction |
| : | Water depth |
| : | Kinematic viscosity |
| : | Density |
| : | Disc thickness |
| : | Change in pressure over the disc. |
| CFD: | Computational fluid dynamics |
| RANS: | Reynolds averaged Navier-Stokes. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.