Computational Fluid Dynamics and Visualisation of Coastal Flows in Tidal Channels Supporting Ocean Energy Development

Flow characteristics in coastal regions are strongly influenced by the topography of the seabed and understanding the fluid dynamics is necessary before installation of tidal stream turbines (TST). In this paper, the bathymetry of a potential TST deployment site is used in the development of the a CFD (Computational Fluid Dynamics) model. The steady state k-and transient Large Eddy Simulation (LES) turbulence methods are employed and compared. The simulations are conducted with a fixed representation of the ocean surface, i.e., a rigid lid representation. In the vicinity of Horse Rock a study of the pressure difference shows that the small change in height of the water column is negligible, providing confidence in the simulation results. The stream surface method employed to visualise the results has important inherent characteristics that can enhance the visual perception of complex flow structures. The results of all cases are compared with the flow data transect gathered by an Acoustic Doppler Current Profiler (ADCP). It has been understood that the k-method can predict the flow pattern relatively well near the main features of the domain and the LES model has the ability to simulate some important flow patterns caused by the bathymetry.


Introduction
Modelling of the flow in an estuary or a channel depends on several parameters such as the bed friction factor and topography of the seabed. Other features like boulders and pinnacles can also affect the shape of the velocity profile and therefore affect the overall performance of a TST. Any model which is used to characterise the flow in a proposed TST deployment site has to have the ability to capture small scale turbulence. Turbulence effects on the structural integrity of turbine blades have been of interest for both wind and tidal energy [2][3][4]. In order to overcome this problem, ADCPs have been used to provide a better understanding of the flow properties [5][6][7].
Research has also been done to use ADCP data to generate turbulence synthetically to be used as an input to a Blade Element Momentum Theory model [8] to assess TST performance in very turbulent conditions [9]. In previous studies used to assess tidal energy resource, depth and velocity averaged data is extracted from admiralty navigational information such as [10][11][12][13]. These assumptions have been criticised for incorrectly predicting the available resources. There are several issues such as "turbulence quantities", "the rate at which a wake mixes with the free stream" and "velocity profile across the rotor" which must be considered [14]. A recent ADCP study which measured the properties of complicated flow in Ramsey Sound (West Wales) confirms that flow fields are spatially heterogeneous in three dimensions. Thus using depth-averaged velocity data to compute TST power extraction could provide inaccurate estimations [15], highlighting the necessity of detailed flow modelling for TST deployment sites.
Over the past few decades, large scale ocean modelling has been conducted and codes like TELEMAC [16,17] have been developed and are currently in practice to model free surface flow and waves in coastal regions. Research has been done in order to model the hydrology of ocean for the tidal sites [18,19]. All of this research is either based on using Shallow Water Equations (SWE) or applying Finite-Volume primitive equation Community Ocean Model (FVCOM). The SWE method is used when the ratio of depth to the length of the domain is very small and it is assumed that there is not any vertical shear on the flow. FVCOM uses the finite volume method so it is better at satisfying the mass conservation in comparison to methods which are using the finite element method (FEM). The FVCOM model consists of momentum, continuity, temperature, salinity, and density equations [20]. A 1D turbulence model which is called the General Ocean Turbulent Model (GOTM) has been added to the FVCOM model as well to solve the vertical mixing of the flow [21]. In all of the above-mentioned methods, an oceanographic grid was used to only model the general patterns of the flow and small scale turbulence (with a scale smaller than the grid cells) was neglected.
Selection of an appropriate modelling technique for a particular computational model is dependent on the requirements of the problem. Since the generation of eddies due to the seabed geometry is the area of interest for this research, turbulence has to be modelled. To model the small scale turbulence, oceanographic models are not suitable as they only model the general patterns of the flow and small scale turbulence (with a scale smaller than the grid cells) is neglected. It is necessary to employ turbulence models that can quantify the flow fluctuation in all directions as well as the important parameters which define the turbulence (turbulence kinetic energy, turbulence dissipation rate, turbulence length scale, etc.).
In general, it is possible to categorise turbulence models as either steady state or transient. Among the steady state frameworks, the two-equation turbulence methods (k-and k-ω) are the most commonly used [22]. For transient simulations, LES methods provide good results along with a lower computational cost in comparison to direct numerical simulation (DNS). Over the years a large amount of research has been done to model small scale sub-grid turbulence [23][24][25][26]. Since resolving turbulence is computationally costly, the majority of turbulence models, particularly the transient ones, have been applied to small scale channel test cases and their results have been validated against the laboratory tests or benchmarked models. Almost all of these test cases did not consider the effect of real bathymetry in their efforts for small scale turbulence modelling. This research will concentrate on modelling flow over a rigid surface with a consideration of the seabed geometry and without considering any vegetation and sediment transport.
In addition to the novel flow analysis, this paper describes the application of a state-of-the-art visualisation system so that a greater understanding of the flow is achieved. Stream surfaces [1,27] have important inherent characteristics that can enhance the visual perception of complex flow structures. Lighting and shading reinforce the perception of shape and depth, images or textures can be mapped to the surface primitives providing additional visual information, colour and transparency can be used to convey additional data attributes. Stream surfaces are able to capture the features within the flow, and also have the inherent ability to convey further information about the local attributes of the flow. This, combined with the reduction in visual clutter when compared to using glyphs or streamlines, significantly enhances their utility for practitioners.

Theory and Method
In this research the turbulence in the flow is modelled using the kmethod for steady state problems and the LES method for transient problems. The reason for choosing the kmethod is that it is proven to provide accurate results in most cases and this makes it by far the most popular two-equation turbulence model [22]. Large eddy simulation is used as it can model the instantaneous turbulence properties of the flow with good approximation without having the extensive computational costs of direst numerical simulations (DNS) and the near solid boundary simplifications of detached eddy simulation (DES) [22]. These methods will be discussed in detail in the following sections.

Standard k-Turbulence Model
One of the most commonly used RANS (Reynolds-averaged Navier-Stokes) turbulence models is the standard kturbulence model. The time averaged form of Navier-Stokes equations can be written as: where ρ is density, p is the pressure, v is the velocity vector, v i is the ith component of the velocity vector, µ l and µ t are the laminar and turbulent dynamic viscosities respectively, and S i includes any additional sources [28,29].
The main focus of the kmodel is on the mechanism of turbulent kinetic energy and can be expressed by two transport equations. The main parameters are k which is the turbulence kinetic energy and which is the turbulent dissipation rate: The values, C 1 , C 2 , σ k , and σ , are the closure constants, which are found empirically and are 1.44, 1.92, 1.0 and 1.3 respectively [29,30]. The turbulent generation rate G is defined by: The dynamic turbulent viscosity is calculated using [31]: The parameter C µ is a constant with the value of 0.09. The coupling between turbulence and momentum equations is achieved through the dynamic viscosity which is a sum of the laminar and turbulent values. In all of the equations, since the flow is incompressible, the derivative of density over time is zero [32]. The turbulence length scale gives a better understanding of turbulent flow. This is representative of the large scale turbulence and is defined as "a physical quantity describing the size of the large energy-containing eddies in a turbulent flow" [33]. The following equation is used to calculate this parameter which is used in this work to obtain the inlet boundary conditions for a CFD simulation:

Large Eddy Simulation (LES)
To solve the large eddies explicitly a filtered approach to the Navier-Stokes equation is used in which the filter is defined by the grid size. The small eddies are modeled implicitly by using a subgrid-scale model (SGS model) [34].
The filtered continuity and momentum equations use filtered variables and in which τ ij is the filtered stress tensor and σ ij is the subgrid-scale Reynolds stresses. Physically, the dynamic coupling between large and small scales in turbulence is described by the SGS stress. Dimensionally, it scales quadratically with turbulent velocity differences at scales of order ∆ [35].
To solve the SGS model one of the most commonly used approaches is the Smagorinsky [36] model.
where µ t is the turbulent viscosity, ∆ is the length scale of the filter, S ij is the filtered strain-rate tensor, C s is the Smagorinsky constant. The length scale of the filter, ∆, is usually considered to be the cube root of the volume of the cell [34]. Choosing the right inlet boundary condition is important for LES because it affects the quality of the results. In the steady-state models, a uniform inflow speed along with the turbulence parameters is defined whereas for the LES method a turbulent unsteady inlet condition is needed [37]. The Vortex Method has been developed for generating inlet velocity fluctuations [38]. This method, instead of applying random noise at the inlet, provides some spatial correlations for the fluctuations. In this method random vortices are created on the inlet flow plane (normal to the streamwise velocity) in a way that the wall-normal components gives a spatial correlation, and the Lagevin equation [39] is used to provide a temporal correlation not only for the generation of the streamwise component but also for the other two cross-stream components [34]. The vortex method is based on a Lagrangian formulation of Navier-Stokes. The centres of the vortices are transported and the velocity is given a certain distribution [34].
In the Vortex Method the vorticity is calculated by defining a circulation parameter, the intensity of which is quantified as a function of the turbulent kinetic energy that can be calculated from a RANS calculation [40]. For each vortex, a characteristic time (which is also known as the turbulent time scale = k/ ) is defined which represents the life time of a vortex. If the time exceeds this characteristic time, then it is destroyed and another random vortex is created on the 2D plane in the inlet. The diameter of the vortices can be set to a fixed number or can be calculated by using a turbulent length scale approach l = k 2 3 / . Then the vortices carry a random walk in the inlet plane to develop fluctuation in time [34].

Visualisation Technique
A stream surface [1] is the integration of a one dimensional curve through 3D steady flow. The resulting surface is everywhere tangential to the local flow. Since there is no normal component of the velocity along stream surfaces, they are useful for separating distinct regions of similar flow behaviour. These surfaces must represent an accurate approximation of the underlying simulation. Adequate sampling must be maintained while reducing the unnecessary computational overhead associated with over sampling. In practical applications a discretised approximation of the stream surface is constructed by integrating discretised seeding curves through the vector field.
To describe the stream surface [1], a streamline is first defined. A streamline is a curve which is tangential to the velocity field at every point along its length. Streamlines show the direction fluid flows within a steady state flow domain. If v( p) is a three dimensional vector field, the streamline through a point, p 0 , is the solution I s ( p 0 , t) to the differential equation: where, in the case of a steady state flow field, t is time, and the initial condition I( p 0 , 0) = p 0 . A stream surface [1] is the trace of a one dimensional seeding curve, C, through the flow. The resulting surface is everywhere tangential to the local flow. A stream surface, S, is defined by: where S is the union or continuum of integral curves passing through the seeding curve C. S(s, −) coincides with an individual integral curve, and S(−, t) coincides with individual time lines.
With the use of surfaces problems such as occlusion can arises. This may stem from multiple surfaces that occlude one another, a large surface that results in self occlusion, or a combination of both. The use of transparency is a generalised solution to this problem. With stream surfaces additional options are available, for example, illustrative techniques can by utilised to improve the visual cues, and thus the visual perception of the underlying data [27]. Figure 1 highlights these characteristics, and demonstrates the use of a visualisation framework [41] for prototyping visualisation algorithms. Figure 1. This illustration demonstrates a simple stream surfaces seeded in front of a cuboid, and capturing the flow characteristics, i.e., vortices, emanating downstream from the cuboid. This visualisation is constructed using the framework described in [41].

Domain
One suitable site in the UK for deployment of TSTs is Ramsey Sound in Pembrokeshire, Wales. The Sound is confined by the Irish Sea to the north and south and it is approximately 3 km long. The width of the channel varies between about 500 m to approximately 1600 m [42]. Kinetic energy of the water has been analysed with a regional coastal model and an estimate of the amount of power available in Ramsey Sound is 74.9 GWh/y [43]. It is going to be the site of DeltaStream, Wales' first tidal stream turbine, which has consent to install a 1.2 MW pre-commercial prototype device for 12 months [44,45].
For flow modelling, an area in the vicinity of the installation site was chosen which is located between Ramsey Island and the mainland with an area of approximately 800 × 1800 m 2 . The bathymetry data was collected during the "Celtic Odyssey Operation", using an echosounder [44] and then MATLAB V4 grid data method which uses the Green functions of the biharmonic operator to find the minimum curvature interpolation of irregularly spaced data points was employed to generate the seabed geometry [46]. The quality of data was not assured in all locations as some of the features of the seabed are not visible in the boat data. This is due to the tide and sea conditions preventing the boat from approaching some areas of the sea. Thus it was decided to use two additional datasets to recreate the bathymetry. One of the datasets was provided by SeaZone . This satellite data has a resolution of 30 m which means that the surface that is created by this dataset is very smooth and does not include the pinnacles and small rocks. The second dataset used was a historic dataset made available by the Royal Navy and was measured using a lead line in the early 20th century. In the regions where less information was available from the other two datasets, this was used to confirm the geometry of the seabed.
To create the final geometry the whole area was divided into 4 regions. For each region the quality of information in each dataset for that zone was investigated. For example, the quality of echosounder data was not acceptable for the deep valley and the Royal Navy dataset was not very accurate for the same region as well, so it was decided to use the SeaZone dataset to construct that region. The same evaluation was done for the other regions and the final geometry of the seabed was constructed. Afterwards the proposed geometry of the seabed was discussed with the locals and boat skippers who were familiar with the area and the bathymetry was finalised. The result gives us a bathymetry which is not only relatively simple to mesh but also resembles the key features of the region accurately, see Figure 2. The elevation for the selected area varies between +5 m and −74 m with the zero altitude datum set to be the highest peak on the seabed. The decision to have 5 m of water above the peak is based on the tidal range measured by the nearby Milford Haven tide gauge and gives the correct height of the high tide inside the sound. This distance was also discussed with local vessel operators who agreed with this assessment of the high tide clearance. In order to capture the eddies generated by the seabed and Horse Rock a fine discretised domain is considered that consists of more than 36.3 million elements (combination of tetrahedral, brick, pyramid and surface triangular). In this mesh, the vertical size of the first four layers of elements above the seabed is set to be 25 cm and then doubled for each level above up to 5 m height (i.e., fifth layer 50 cm, sixth layer 1 m, etc.). The smallest element in the flow direction is nearly 0.5 m and the largest element is slightly over 10.8 m. However, the mesh has been constructed so that over 80% of the elements are restricted to be below 2.1 m in length. This geometry is used for both the kand LES turbulence models.
The simulation uses the commercial code ANSYS FLUENT . The inlet velocity for the kcase is set to be 1.5 m/s with a turbulent intensity of 20% and turbulent viscosity ratio of 100. For the bed boundary condition a no slip wall with roughness height (k s ) of 5.4 m and roughness constant of 0.9 is selected. In ANSYS FLUENT the roughness constant is a value to define how uniformly the sand grains are distributed on the surface. A suitable value is usually between 0.5 and 1.0 and the smaller number means that the sand grains (here rocks) are placed on the surface more uniformly [47]. ANSYS FLUENT uses the following equation to calculate the roughness height: Roughness Height = (9.793 × Roughness Length)/Roughness Constant. In which the roughness constant is a value to define how uniformly the sand grains are distributed on the surface, the roughness length is calculated based on the Nikuradse Sand Roughness. Based on this equation the size of the grains on the seabed is averaged to be roughly 0.5 m.
ANSYS FLUENT uses the logarithmic law of wall to solve the wall-bounded turbulent flows [33]. The reason for selecting these values (5.4 for roughness height and 0.9 for roughness constant) is that the model has simplifications in terms of the exact shape of the bed and also in terms of the inlet boundary condition. In the real world, the inlet has a very turbulent boundary condition because of the big rocks located in the south of the sound (known as the bitches).
To model the free surface, a flat plane symmetry boundary condition (i.e., zero shear non slip wall) is applied to the upper surface of the domain. This rigid lid approach has been used in various research papers to model the turbulent flow in channels [48][49][50]. It is also possible to model channels with a rigid but sloping lid [51,52], where the slope corresponds to a drop of 1/1000, for example, as observed in the field. This has not been applied here due to the complexity of the bed forms.
To model the subgrid-scale turbulence for the LES, the Smagorinsky-Lilly method is selected. Again the same inlet velocity as for the kmodel has been chosen for the flow. By using the Vortex Method to create random fluctuations at the inlet, 1000 vortices are applied at the inlet. A turbulent intensity of 20% along with the turbulent length scale of 1m is chosen. Other boundary conditions are identical to the kcase.
A constant time step size of 6 seconds is selected. This value is calculated by multiplying the time that it takes for the flow to pass through the smallest element in the streamwise direction by a factor of 10. The maximum number of iterations for each time step is set to be 40. To run this simulation the HPC Wales platform is equipped with Intel Sandy Bridge E5-2690 at 2.9 GHz and 4 GB memory per core. The total problem is run for 5000 time steps and it takes about 466 h to run the LES job on 32 cores which is the maximum number of cores available to run the simulations. Figure 3 shows the results of velocity contours for the kmethod. It clearly shows the influence of Horse Rock and Pony Rock forming a wake in the flow field. The largest velocity reduction occurs behind the main feature of the sound, Horse Rock (775 m from the inlet), and this can be see in the data collected by an ADCP transect as shown in Figure 4. The ADCP data shown here is from a single transect and was collected from a moving vessel at a low sample rate of 1 Hz. The boat data only shows the velocity for only one transect across the channel. This data is noisy due to the condition of the sea on that particular date. Averaging several transects at a given location would provide data with improved noise characteristics. The averaged data can provide a better understanding about the state of the sea, and deliver a better comparison. Due to the transient nature of the model, vortices and eddies are constantly generated and dissipated, thus the results of consecutive iterations could show more fluctuations.

Results
The LES result ( Figure 5) shows that at least one eddy is forming behind Horse Rock. Figure 4 also shows two bodies of water coming from each side of the Rock and circling around the generated eddy before joining each other again further downstream. In Figure 6 it can be seen that the eddy is very flat. It starts its movement upward but then it seems that because of the fast water flowing overhead, it is crushed.    The simulations are conducted with a rigid lid representation of the ocean surface. Changes in sea surface level in the vicinity of some underwater obstruction, in this case Horse rock, effect flow characteristics in this region. The quantitative effect of this phenomena, and thus the accuracy of the prediction using a fixed lid, is assessed by studying the pressure change across Horse rock. The pressure field is sampled at the fixed lid at locations of maximum pressure upstream and minimum pressure downstream of the subsurface obstruction, as shown in Figure 8. Change in height resulting from the pressure difference can be predicted using the hydrostatic pressure relation between fluid height and pressure: p = ρgh, where ρ is the density of the fluid, g is gravity, and h is height. Rearranging for h using values of 9.81 m/s for gravity, 1000 kg/m 3 for density and 9763 N/m 2 for the pressure drop, h is calculated to be 0.995 m at a distance of 61 m downstream of Horse Rock. The samples are taken at x = 449, y = 5, z = 361, where p = −9014 and x = 449, y = 5, z = 300, where p = 749. Given the scale of the model and the small amount of blockage we can assume that the results are not significantly affected by the rigid lid. Therefore, the significant computational overhead of a free surface is not required in this case.

Conclusions
In this research the flow modelling and visualisation of a real bathymetry was investigated. The purpose of this model is to apply these methods to give a realistic view about the behaviour of flow in coastal regions before the deployment of tidal turbines. Two different turbulence models have been investigated for a possible site deployment of TST and the model results were compared with the real flow data acquired from the location. It has been understood that the kmethod can predict the flow pattern with relatively good accuracy near the main features of the domain and although the LES model has the ability to demonstrate some expected flow patterns because of the bathymetry, it still needs more investigation before considering it for a possible robust method for solving these kinds of problems. Due to its nature, the LES model could provide more precise results in this current problem if there had been more knowledge about the inlet boundary conditions. In ocean modelling the target area is very extensive compared to a usual application of CFD and the vortices created by the vortex method at the inlet are not propagating all along the length of the domain, therefore it is suggested that a direction for future work would be the definition of a suitable bed boundary condition that can generate such turbulent structures in the flow.
Having the free surface feature helps the model to behave more realistically near the surface. It was shown in Figure 7 that the forming eddy is propagating towards the surface and because in the simulations presented here a rigid lid was placed on the top surface, the propagation of it is restricted. With a free surface model included the eddy may continue to develop toward the sea surface and emerge as a swirl of the surface. However the analysis of the pressure difference across the subsurface obstruction demonstrates a small change in height relative to the distance downstream which is not significantly affected by the rigid lid assumption. This demonstrates a reasonable level of confidence in the simulation results. The simulations showed that although using traditional kand LES models can give reasonable results, in order to have a more precise model that can cover all the aspects of coastal flows it is necessary to customise the bed roughness parameters in a way that can compensate the effects of simplification of the bathymetry.