3D NUMERICAL MODELLING OF LARGE SCALE IMPACTS OF TIDAL TURBINE ARRAYS USING AN OCEANOGRAPHIC MODEL

It is still challenging to predict the effects of large scale tidal turbine arrays on its surrounding hydrodynamic environment and the sediment transport process, especially when a realistic method has yet to be fully developed for the representation of tidal stream device in the existing oceanographic models. Generally, the commonly used regional oceanographic models are developed for near-horizontal flows, which make them inappropriate to simulate the complex 3D flows around the tidal turbine structure directly. Therefore, the present study aims to develop a threedimensional parameterization of a tidal turbine that can be used in a large scale oceanographic model, such as Finite-Volume, primitive equation Community Ocean Model (FVCOM). The additional retarding force method is extended in 3D flow conditions and applied in FVCOM to represent the tidal stream devices. Detailed laboratory measurements and computational fluid dynamics, CFD, calculated results are used to verify the model’s accuracy in prediction of hydrodynamics around the devices.


INTRODUCTION
The importance of getting started on utilizing renewable energy has been widely recognized in the face of climate change and reduction in natural energy sources.Tidal stream energy, which is the kinematic energy carried by ocean flows, can be harnessed in a way similar to that has been used in the wind energy industry for a while.Many researches have been carried out to illustrate ways a single turbine/a turbine array affects the near field hydrodynamic environment.For example, laboratory experiments with turbines simplified as porous disc were carried out by Myers & Bahaj (2007, 2010, 2012) and Maganga et al (2009) and more recent by Tedds et al (2014).Along with laboratory work, CFD modelling is another common way to study behaviours of a turbine.Sun et al (2008), Harrison et al (2009) and Bai et al (2009) approximated the turbine as a porous disc in their works using CFD software packages.More recently, works applying open source CFD codes with geometry of a turbine resolved in the calculating mesh are published (Bai et al 2014, Malki et al 2014, Goude et al 2014).In general, the impact is a chain effect that the flow rate behind the turbine will be reduced due to energy loss leading to shadow area where the turbulence eddy will be generated.The altered turbulence condition enhances the sediment pick-up in the shadow area which shows as bed scouring in the near field and traps more sediment into the water which will travel further downstream.The far field large-scale morphology is therefore very likely to be affected as more suspended sediments are available to settle down further downstream.
To parameterize the energy extraction process into a large scale oceanographic model to study the far field hydrodynamic or sediment transport changes, the most commonly applied scheme by far is to add a flow rate related sink into the momentum equations of the model.This concept was initially proposed by Bryden and Couch (2006) and then be adopted by Defne et al (2011) and Ahmadian et al (2012) with the additional sink term presented in different ways.But a coefficient deciding the strength of the sink term with different names in different presences is always involved.Defne et al (2011) and Ahmadian et al (2012) used a 2D model and some sensible results are produced.A 2D model requires only a depth-averaged coefficient without worrying about the vertical structure of the coefficient across the water depth.However, results predicted by a 2D model are largely limited to velocity and energy analysis.Shear stress which determines the behaviour of the particles carried by the water depends highly on the vertical velocity structure along the water depth.Unfortunately, the vertical velocity structure is missing in a 2D model, which makes it unreliable in predicting sediment transport patterns and, in the long run, seabed transformations.Consequently, a 3D model is used in this work to take in account the practical motions and transport.
It is noted that there are many different types of tidal stream energy devices currently being proposed and tested, including the horizontal axis tidal turbines (HATT) and vertical axis tidal turbines (VATT).They influence the passing flow in different ways.The swiping area of a HATT is perpendicular to the coming flow.It uses the lift force of the passing flow to generate the torque.A 'w' shaped x-direction velocity profile can be observed closely behind the turbine as the hub of the turbine slows down the flow dramatically and the flow rate increases gradually along the blade until the tip.A ring of increased flow rate compared to the inlet flow rate around the swiping area is expected.On the other hand, blades of a VATT are parallel to the turbine axis which is perpendicular to the seabed and are pushed by the flow to rotate around the vertical axis.Therefore the swiping area of a VATT is a cylinder which exerts an influence on the passing flow rate identically along its height.To parameterize the energy extraction process of a HATT three dimensionally, layers through which the blades swipe should be treated differently in a way that different coefficient should be assigned to different layers.While for a VATT, the coefficient can be the same along the blades.This work aims to parameterize HATTs.The turbine housing is not considered in this research.
An idealized channel is created in this work to apply the model.A single turbine is represented in the model as an additional retarding force working on the flow.Several tests are carried out to look into the relation among water flow rate, retarding force strength and energy coefficient.The model computed results are compared against CFD predictions to ensure the reliability of the parameterization method.

THE NUMERICAL MODEL FVCOM
FVCOM, which was originally developed by Chen et al. (2003), is based on 3D shallow water equations and discretises the domain using unstructured triangle mesh to fit the complex land boundaries in the horizontal plane.Along the vertical direction across the water depth, a terrain-following (sigma) coordinate transformation is used to fit with the bottom bathymetry and free surface.
In Cartesian coordinates, without considering of snow and ice, one of the momentum equations of FVCOM is: where u, v, and w are the x, y, z velocity components; ρ is the density; P a is the air pressure at sea surface; P H is the hydrostatic pressure; q is the non-hydrostatic pressure; f is the Coriolis parameter and K M is the vertical eddy viscosity coefficient.F u represents the horizontal momentum diffusion terms.

Parameterization method
The ocean dynamic is a complex system, involves tide, waves, turbulence, sediment transport as well as biological and ecological processes.An integral parameterization considering the impact of turbines is supposed to have all the physical process modules modified.However most of the existing studies concentrate on modification of the hydrodynamics only.There is generally a lack of 3D study and reliable parameterization method.Therefore this manuscript focuses on 3D current module parameterization with a preliminary analysis of the model reliability.
The concept proposed by Bryden and Couch ( 2006) is followed herein.The presence of the additional sink term in the momentum equations are that in Defne et al (2011) as: where F u and F v are the additional sink term components per unit area.C ext is the energy extraction coefficient determines the strength of the sink term.Across the water depth, different values for coefficient C ext are assigned to each vertical layer of the chosen triangles according to the turbine configuration and operation.

CFD Model setup
Facing the scarcity of on-site experimental data, a CFD model which is built from ANSYS FLUENT is used to provide validation data for the FVCOM results.The CFD model is firstly applied to a water flume with dimensions of 110m long, 42m wide and 50m deep.A HATT with a diameter of 15m is represented in the CFD model using Virtual Blade Method (VBM).The turbine is located at 10m downstream of the inlet and the hub of the turbine is set at the middle of the flume width and 20m above the bed leaving a distance of 12.5m between the lower diameter of the turbine and the flume bed.Multi-phase is activated to enable the water-air interaction at the water surface.The inlet flow is a plug flow with a flow rate of 1m/s.The rotation speed of the turbine in terms of rotation per minute (rpm) is 8.

FVCOM Model setup
To be comparable with the CFD model results, ideally, the dimension of the FVCOM model domain should be kept identical as the above open channel test in the CFD simulation.However, as FVCOM is based on shallow water equations which require the horizontal scale of the model to be much larger than the vertical scale, the water channel created with FVCOM is widened and lengthened into 1km in width and 10km in length with the water depth kept as 50m.Fig. 1 shows the mesh of the model, starts with a coarsest spatial resolution of 100m at the two ends symmetrically, and then it refines gradually towards the middle of the channel until it reaches the finest mesh size of 15m which is consistent with the turbine diameter in CFD model.In total 21 layers are used in the vertical direction across the water depth.The model is driven by surface water elevation time series.Both two ends of the channel are selected as open boundary.Through varying the water level difference of the two boundaries, the velocity of the water which flows from the left end to the right end can be changed.The bottom friction of the channel is reduced to a reasonable level to produce a plug flow at the inlet and maintain it across the whole length of the water channel.The turbulence module applied in this research is the Mellor and Yamada (1982) level 2.5 (MY-2.5)turbulent closure.

Energy extraction coefficient
The energy extraction coefficient needs to be selected before the model can be used for any predictions with turbine implemented.In Defne et al (2011), the model was optimized by changing energy extraction coefficient until the right amount of energy is extracted from the flow.The rated power is usually known for a given turbine design making it a reasonable indicator to judge the reliability of the coefficient set.It is understood that the sea condition changes which could affect the performance of the turbine in terms of power output.For example, it is mentioned in Ahmadian et al (2012) that the turbine is not always normal to the flow during a tidal cycle when the turbine is designed to rotate only 180° and hence the effectiveness of the turbine is reduced.But this process is not considered in this research mainly for simplification reason.It is an acceptable simplification as technology making sure the turbine is always facing the tide is available.Also, the velocity of the tide current surely changes over tidal cycles, which leads to the variation of the rotation speed and hence the energy output of the turbine.However the tip speed ratio (TSR) of the turbine is expected to be kept at its optimum regardless of the changing of the flow rate.In Allan et al ( 2013) the power coefficient (C p ) which is the ratio between the energy extracted by the turbine and the total available energy carried by the flow is found to be unaffected by the water velocity as long as a certain TSR is maintained.Thus, in the present study, a constant power coefficient (C p ) is assumed to lead to a constant set of energy extraction coefficient (C ext ), despite the varying flow rate.
Three tests with different inlet flow rate, 1.25m/s, 1m/s and 0.75m/s, are run with a fixed set of coefficient.One triangle at the right middle of the mesh is chosen to exert the coefficient set.The mesh size here (15m) is fine enough to represent a single turbine.The vertical structure of the coefficient set is shown in Fig. 2. Fig. 3 shows the normalized x-direction velocity component deficit (Δu) profile along the water depth of the three cases at the turbine location.It can be seen that profiles from the three cases almost merge together.Therefore, the power coefficient is the same for these three cases as energy is proportional to velocity.In this case, as a turbine is anticipated to work under its optimum TSR which decides a constant power coefficient for the turbine, a turbine design has its own corresponding set of energy extraction coefficient.Without considering the efficiency of the device, the power coefficient is 88.7% in this case, and the corresponding power produced is 208253.47,107901.68 and 45267.94Wfor cases with inlet flow rate of 1.25m/s, 1m/s and 0.75m/s, respectively.Another test is conducted by keeping the inlet velocity as 1m/s and varying the energy extraction coefficient set proportionally.Table 1 lists the energy extraction coefficient set scaling factors and the corresponding power coefficients and power outputs.The normalized x-direction velocity component deficit (Δu) profiles of these cases are drawn in Fig. 4. The profiles are in the same shape.But as expected, the velocity deficit becomes larger when the scaling factor increases.Curve shows the relation between the power coefficient and the scaling factor is drawn in Fig. 5. Following a rapid increase of the power coefficient up to the scaling factor of 0.2, the rate of the increase becomes steadier, and when the scaling factor falls in the range of 0.5-1.0, the power coefficient increases linearly but not proportionally.The linear correlation between the scaling factor (x f ) and the power coefficient, in this case, is: = 14.8  + 73.933 (0.5 ≤   ≤ 1) It should be noted, however, a high sensitivity of this equation is expected, which suggests that this result will not be applicable universally due to the limited validation and test conditions in the present study.

Model reliability analysis
In this study, the indicator chosen for the discussion of model reliability is the vertical structure of the x-direction velocity behind the turbine which depends highly on the vertical distribution of the energy extraction coefficient set.The CFD calculated result is used to compare with the FVCOM predicted velocity structure at several locations.The energy extraction coefficient set used in this case is the one shown in Fig. 2 and the inlet flow rate is 1m/s.CFD and FVCOM produced xdirection velocity profiles at locations of 4D, 5D and 6D are drawn against each other in Fig. 6.It can be seen that the FVCOM simulated wake is recovered much quicker than the CFD calculated one in the near field.However the turbulence behind the turbine is at the same level in these two cases in terms of turbulence intensity (TI) (see Fig. 7).But the vortex which has a length of about 2.5D behind the turbine caught by CFD is absent in the FVCOM case.This could be the reason resulting in the overestimation of the flow rate in the near field in FVCOM.The vortex itself also makes the results before 4D obtained through the two models incomparable.Fig. 6 shows clearly the way the parameterization method alters the flow structure on behalf of a real turbine.The velocity at layers with which the turbine intercepts is reduced due to energy loss as well as the blockage effect of the turbine.The hindered water finds its way to push further downstream through the neighboring layers, which results in the slight speed up of the water of these layers.The velocity of the surface layers is also decreased due to the friction between the water and the air.As the water flume created with CFD is only 110m long yet the diameter of turbine applied is 15m, it is unable to provide data further than 6D downstream the turbine.To recognize until where downstream the turbine the wake is dominant in the FVCOM predicted flow, two normalized depth-averaged wake velocity profiles are drawn in Fig. 8 based on data obtained through case with 1m/s inlet flow rate and original power extraction coefficient set.One is the velocity averaged across the whole water depth and the other one averaged at layers the turbine intercepts with.Location 0D is where the turbine is and location -1D represents the background state.The wake velocity profile averaged at layers strongly involve the wake shows a slower flow rate along the wake length in general.The recover process the turbine layer focused profile indicates is more dramatic at the first several diameters downstream.However the difference between these two profiles becomes less noticeable after a 5D distance behind the turbine and almost undistinguishable after a 10D distance.The wake recovery has reached 90% at the location of 5D as both profiles indicate.And the wake recovery reaches 95% eventually by 25D downstream the turbine.Malki et al (2014) reported a 92% recovery after 40D downstream the turbine based on CFD modeling, which could be quicker as the auther pointed out if the water tank was less smoother.
The unsatisfying wake recovery results recommend the need of looking into the mixing behind the turbine.The wake behind the turbine could be largely affected by the mixing.While the original turbulence closure module in FVCOM is not able to produce the desired wake.Therefore, in order to simulate the velocity reduction in the wake with better agreement, the turbulence closure module should be modified to take the presence of turbine and its operation into account.The incapability of FVCOM in catching the vortex just behind the turbine indicates that FVCOM may not be reliable when it comes to predict the scouring effect brought on by the turbine.This makes the modification of the turbulence closure module more justifiable.Also, the swirling effect which may carry the sediment for a longer distance and hence result in different downstream morphology alternation is not considered in the model yet.

DISCUSSION AND CONCLUSIONS
In the current study, a retarding force method has been extended to 3D flows and applied to the FVCOM model to examine the impacts of tidal stream devices on the flow hydrodynamics.Comparing with available validation data, the results shows good agreement which suggests the current method for representing tidal turbine is realistic (see Fig. 3).However, further tests involve the turbulence calculating scheme and unsteady flows need to be carried out to improve the method's reliability and applicability.

Figure 1 .
Figure 1.Triangular grid mesh of the water channel in FVCOM simulation.

Figure 2 .
Figure 2. Distribution of Cext across the water depth.

Figure 3 .
Figure 3. Distribution of computed flow velocity reduction rate across the water depth.

Figure 4 .Figure 5 .
Figure 4. Distribution of computed velocity reduction across the depth with different Cext coefficients.

Figure 6 .
Figure 6.Comparison of computed x-velocity from FVCOM with that from CFD.

Figure 7 .
Figure 7.Comparison of computed TKE from FVCOM with that from CFD.

Figure 8 .
Figure 8. Computed flow velocity distribution behind the turbine at the device level and averaged over the depth from FVCOM.