## Abstract

The modelling of tidal turbines and the hydrodynamic effects of tidal power extraction represents a relatively new challenge in the field of computational fluid dynamics. Many different methods of defining flow and boundary conditions have been postulated and examined to determine how accurately they replicate the many parameters associated with tidal power extraction. This paper outlines the results of numerical modelling analysis carried out to investigate different methods of defining the inflow velocity boundary condition. This work is part of a wider research programme investigating flow effects in tidal turbine arrays. Results of this numerical analysis were benchmarked against previous experimental work conducted at the University of Southampton Chilworth hydraulics laboratory. Results show significant differences between certain methods of defining inflow velocities. However, certain methods do show good correlation with experimental results. This correlation would appear to justify the use of these velocity inflow definition methods in future numerical modelling of the far-field flow effects of tidal turbine arrays.

## 1. Introduction

The development of tidal turbine arrays is gradually moving towards the development of full-scale demonstration farms [1] and fully commercial arrays [2]. Critical to the successful installation and operation of these is the optimal positioning of turbines, both within the farm or array itself and also in its respective channel. This will ensure both cost effectiveness as well as minimal environmental impact, which are likely to be critical to the success of any tidal farm or array project. Much research has examined flow effects in marine current turbines, with the aim of contributing to the understanding of how they affect the flow around them and perform under certain conditions. Garrett & Cummins [3] examined issues such as the constraining of flow causing localized increase in velocity around turbines. This work demonstrated that power output and overall efficiency were functions of flow velocity in different regions, and hence were dependent on array size and proximity of arrays to flow boundaries. This analysis was expanded further by Vennell [4], who looked at arranging tidal turbine arrays such that the optimal flow velocities were achieved in these important regions. This was done for different conditions including fixed numbers of rows of turbines and fixed cross-sectional area blocked by turbines. The sensitivity of power output to flow velocities in different regions demonstrated by these studies shows local flow magnitude increases and wake effects to be of critical importance in tidal farm development.

Experimentation with tidal turbine rotors is the ideal method of investigating the many parameters associated with their operation. However, the high costs involved make this method unacceptable for smaller studies with relatively modest resources. Validated numerical modelling therefore plays a crucial role in researching tidal turbine arrays. However, the many possible variables associated with modern numerical modelling techniques means that a thorough investigation is required to ensure modelled results represent real-life effects as accurately as possible.

This paper outlines an investigation into the numerical modelling of actuator fences. These have been used extensively in experimental work and numerical simulations to induce the momentum changes in a flow which occur with tidal and wind turbines. Examples include an investigation of tidal farm layouts (figure 1) in the study of Myers & Bahaj [5] and an analysis of wind turbine wakes in the study of Sforza *et al.* [6]. Experiments were previously carried out in the Chilworth hydraulics laboratory at the University of Southampton to examine wake and local flow constraining effects associated with different positions of tidal turbine arrays in an open channel [7]. Replicating experimental conditions using Ansys CFX software required the construction of a mesh of the geometry of the experimental flume using the ICEM mesh generator. Subsequently, a number of boundary conditions needed to be specified for different regions of the mesh, including no-slip walls for the channel sidewalls of the flume and a free-slip wall to represent the free surface of the flow. The inlet of the flume, where water is pumped into the channel from the sumps, also needed to be specified, along with some mathematical expression to represent the flow velocity magnitude throughout the depth of the flow at the inlet. A number of different mathematical expressions were used to define this flow velocity magnitude. This was to examine how the development of the flow changed with each expression and to examine which method most closely replicated experimental conditions. Determining which method for inflow velocity definition is most accurate can justify its continued use in future numerical modelling.

## 2. Literature review

The use of momentum sink methods to simulate tidal turbines in numerical modelling is a concept that has greatly simplified analysis of flow effects. Extensive work was carried out using these techniques in the study of Harrison *et al.* [8]. The aim of the study was to compare the characteristics of the wakes of actuator discs measured experimentally with those found using computational fluid dynamics (CFD) analysis and the Reynolds-averaged Navier–Stokes (RANS) equations. The experimental method involved measuring flow velocity magnitudes in wakes of actuator discs in a circulating water channel. The numerical method involved using proposed theory to calculate roughness coefficient and momentum sink for each actuator disc and applying this to the regions of the flow in CFD in which the actuator disc was present. This momentum sink method gave very similar far-wake properties in the numerical model to those experienced in the flume during experiments. Here, the far wake is generally accepted to be the region five actuator disc diameters or more downstream of the disc. For the actuator fence used in the analysis presented in this paper, this diameter is represented by the actual height of the porous rectangular fence.

The *k*–*ω* shear stress transport (SST) turbulence model was used in the simulation. The justification given by the authors was that it performed better in applications where adverse pressure gradients were present, and was suited to simulating a water flume due to the development of boundary layers on the flume walls. The findings of this analysis showed that velocity and turbulence parameters from CFD agreed well with those of experimental results. Also as expected, the near-wake characteristics did not agree well between experiments and the numerical model; however, in the far wake both methods correlated quite well. This is expected, as momentum sink methods such as actuator discs and fences are not specifically designed for the near-wake region. Device-generated turbulence and the initial velocity reduction in the near wake of tidal turbines are not equivalent between actuator discs/fences and mechanical rotors. However, a crucial factor mentioned in Harrison *et al.* [8] which justifies this method is that the downstream spacing between tidal turbines in a farm is likely to be seven turbine diameters or more, as with many wind turbine farms to date.

Momentum sinks were also used in Ansys CFX software in Lartiga & Crawford [9]. This analysis aimed to examine methods of correcting for the wall blockage effects in wind and water tunnels. These corrections are needed as blockage effects lead to inaccurate non-dimensional turbine performance parameters, such as power coefficient (*C*_{P}) and tip speed ratio (*λ*). The specific analytical correction methods that were analysed here were derived in Bahaj *et al.* [10]. The CFD analysis involved modelling both an unbounded domain with no blockage effects, and a bounded water tunnel with various levels of blockage. Results showed that the analytical expressions in the study of Bahaj *et al.* [10] agreed well with CFD for blockage of up to 10 per cent, but that theory underestimated *C*_{P} for higher blockage ratios.

The CFD software Ansys FLUENT was used in the study of Bai *et al.* [11] to investigate different spacing of turbines in a tidal farm. A number of different turbines were added to a farm, and the resulting effects on the central tidal turbine were examined. In the second instance, the number of turbines in the farm was fixed at seven, but with different latitudinal and longitudinal spacing used to examine the performance of the entire array. The computational domain used was almost identical to the basic setup of a conventional water flume, with sidewalls and an inflow, outlet and seabed. As with the work presented in Harrison *et al.* [8] and Lartiga & Crawford [9], horizontal axis turbines were modelled using a momentum sink approximation. Once again the *k*–*ω* SST turbulence model was used, with the authors justifying use of this method by stating that it performs better in regions of adverse pressure gradients.

For the first investigated scenario of adding turbines to the farm, results showed that the addition of turbines directly beside and downstream of the central turbine resulted in increases in the central turbine power output of up to 7 per cent. However, adding turbines upstream resulted in slightly decreased power produced. For the second scenario of varying the spacing between a fixed number of turbines, results showed that spacing increased total array power output by up to 1.3 per cent, but could also lead to decreases of up to 8 per cent. However, as the numerical work was not validated against experimental work, it is difficult to draw any quantitative conclusions.

These studies clearly showed many similarities in certain aspects of their numerical modelling approach. All three used the momentum sink approach to simulate the changes to the flow environment caused by an energy-extracting turbine. Both Harrison *et al.* [8] and Bai *et al.* [11] used the same *k*–*ω* SST turbulence model, citing its ability to cope better with near-wall results than other turbulence models. Lartiga & Crawford [9] also used a *k*–*ω* model, but it is not clear whether this was the same model as used by the authors of the other two studies. However, all three contrasted differed in the method in which they defined their inflow velocity. Both Harrison *et al.* [8] and Bai *et al.* [11] used different boundary laws, whereby the velocity at any point from the seabed or flume bed up to the water surface is some function of the height above the bed. Lartiga & Crawford [9] simply gave a single value for inflow velocity of 2 m s^{−1}. None of the authors gave any discussion of what other possible methods could be used, or any reason for using their particular method over others available. For any study that looks at the hydrodynamic effects of tidal turbines, determining the appropriate inflow velocity definition method could be crucial. It would ensure that the use of these methods in future work to expand on knowledge already gained would be strongly justified.

## 3. Experimental work and numerical modelling

### (a) Previous experimental studies

Previous experimental work in the Chilworth hydraulics laboratory is outlined in detail in Daly *et al.* [7]. A 21 m long and 1.37 m wide circulating water flume was used to carry out the experimental work. Initially, an ambient flow map of the flume was taken using an acoustic Doppler velocimeter (ADV). Flow velocity measurements were taken at several positions downstream of the flume inlet and at lateral positions of 0.27, 0.54, 0.84 and 1.11 m from the sidewall. Actuator fences were then used to represent tidal turbine arrays in a circulating water flume. Wake velocity magnitude measurements were taken downstream of fences. Fences were placed in different lateral positions 12 m downstream from the inflow to the flume to determine how changing proximity to flow boundaries affects the structure of the wake of the array. The 12 m point was chosen, as examination of the ambient flow map of the flume showed the flow had fully developed at this longitudinal location.

Results showed that the wake reduced in intensity and became more laterally extended as the actuator fence was moved closer to the sidewall of the flume. A separate study was carried out to examine the flow field between the actuator fence and sidewall owing to this change in proximity. Three fences of varying porosity (ratio of open to total area of fence) were gradually moved from the centre of the flume towards the sidewall. At each of these positions of the array, the flow velocity magnitude at a fixed point 0.5 fence diameters out from the sidewall of the flume was measured, and compared with the velocity at the same position when no fence was present in the flume. Results of this analysis showed an increase in flow velocity of up to 15 per cent occurred by decreasing the proximity of the array to the channel sidewall compared with that existing upstream. The aim of the numerical modelling outlined here was to validate these experimental results with a view to expanding on them in future research.

In order to compare experimental results to those of a numerical RANS solver, ICEM CFD software was used to recreate the geometry of the Chilworth flume. ICEM CFD is geometry creation software compatible with a number of numerical solvers. In this case, the Ansys CFX solver was used to replicate experimental conditions.

### (b) Geometry creation and meshing

After creating the basic geometry of the Chilworth flume in ICEM CFD, a region of the flow domain, given the name ‘fence’, was also specified as its own part. This region represented the region where the porous actuator fence used in experiments was present. The effect of the actuator fence on the flow was replicated by applying an additional momentum source term to this region, the value of which is directly related to the porosity of the actuator fence [8]. This was included in the RANS equations by Ansys CFX. Isolating this momentum sink region was done using the block splitting function of ICEM. A hexahedral mesh was applied to the geometry, with a finer mesh around more critical areas such as the close to the walls of the flume and around the actuator fence. The momentum source term *S* to be applied to this region is given by the following expression
3.1
Here *U*_{d} represents the flow velocity at the face of the actuator fence and *K* is the drag coefficient. An approximation for *K* proposed by Taylor [12] and used in Harrison *et al.* [8] is given by the expression
3.2
The resistance loss coefficient of the fence is then given by *K*/*Δx*_{i} where *Δx*_{i} is the thickness of the fence, which in this case is 4 mm. The experimental conditions for a 0.38 porosity fence were represented in CFX. Using equations (3.1) and (3.2), this corresponds to a value for resistance loss coefficient of 1481. An example of the ICEM geometry and mesh is displayed in figures 2 and 3.

### (c) Inflow velocity boundary condition

CFX-Pre was the software used for specifying the boundary conditions of the problem. In specifying the inflow boundary conditions to the computational domain, it can be useful to have some actual experimental data that can be used to define flow velocity and turbulence parameters. The plane along which actuator fences were placed during the experiments in Daly *et al.* [7] was 12 m from the inflow of the flume. A high-resolution flow field map was also taken to determine the ambient flow conditions in the flume in the absence of any obstructions. The closest downstream position to this at which data were taken from was 9 m downstream. This 9 m position was therefore taken in CFX analysis as the inflow point to the computational domain.

The following sections outline the different methods for defining inflow velocity magnitudes that were used. The curves of the expressions for the following three methods are shown in figure 4, with the experimental readings from 9 m downstream of the inlet of the flume and 0.54 m from the sidewall of the flume shown for reference. A comparison of the ambient flow predicted by both the experiments and all the methods is shown in figure 5, with the location being 12 m from the inlet of the flume and 0.54 m from the flume sidewall.

#### (i) Power law

A commonly used method of defining inflow velocities, also used in the study of Bai *et al.* [11], is to use a power-law profile. This is originally derived from boundary layer theory. This power law assumes the no-slip condition, which states that there is some region of flow towards the bottom of the channel where flow velocity magnitude equates to zero. The law further assumes that throughout the remainder of the hydraulic depth *D*_{h}, the velocity *U* at any point is a function of the distance from the bottom of the channel (*Z*). The expression for this power law is given by
3.3
Here, *U*_{max} is the maximum flow magnitude throughout the depth and *n* is any arbitrary constant. Although it is possible to use any value for *n* to attempt to replicate flow conditions, a comparison of several values for *n* was considered to be beyond the scope of this study. Therefore, one of the most commonly used values, 8, was inputted to the CFX simulation. Based on experimental results, 0.34 m s^{−1} was chosen as the value for *U*_{max}.

#### (ii) Exponential velocity profile expression

The next method tested was to set the inflow velocity conditions to the CFX simulation as close as possible to the boundary profile that was observed at the inflow point of the flume. The 9 m downstream from the flume inflow point was chosen as the inflow point to the CFX simulation. Data for flow velocity magnitude against vertical height from the bed were entered into Matlab software, and a curve was fitted to these data in an attempt to determine an expression which would give the closest value possible to the observed flow velocity magnitude for any input value of height above the bed. Several possible curve fittings were attempted, and it was found that fitting an exponential curve gave the most accurate results. The generic expression given by this curve-fitting exercise for any function *f*(*x*) whose value is dependent on the value of an independent variable *x* is
3.4
where *a*, *b*, *c* and *d* are arbitrary constants. In this case, *f*(*x*) will be the flow velocity, and the independent variable *x* is the vertical height above the flume bed. At the 9 m point from the inflow during experimentation, the change in flow velocity with height was recorded at four lateral positions along the width of the Chilworth flume. These were at distances from the sidewall of the flume of 0.27, 0.54, 0.84 and 1.11 m. The data recorded were entered into the curve-fitting exercise to determine values for *a*, *b*, *c* and *d*, as well as the value for certain statistical measures of accuracy. Table 1 gives the values returned for each data entry. Here, SSE represents the sum of squared errors, while *R*^{2} represents the coefficient of determination. This has been defined by Ryan *et al.* [13] as the fraction of the variation in the dependent variable that can be explained by the fitted equation. It can therefore be considered a measure of the accuracy of the equation in predicting observed real-life values, with a maximum of 1.

As table 1 clearly demonstrates, the expressions outputted by Matlab very accurately replicate the experimental conditions at the 9 m downstream position, and the differences between each are very small. Given these small differences and the time required, it was decided that using all four expressions for simulation was unnecessary. Given the small differences between them, any one of the lateral positions could have been chosen. In the end, the expression from a lateral position of 0.54 m was chosen.

#### (iii) Dyer boundary law

The third method to be tested was the Dyer boundary law. The theory behind this method is outlined in detail in Dyer [14], and it was also used in the analysis in Harrison *et al.* [8]. The Dyer boundary law states that velocity at any depth *y* (*U*(*y*)) is given as a function of the depth according to the following relationship:
3.5
where *U** represents the friction velocity, *A* is a constant and *ν* represents the kinematic viscosity of water. The values for *U** and *A* were found by curve-fitting experimental data to the above equation. This was done by curve-fitting data from the 9 m downstream point of the Chilworth flume, as was done with the curve fitting for the exponential velocity profile expression. This is done using the curve-fitting tool of Matlab, and using the custom equation creator to specify the above equation as that to which the data are fitted. This was done with all the lateral data found from the Chilworth flume. The values for *U** and *A*, as well as the sum of squared errors and *R*^{2} values found from the curve-fitting exercise, are displayed in table 2. To ensure consistency between all inflow velocity methods, the data for the 0.54 m lateral position fit were used in Ansys CFX.

#### (iv) Flow recirculation

The final method used differs quite substantially from the other methods used and those used in the studies of Harrison *et al.* [8], Lartiga & Crawford [9] and Bai *et al.* [11]. It involved using the interface capabilities of Ansys CFX to constantly recirculate flow from the outlet of the Chilworth flume back to the inflow, in order to allow the flow to continually develop. With the other methods which have been previously discussed, it was not possible to tell whether or not the flow had actually reached a point where it had fully developed by the time it had reached the outlet of the domain. It is possible that if a longer computational domain were constructed, the flow may have continued to develop within the numerically modelled domain. By allowing flow to recirculate from the outlet back to the inflow continually, it was possible to allow the flow to continually develop until it had reached a steady-state condition immediately upstream of the numerical domain.

Applying recirculation like this in Ansys CFX required some different boundary conditions from the other methods. Rather than a domain inflow and outlet, these areas were instead defined as two sides of an interface. Then a constant mass flow rate was defined across the interface. The flow then automatically developed throughout the depth and width of the flume according to other boundary conditions. Also unlike other methods, turbulence boundary conditions did not need to be specified, as the development of turbulence was automatic according to the roughness of the flume bed and the behaviour of the RANS equations. Once the flow had constantly been recirculated and found to have fully developed, a plane section of flow conditions was taken and then used as the inflow condition for the case where an actuator fence is present in the flume.

The mass flow rate through the flume was found from interpolation of the experimental data. The depth of flow was divided into the eight sections of the depth from which flow velocity measurements were taken. The mass flow rate for each section was then found using the expression for mass flow rate :
3.6
where *ρ* is the fluid density and *A* the cross-sectional area of the flume. The mass flow rate was then applied and the flow was allowed to fully develop. Once it had developed fully, the results for one plane of the flume were exported and used as the inflow condition for all meshes where an actuator fence was present.

### (d) Other boundary conditions

Following application of the correct conditions for the inflow velocity, the remaining boundary conditions were specified. The sidewalls of the flume and the bed were set as no-slip walls, while the flow surface, as it was not a physically constrained boundary, was set as a free-slip wall.

As all real walls are not hydraulically smooth, any information on how rough they might be is needed for CFD simulations so that turbulence can develop accurately. This is particularly important in the case of the recirculation method, as turbulence parameters cannot be specified here and so the development of turbulence is entirely driven by wall roughness. CFX requires an equivalent sand grain roughness to be specified for each wall, rather than the actual roughness itself. The relationship between sand grain roughness *k*_{s} and Manning’s coefficient *n* was originally developed by Nikuradse [15] and was expanded on subsequently by various authors. A final expression derived by Webber [16] and also presented in Marriott & Jayaratne [17] is
3.7
Manning’s coefficient for glass of 0.010 was entered into the above equation, and the equivalent sand grain roughness subsequently calculated entered into CFX.

As the method of inflow velocity definition, it was also possible to specify turbulence conditions for the inflow boundary condition. This was necessary for all the methods except the recirculation method for reasons already mentioned. Expressions were entered for the turbulent kinetic energy *k* and eddy dissipation *ε*. These quantities are functions of the flow velocity at any point and are given by the following expressions:
3.8
and
3.9
where *I* represents the turbulence intensity, *U*_{avg} represents the average freestream flow velocity and *D*_{h} represents the water depth, which in this case is 0.3 m. The turbulence intensity entered was given by the following expression:
3.10
where *Re*_{l} is the Reynolds number associated with the characteristic length. If the hydraulic diameter (which in this case is the water depth) of the flow is *L*, then the characteristic length is given by
3.11
The flow velocity used to compute the Reynolds number was 0.34 m s^{−1}, which, as with the 1/8th power law in numerical modelling, was chosen on the basis of observed experimental results. Finally, the *k*–*ω* SST turbulence model was selected. It was decided its good performance with adverse pressure gradients and separated flow meant it was ideal for the current flow problem.

### (e) Mesh independence study

Because of approximations and averaging used in the RANS equations, the number of cells in any given CFD mesh can directly affect the results after processing. A greater number of cells results in a more accurate solution to a problem for a given set of boundary conditions and solver settings. However, increasing mesh density also requires greater processing time. For any given problem, the optimal mesh gives sufficiently accurate results for certain important flow parameters in the shortest possible processing time. In this particular study, those important parameters are the flow velocity magnitudes throughout the computational domain. To verify the robustness of the mesh used, the flow velocity magnitude in the direction of flow at six different points was gathered for the scenario where an actuator fence was placed in the centre of the flume. The points chosen were at the mid-height and mid-lateral position 9 m from the inflow to the Chilworth flume (the inflow point to the CFX domain for reasons explained previously), 12.7, 15, 18, 21 m and the point 50 mm from the sidewall of the flume at which increases in flow velocity magnitude were observed and measured during experimentation. The 12.7 m position was examined as this lies in the immediate far wake of the actuator fence when it was in the flume centre, which was placed at 12 m from the inflow to the Chilworth flume. As replicating far-wake effects is what the actuator fence method is specifically designed for, ensuring robustness in this region of flow is of critical importance. The method used is a standard method recommended by the American Society of Mechanical Engineers (ASME), outlined in Celik *et al.* [18] and also used in Harrison *et al.* [8]. The authors describe a grid convergence method, which is based on the Richardson extrapolation method and has been used and tested in hundreds of CFD studies. The interested reader is referred to Celik *et al.* [18] for a more in-depth description of this method. The results are displayed in table 3.

## 4. Results and discussion

### (a) Mesh independence study results

The results of the mesh independence study are given in table 3. This was carried out for the actuator fence positioned in the centre of the flume. The inflow velocity definition method used for this was the 1/8th power law. The parameters across the top row are the *X*-direction flow velocity magnitudes in the centre of the flume at various distances downstream of the inflow. The last column records flow velocity magnitude at the point where flow localized acceleration was recorded in experiments. The parameters down the first column are some of the calculated values which [18] instructs to report. The coarse, medium and fine meshes had element numbers of 1 626 570, 3 107 982 and 6 339 168, respectively. In most of the parameters the errors are almost negligibly small. The largest errors occur at the point of the flume where velocity measurements were taken, and also at the 12.7 m point, which is in the far wake of the fence. The most probable reason for this point having higher deviation between meshes is that it is in an area where there are greater fluctuations in velocity magnitude between cells in all directions. Even these errors are not higher than 2 per cent. Overall, the results appear to support the view that the mesh with 3 107 982 cells gives sufficient accuracy to allow it to be used more extensively.

### (b) Constrained flow analysis

Figure 6 gives a plan view description of the experimental domain. An ADV was placed at a fixed point 50 mm from the sidewall of the flume and three actuator fence heights downstream of the fence. The fence was initially placed in the centre of the flume, and then gradually moved closer to the channel sidewall. For each of the positions, the flow velocity magnitude is measured by the ADV. This procedure was replicated in Ansys CFX, with the inflow velocity being defined using all the methods previously outlined. The results of each of these methods compared with experimental results are shown in figures 7 and 8. Figure 7 gives the flow velocity magnitudes recorded by the ADV at the point specified in figure 6. Figure 8 gives these ADV readings normalized against the natural flow velocity present at the same point with no actuator fence present.

Both figures show the same general trend. The flow velocity magnitude and increase compared with ambient conditions at the fixed point is highest when the actuator fence is closest to the sidewall of the flume, with a gradual reduction as the fence is moved further towards the centreline of the channel. This is an expected result and can simply be explained by the continuity equation. If the flow has less area through which to flow, the flow velocity needs to increase so as to ensure continuity of mass, with the effect becoming less pronounced as the area expands. However, both figures show quite contrasting performance of each numerical method in terms of both how they compare with measurements and the percentage increase in velocity due to constraining.

If figure 7 is examined independently, it can be seen that the recirculation method and 1/8th power law have results which replicate those of experiments most closely. Significant deviations from experimental results can be seen with both the Dyer boundary law and exponential curve fit. One reason postulated for the poor accuracy of these particular methods is that they may not replicate the natural boundary layer development off the bed of the Chilworth flume. Experimental work has shown that there is a gradual growth in the boundary layer from the inlet of the flume until the vertical velocity profile stabilizes at approximately 9–12 m downstream. This is analogous to the development of flow over a flat plate, as explained by Hamill [19]. In general, the boundary layer thickness settles to a constant value at some point downstream where the energy loss disruption to flow is replaced by momentum being transferred from faster moving upper layers of fluid to the slower lower layers. With the Dyer boundary law, Ryan *et al.* [13] claim that flow follows this profile in a region which begins just outside the viscous sublayer and finishes at 10–20% of the overall hydraulic depth. But Ryan *et al.* [13] also suggest that power-law profiles have been shown to be more appropriate for describing flow development in this outer region. Also crucial is that the 1/8th power law has also been extensively used in pipe flow, where the development of boundary layers is similar to that which occurs with flat plates. When considering the performance of the recirculation method, the no-slip condition is an important consideration. As this method does not strictly define flow velocity throughout the depth, the application of the no-slip boundary condition in CFX is one which will dictate flow development, along with wall roughness and turbulence parameters. This is therefore also more likely to replicate the boundary layer development than the exponential fit and Dyer boundary layer profile.

Figure 8 shows some interesting results for the prediction of increase in flow velocity compared with ambient conditions. When comparing results all numerical methods predict the measured localized increases in velocity owing to the constrained nature of the flow to within 5 per cent. In fact, the exponential curve fit, 1/8th power law and Dyer curve fit all have almost identical results to each other. This is an interesting result, as it shows that the diversion of flow around the actuator fence and the subsequent increase in local flow velocities is a phenomenon which is replicated equally by all methods. This diversion of flow appears to be handled quite differently by the recirculation method. Although in some instances it predicts magnitude of velocity increase close to the values of experiments, it does give much higher values at certain points. The fact that it does not have identical results to the other inflow methods is another observation which can most probably be attributed to the lack of strict rules on flow velocity through the depth.

The results of the above analysis would appear to suggest that the 1/8th power law is overall the method whose results correlate closest to those of experiments, and appear to suggest that it may be the most suitable for future numerical modelling of constrained flow. The next section will outline how all the inlet velocity definition methods compare when analysing the wake effects associated with tidal turbine arrays.

### (c) Wake modelling results

Figure 9 shows a plan view of the setup for the wake modelling of the actuator fence. The fence was initially placed in the centre of the flume, and the ADV was used to take flow velocity measurements at several points downstream. This procedure was then repeated for cases when the fence was moved closer to the sidewall. These cases were 200 mm, or two actuator fence heights from the centre, and 400 mm or four actuator fence diameters from the centre. Results of this wake mapping outlined in Daly *et al.* [7] showed a smaller reduction of velocity as the fence was moved closer to flow boundaries, and also the wake became more laterally stretched. Figures 10–12 show the velocity profiles downstream of the fence along a vertical plane intersecting the centre of the fence. Both experimental and numerical results (using CFX) are presented.

These results show that each numerical method displays the same general trend, with a quite dramatic decrease in velocity in the region directly behind the actuator fence. The low velocity magnitudes towards the mid-height of the wake are simply explained by the dramatic and sudden drop in momentum owing to the presence of the fence, while the higher values both below and above the fence are simply because of a gradual recovery towards freestream conditions owing to flow being diverted around the edges of the fence. When comparing the correlation of each velocity inflow method with experimental results, a very noticeable difference between each method is in the lower 10 per cent of the depth. The 1/8th power law fits quite well in this region, while other methods deviate quite substantially from experimental results. This once again supports the postulation made in the previous section that power laws perform better in regions well outside the viscous sublayer, further than the initial 10 per cent of the depth. It also supports the postulation that the Dyer law and exponential curve fit are incapable of replicating the natural boundary layer development in the Chilworth flume. As for the recirculation method, the deviation of this method from experimental results could be attributed to the relative lack of control of flow development associated with it. While the other methods quite strictly control the flow velocity magnitude throughout the depth, the recirculation method relies solely on the inputted mass flow rate, wall roughness and turbulence parameters. It is not clear how sensitive the results of the analysis are to each of these parameters, but getting values for each which are exactly the same as in the Chilworth flume may be quite difficult. As an example, the mass flow rate was calculated using interpolation from experimental results taken at a lateral distance of 0.54 m from the sidewall of the flume. It is most probable that the mass flow rate would have been different if experimental results were taken from a different position, leading to different wake results.

There are also quite considerable deviations in other regions of the depth, in particular, at depths furthest away from the centreline where simulations under-predict wake velocity deficits. One possible reason for this particular disagreement between modelled and experimental results may be a limitation of the use of momentum loss coefficients in CFX. The porosity of the actuator fence used was 0.38; however, there is no definite rule for how the hole pattern of the fence should be defined in order to achieve this. In operation, the highest induction factor and hence momentum drop are close to the blade tips for a conventional horizontal axis turbine. Modelling and simulating a membrane that has varying porosity might go some way towards achieving better replication of full-scale conditions and closer agreement between experimental models and numerical simulations.

However, it is not possible to account for this when applying the momentum loss coefficient calculated directly from the porosity in CFX using equations (3.1) and (3.2). Instead the momentum loss has to be applied either uniformly across the fence region, or by specifying individual losses in the *X*-, *Y* - and *Z*-directions. The uniform loss was chosen in this instance, but neither of these would be able to account for different fence hole patterns. Manually introducing hole patterns would undoubtedly lead to significantly greater computational time but could be a potential route to reducing the disparity between observed and simulated results.

The results of both the wake and constrained flow analysis appear to definitively rule out the continued use of both the Dyer boundary law and exponential curve fits. While they may be applicable to some particular instances, the above results strongly suggest that they are incapable of replicating natural flow development of the flume at the Chilworth hydraulics laboratory, and also do not perform as well as other suggested methods in replicating the wake of the actuator fence. The only method which appears to perform consistently well throughout the flow depth and as the fence is moved closer to flow boundaries is the 1/8th power law. The initial curve fit in figure 4 suggested that it may not be the most accurate method. However, its strong performance in the bottom 10 per cent of the flow depth is very significant, as it strongly suggests that it is able to replicate the natural development of the flow in the Chilworth flume, as already outlined by Hamill [19]. This could possibly be attributed to the fact that this method was originally developed from pipe flow, where the gradual transition of the flow velocity magnitude from no slip at the wall to freestream conditions is identical to that of a flat plate. The correlation of the results of this method with experimental results in analysing both wake and flow acceleration effects strongly supports its continued use in future numerical modelling of actuator fences.

## 5. Conclusions and future work

A study to investigate varying methods of defining the vertical velocity profiles at the inlet of a numerical model simulating a channel containing porous surfaces has been presented. Numerical modelling was carried out using Ansys CFX software and for each inlet definition the structure of the resultant flow field driven by the flume characteristics and the presence and position of actuator fences in the channel were compared with experimental results. Numerical results appear to show that the 1/8th power law predicts constrained flow effects towards the boundaries of the channel better than defining inflow velocities as following exponential and Dyer boundary law curve fits to experimental data. It also appears to have a superior performance to a method of defining a mass flow rate through the domain and allowing flow to recirculate until it has reached constant steady-state conditions. As well as constraining flow effects, the 1/8th power-law method is also seen to be most accurate at replicating the development of the wake behind a tidal turbine array. This appears to suggest that it is the method which is most capable of replicating the natural development of flow from the bed of the channel through to the freestream conditions. Future offshore flow field measurements around full-scale tidal turbines will require a study similar to that presented herein. There might be reasonable variation in the vertical velocity profiles between different tidal energy sites depending upon depth, flow speed and bathymetric profile. The inclusion of tidal turbines will again alter the flow field significantly and the inlet velocity profile of any numerical modelling will be driven by a myriad of factors including the number of turbines, the lateral blockage ratio and the proximity of devices to each other and also to solid boundaries.

This analysis has assumed a fixed size of actuator fence and hence tidal turbine array, and has examined the most suitable numerical modelling methods. Future analysis will use these methods to examine other sizes of fences and will examine the relationship between tidal array size and flow acceleration magnitude in the surrounding hydrodynamic environment, as well as the possible implications for the power output of other turbine arrays in the channel.

## Acknowledgements

The authors would like to thank the Engineering and Physical Sciences Research Council (EPSRC) for providing the funding for a PhD studentship.

## Footnotes

One contribution of 14 to a Theme Issue ‘New research in tidal current energy’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.