Simulating tidal turbines with multi-scale mesh optimisation techniques

Embedding tidal turbines within simulations of realistic large-scale tidal ﬂows is a highly multi-scale problem that poses signiﬁcant computational challenges. Here this problem is tackled using actuator disc momentum (ADM) theory and Reynolds-averaged Navier-Stokes (RANS) with, for the ﬁrst time, dynamically adaptive mesh optimisation techniques. Both k − ω and k − ω SST RANS models have been developed within the Fluidity framework, an adaptive mesh CFD solver, and the model is validated against two sets of experimental ﬂume test results. A brief comparison against a similar OpenFOAM model is presented to portray the beneﬁts of the ﬁnite element discretisation scheme employed in the Fluidity ADM model. This model has been developed with the aim that it will be seamlessly combined with larger numerical models simulating tidal ﬂows in realistic domains. This is where the mesh optimisation capability is a major advantage as it enables the mesh to be reﬁned dynamically in time and only in the locations required, thus making optimal use of limited computational resources.


Introduction
This study focuses on the extraction of tidal stream energy from coastal waters via horizontal axis tidal turbines which are currently the favoured approach to efficiently harness the vast and reliably predictable tidal resource. The deployment of tidal turbines is a complex and expensive operation and this makes the task of locating the optimal position for such turbines even more important. Maximising the power output of arrays of turbines is essential, but the environmental impacts must also be studied and modelled in depth as it is vital to ensure that the efforts to reduce carbon emissions do not result in new environmental concerns. Previous studies have shown that in order to correctly assess the power extraction from tidal turbine arrays, an undisturbed flow approach, also termed the flux method [1], does not suffice and the hydrodynamic influences of the turbines and their wake interactions must be accounted for [2,3,4,5]. Therefore, a numerical model that aims to examine the power output and environmental impacts of tidal turbine arrays must be able to capture these features.
Currently many large-scale marine hydrodynamic models employed to study marine energy use the depth-averaged shallow water equations. In order to numerically simulate arrays of turbines, these models usually adopt a discrete approach where the turbines are represented as a region of increased bottom drag [6] and the low complexity of such an approach allows for adjoint-based optimisation techniques to be used to improve turbine positions [7]. More recently, Funke et al. [8] extended this to find the optimal turbine density, i.e. the number of turbines per unit area represented as a continuous field. These models benefit from relatively low computational costs that allow for large scale realistic simulations, such as modelling arrays of tidal turbines in the Pentland Firth [9]. However, the main shortcoming of such an approach is that it can fail to account for important turbulent physics and 3D effects, e.g. since the flow passing below and above the turbine is not modelled. Alternatively a higher fidelity 3D model coupled with an appropriate turbulence model can be used. A number of these models have been developed and validated against experimental flume tests with promising results [10,11], but their high computational expense has generally prevented their application in larger scale regional simulations.
In a recent study, Masters et al. [12] investigated the performance of a range of computational models of horizontal axis tidal turbines at different spatial scales. As part of that study, both a high resolution blade element momentum (BEM) model and a large scale coastal model were used to simulate the flow past a small tidal turbine array at an idealised headland. It was demonstrated that while the flow velocities in the far upstream were very similar, substantial differences exist in the wake profiles of the two models. This is not surprising given the differences in spatial and temporal scales used, as well as the different treatment of turbulence in the two approaches. Ultimately, it was concluded that the choice of model will depend on the physical scale of interest and the computational resources available [12].
Mesh optimisation techniques can help bridge the gap and improve the accuracy of large scale simulations without the need for excessive computational power. Creech et al. [13] utilised dynamic mesh optimisation to develop a high fidelity ADM-LES model to accurately model the flow past wind turbines. Furthermore, Divett et al. [6] developed a 2D depth-averaged model coupled with an adaptive mesh flow solver to demonstrate the greater energy extraction that can be achieved from turbines arranged in staggered layouts.
The model presented here is based on the actuator disc momentum (ADM) theory, where the turbines are represented as momentum sink terms. This 3D approach is coupled with a mesh optimisation algorithm that employs a fine spatial resolution only in regions of interest. This allows for an accurate turbine wake characterisation whilst maintaining a relatively low computational cost. An ultimate aim will be to use this model to assess the effects of deploying tidal turbine arrays in realistic domains, such as the Inner Sound of the Pentland Firth where MeyGen Ltd plan to deliver a fully operational 398MW renewable energy plant powered purely by the tide [14].
An alternative approach to ADM theory would be to apply BEM theory, whereby radially varying forces in both axial and azimuthal directions are applied. This method better represents the performance of a real turbine as it introduces a swirl component in the wake profile. Stallard et al. [15] investigated the mean wake properties behind a single three-bladed rotor and demonstrated that between 0.5-2 diameters downstream of the turbine, BEM theory can account for the near wake properties reasonably well. However, Batten et al. [16] demonstrated that the swirl component of the wake dissipates quickly in the streamwise direction and the far wake profiles produced by ADM and BEM methods are very similar. Given that the aim of the current study is to develop a tool for array design, as long as the turbines are not positioned in very close proximity to each other, neglecting swirl remains a reasonable assumption.
The paper is organised as follows. First in section 2 the Fluidity ADM model is presented along with a description of the turbulence models used and the mesh optimisation techniques employed. Then in section 3 verification against the Conway solution for flow past an actuator disc [17] as well as a comparison against a similar OpenFOAM model is presented. This is then followed by two more realistic test cases, sections 4-5, where the Fluidity ADM model is validated against experimental flume test results. The paper concludes with a general overview of the results and a discussion on the benefits of the model presented.

Methodology
The numerical model presented has been developed within the Fluidity framework, an open source finite element CFD code with 3D mesh optimisation capabilities [18]. The main feature of Fluidity that motivates this study is its ability to adapt the mesh dynamically in time and only in the locations of interest. The code is also highly parallelisable which makes it attractive for larger scale computationally challenging applications.

Governing equations
One of the main challenges when attempting to simulate the flow past a turbine is the ability to correctly account for the turbulence within the flow. This is vital, given that the ambient turbulence intensity has a significant effect on the structure of the turbine wake and its recovery as demonstrated by Blackmore et al. [19] and Mycek et al. [20,21]. For the purpose of this study it has been decided to incorporate turbulence models based upon the Reynolds-averaged Navier-Stokes (RANS) approach, given that their computational cost is significantly lower than that required with large eddy simulations (LES). These models are based on the RANS equations, in which the velocity is decomposed into mean and fluctuating (turbulent) components where ⊗ denotes the outer product, u is the mean velocity, u is the fluctuating velocity, p is the mean pressure, µ is dynamic viscosity and S u is the momentum sink term included here to account for the presence of the turbine. The third term on the right hand side, containing the Reynolds stress tensor −ρ u ⊗ u , represents the effect of turbulent fluctuations on the mean flow and for incompressible flows is defined as where k = (u · u )/2 is the turbulent kinetic energy and µ T is the dynamic eddy viscosity.
The 3D Fluidity model presented uses a P1 DG − P2 finite element pair [22] to discretise the RANS equations. This scheme uses the space of discontinuous piecewise linear functions (P1 DG ) to represent velocity and the space of continuous piecewise quadratic functions (P2) for pressure. This is a stable element pair on tetrahedra which is compatible with the geophysical balances, hydrostatic and geostrophic, [22] that become important at the larger scale where the turbine simulations will be embedded in the future.
The Fluidity model uses a Crank-Nicolson time discretisation which is second-order accurate in time. However, in order to achieve a stable and bounded solution with the discontinuous Galerkin (DG) discretisation, a slope limiter is used in conjunction with an explicit treatment of the advection term. Therefore, the RANS equations are considered in two stages where initially advection is considered, and then the other terms. For advection, explicit subcycles utilising adaptive sub-timesteps are used in order to satisfy a specified Courant-Friedrichs-Lewy (CFL) condition. A larger overall timestep and Crank-Nicolson discretisation are used for the remaining terms including diffusion, pressure gradient and sources. For further details, the reader is referred to [23,24].

Actuator Disc Momentum theory
The 3D numerical model developed and validated here incorporates turbines which are parametrised based on the ADM theory outlined by Houlsby et al. [25]. ADM theory is based on the assumptions that the flow is inviscid and incompressible with uniform inflow. The disc is infinitely thin and the thrust loading on the disc is uniformly spread. Given that many of these assumptions do not hold in the numerical model, or the real world, it is vital to verify the implementation and to validate the outcome from the model against laboratory and real world data.
In the current model the circular disc has a small finite thickness and is represented as a scalar turbine field which is unity at the location of the disc and zero everywhere else in the domain. This scalar field is discretised using the same space as the velocity field; i.e. P1 DG . This ensures that the disc is present only in the intended region and does not spread into the surrounding elements, which would be the case if a continuous Galerkin (CG) discretisation was used. In order to set the appropriate loading on the disc, the Fluidity model uses the established definition for thrust coefficient, C T , to compute the magnitude of thrust loading that should be applied at the disc. This is uniformly spread across the volume of the disc and is implemented as a momentum sink term where S u is the momentum sink applied only at the location of the disc, ρ is the fluid density, A disc is the cross-sectional area of the disc, V disc is the volume of the disc and u 0 is the unperturbed upstream streamwise component of velocity. In the context of the numerical simulations, it is important to ensure that u 0 is predicted accurately; this being discussed in detail as part of the model validation in sections 4-5.

The k − ω model
One of the most widely used RANS models is the k − ω model where the momentum equations are closed by solving transport equations for k and the turbulent frequency, ω, where u is the mean velocity (with the overbar dropped for convenience from here onwards), µ T = ρk/ω and with σ k = 0.5, σ ω = 0.5, β = 0.075, β * = 0.09 and α = 5/9 [26]. S k and S ω represent additional source terms which are not present in the original k − ω model, but are included here to account for the presence of the turbine.
Alternatively, the k − ε RANS model is a popular approach which can be used where the momentum equations are closed by solving transport equations for k and the turbulent dissipation, ε, where ε = β * ωk. However, both of the k − ω and k − ε models tend to predict faster wake recoveries in comparison with experimental data, as demonstrated by Roc et al. [10]. The reason for this lies in the fact that they fail to account for the energy transfer rate from large-scale turbulence to small-scale turbulence in the near-wake region. This is known as the short circuiting of the turbulence cascade due to the presence of the turbine. To overcome this, turbulence correction terms have been suggested in both the tidal and wind turbine literature [10,27,28]. Rethore et al. [28] and Roc et al. [10], include additional terms for the turbulent kinetic energy (TKE) equation (4). They suggest including an additional source term S k,p which accounts for the production of wake turbulence and also include an additional sink term S k,d to account for the short circuiting of the turbulence cascade [28]. In this study the approach by Rethore et al. has been followed where the source and sink terms are scaled with C T and applied only at the disc. These are calculated via with and a = 1 where a is the axial induction factor [29], β p = 0.05 and β d = 1.5 [28]. El Kasmi et al. [27] suggest an additional production term for the k − ε model which is proportional to the quadratic production of TKE by shear and is applied at the disc ±0.25 diameters directly upstream and downstream of the disc. Rados et al. [30] extended this for the original k − ω model and compute an additional source term, S ω , via with C ω = 4 [30]. This is the definition of S ω used in this study.

The k − ω SST model
It has been suggested that a k − ω model is better suited for modelling separating flows compared to the k − ε model due to the latter's inability to capture turbulence correctly in near-wall regions [26]. However, despite the k − ω model's superior performance in the near wall region, its strong sensitivity to the freestream ω values has prevented it from becoming the standard model of choice for RANS modelling [31]. This shortcoming was one of the main motivations for the development of the k−ω shear stress transport (SST) model by Menter [32] whereby blending functions are introduced so that the k −ω model is used in the inner region of the boundary layer and the k − ε model is used in the outer region and in free shear flows. Furthermore, the definition of the eddy viscosity is modified to account for the effect of the transport of the principal turbulent shear stress [32]. In the k − ω SST model the transport equations are modified to where the blending function F 1 is used to switch from the k − ω model, when F 1 = 1, to a k − ε model, when F 1 = 0 and the limiting is applied in order to prevent the build-up of turbulence in stagnation regions [33]. Furthermore, all the coefficients in the transport equations (10) and (11), are computed using a blend from the coefficients of the original k − ω and k − ε RANS models. For example with σ k1 = 0.85, σ k2 = 1, σ ω 1 = 0.5, σ ω 2 = 0.856, α 1 = 5/9, α 2 = 0.44, β 1 = 0.075, β 2 = 0.0828 and β * = 0.09 [33]. For completeness detailed descriptions of the blending functions and the eddy viscosity calculation have been provided in the appendix.

Mesh optimisation
The Fluidity ADM model is capable of dynamic mesh optimisation which is used to help reduce discretisation errors by refining the mesh in locations of numerical complexity or specific interest, e.g. regions with high velocity shear. This is achieved via the construction of a metric tensor field based upon the Hessians of solution fields, user-defined error bounds, and maximum and minimum allowed edge lengths (l e ). Thereafter, using this metric, a functional of the mesh elements is formed whose minimisation through heuristic topological operations ensures a mesh with elements which are appropriately shaped and sized. Consequently, all the field data from the previous mesh is projected onto the new mesh by interpolation. As a result of this the mesh is not only refined in regions of interest, it is also coarsened in regions where high resolution is not required. This helps reduce the computational cost over a tidal simulation where the regions of interest within the domain continuously change.
In the interest of brevity, a detailed description of the mesh optimisation techniques is omitted here, and reference is made to the extensive work by Piggott et al. [18] and Pain et al. [34]. For the purpose of this work the mesh is refined in regions of high curvature in the velocity, k and ω fields, motivated by the desire to correctly capture the re-energisation of the wake downstream of the turbines. While mesh optimisation techniques by [18] and [34] are now well established, they have never been demonstrated with the turbulence models described under sections 2.3-2.4. Indeed, the combined application of accurate turbulence modelling and efficient mesh optimisation forms one of the key advances of the present work, and is central to the description of the flow fields in the wake of tidal turbines.

Semi-analytical solution
The implementation of the ADM model in Fluidity is novel and must be verified. Conway [17] suggests a semi-analytical form for the velocity profile of the wake behind an actuator disc and this will be used to help verify the numerical ADM implementations; the approach adopted being similar to that by Viré et al. [35]. Assuming incompressible and inviscid flow, the velocity profile takes the form if r < D/2 and where u x is the velocity in the streamwise direction, Q −1/2 (w) denotes the Legendre function of the second kind, Λ 0 (b, h) denotes Heuman's Lambda function, x is the streamwise direction, r is the axial direction and D is the disc diameter. The only unknown in the solution is u wake , the streamwise component of velocity in the wake which is sufficiently far downstream from the turbine that the pressure can again be treated as uniform. In this study u wake has been computed using [25]: The other parameters are computed from x, r and D as follows:

OpenFOAM
To provide an additional means of comparison, the same ADM implementation described above has also been applied to the PISO (Pressure Implicit with Splitting of Operators) solver in OpenFOAM, an open source finite volume CFD package [36]. This was achieved following a similar approach to that described in [37]. OpenFOAM is a well-established CFD package and the OpenFOAM ADM model will help verify the implementation in Fluidity. It will also aid in identifying the potential benefits of the finite element approach used in Fluidity. As with the Fluidity model, the disc has a small finite thickness and is represented as a scalar field which is unity at the location of the disc and zero everywhere else in the domain, with the momentum sink term, (3), applied at the disc only. Furthermore, the finite volume discretisation used in OpenFOAM ensures that the turbine scalar field is only present where intended, similar to the P1 DG finite element discretisation used in the Fluidity model described in section 2.2.

Comparison
In order to computationally reproduce the inviscid flow assumption of the Conway solution, a low value for kinematic viscosity was used in both Fluidity and OpenFOAM, ν = 1 × 10 −6 m 2 s −1 . The two numerical models have used different, yet similar size meshes where a tetrahedral mesh was used for Fluidity and a hexahedral mesh was used for OpenFOAM. The domain extent and the disc size were kept the same with a 0.25% difference in disc volume between the two meshes. The total number of elements was also very similar with the Fluidity mesh containing 8.47 × 10 5 elements and the OpenFOAM mesh containing 8.50 × 10 5 elements. Cross-sections of the two meshes are shown in Fig. 1. A uniform inlet velocity, u in , is applied with the side walls, top and bottom surfaces all set to free slip and a zero pressure boundary condition applied at the outlet. Furthermore, u 0 is simply set equal to u in in order to ensure that the magnitude of S u , Eq. (3), set in both numerical models is the same.    However, adaptive sub-timesteps were chosen to ensure a maximum Courant number of ∼ 0.1 for the advection term, as described in section 2.1, which   Fig. 2. The velocity deficits are well predicted by both the Fluidity and the OpenFOAM models and the plots demonstrate that using a higher C T leads to a greater velocity deficit, as expected. Along r = 0 the Fluidity model provides a better match to the semi-analytical solution and a closer look at Fig. 3a reveals that the Fluidity model performs better at capturing the sharp velocity variation present at the disc location, whereas the OpenFOAM model exhibits some fluctuations. Along r = D the simulated results diverge from the semi-analytical solution and this is due to the infinite domain assumed in the Conway solution, whereas the wall effect in the numerical simulations causes the velocity outside of the wake, in the bypass region, to accelerate. This effect can be reduced by moving the walls further away from the disc; i.e. using a larger computational domain and a smaller blockage ratio. This is demonstrated in Fig. 3b where a larger numerical domain has been used. Increasing the domain site indeed leads to an improved agreement due to reduced blockage. The remaining departures are of order < 0.01u 0 , which is considered acceptable. In the radial direction, although both models struggle to provide a close match to the Conway solution just behind the disc, Fig. 4a, these discrepancies are drastically reduced two diameters downstream, Fig. 4b, and both models follow the analytical solution with the Fluidity results providing a closer match.
Overall, considering the assumptions used in the semi-analytical solution, both numerical models simulate the inviscid flow past the actuator disc accurately. Moreover, the P1 DG − P2 finite element discretisation employed in the Fluidity model has helped in providing a closer match to the Conway solution compared to the OpenFOAM model.

Stallard et al. test case
Having considered the model performance based upon a comparison to a simple and well-established reference solution, a much more realistic test case is considered next. For this purpose, a comparison is made to the experimental data from a study of a group of three-bladed rotors carried out by Stallard et al. [38]. These experiments were conducted in the University of Manchester wide flume, which is 5 m wide with a 12 m long test section.
The tests utilised fixed pitch rotors with a diameter D = 0.27 m which were located at mid-depth and positioned 6 m away from the inlet. The Fluidity ADM model has been used to simulate the flow past the three scaled turbines positioned in-line, with the numerical domain used shown in Fig. 5. Previously, Mungar [39] used a turbine model consisting of a momentum sink term developed within the Delft3D framework [40], Shives et al. [41] used a BEM model implemented in Ansys CFX [42] and Olczak et al. [43] implemented a RANS-BEM model in the commercial CFD code StarCCM+ to simulate these experiments and their results are briefly discussed later.

Boundary conditions
At the inlet, a Dirichlet boundary condition with constant inlet values, u in , k in and ω in is applied and the side walls and top surface are set to free slip. A zero flux boundary condition has been applied at the walls for both k and ω and a zero pressure boundary condition has been applied at the outlet. In order to simulate the flow inside a flume, it is important to capture the vertical asymmetry caused by the slower moving fluid near the flume bed. Hence, in the velocity field a quadratic drag boundary condition is applied to the bottom surface. In the Fluidity finite element discretisation, integration by parts is applied to obtain the weak form of the RANS equations, Eq. (1). This leads to a boundary integral for the stress term. For the quadratic drag boundary condition, the tangential component of this integral is set equal to C D ||u||u, whereas for the free slip boundary condition, this is simply set equal to zero [24]. This is similar to the bottom boundary condition used in [10]. The bottom drag coefficient, C D , value commonly used in marine simulations is C D = 0.0025 [4,6,7,9] and this is the value chosen for this study.
Furthermore, since upstream quantities are not always readily available, u 0 has been computed here using the velocity at the disc, u disc , and the axial induction factor, a, Eq. (8), where u disc is calculated by computing the average streamwise component of velocity of the elements making up the disc. Due to the assumptions made in the derivation of Eq. (8), this is applicable for low blockage cases only. The scenarios considered here have blockage ratios of 2.5% (single turbine) and 7.6% (three turbines), and therefore satisfy this assumption. Having said that, the local blockage ratios are much higher, especially for the centre disc, and this may affect the C T [5] which will influence the wake interactions downstream of the turbines [38].

Ambient conditions
The inlet values must be chosen carefully in order to ensure the ambient velocity and turbulence properties of the experimental flume are accurately represented via the numerical model. During the experiment, an inflow velocity of 0.47 m s −1 was used and the turbulence intensity (I) at the inlet was reported to be around 11% [38]. I is computed via and the values at the inlet were set to u in = 0.47 m s −1 and k in = 4 × 10 −3 m 2 s −2 .
In the absence of an appropriate bottom boundary condition, the inlet turbulence intensity, I, decays downstream and is dependent on the inlet ω value. Fig 6 displays the variation of I in the streamwise direction for different ω in values. It is evident that the quadratic drag boundary condition leads to minimal bottom turbulence generation and doubling C D does not do enough to help maintain I in the streamwise direction. A better representation can be achieved via the use of wall functions for the k and ω fields similar to those available in commercial CFD packages, e.g. Ansys CFX, and used in [41,44]. These are currently unavailable in the Fluidity model, but are being actively pursued.
Hence, for this study, in order to ensure that the rate of decay of k is similar to that observed in the experiments, ω in = 1.0 s −1 has been chosen. In choosing the ω in value, one should also bear in mind the implications that this will have on the turbulence length scale, l. This can be approximated for both the k − ω and the k − ω SST RANS models via Therefore, ω in = 0.5 s −1 has been avoided as it suggests l = 1.41 m which is more than three times the depth of the flume.   The experimental data is compared against numerical results from the Fluidity ADM model using both the k − ω and the k − ω SST RANS models ran on fixed meshes. The numerical data matches the velocity profiles well, but slightly underestimates the turbulence intensity. Both RANS models produce very similar results for the empty flume simulations. This is not surprising because in the absence of a no-slip wall boundary condition, which is the case here, and given that there is no pressure discontinuity in the flow, due to the absence of any turbines, the F 1 blending function in the k −ω SST model equals unity almost everywhere in the domain. Consequently, the k − ω SST model behaves very similar to the original k − ω model.

Single turbine
Having established a match with the ambient data, initially the flow past a single turbine is considered. In the Fluidity model C T = 0.92 which corresponds to the upper bound of the C T value given in Stallard et al. [38]. Stallard et al. [15] investigated the mean wake properties behind a single three-bladed rotor identical to those considered in [38] and demonstrated that for distances greater than 8D downstream the velocity deficit becomes two-dimensional and self-similar. The self-similar form follows a Gaussian profile ∆u ∆u max = exp −ln(2) where ∆u = u 0 − u x and ∆u max u 0 = 0.864(x/D) −1/2 − 0.126 , where R is the turbine radius [15]. Recently, Stansby et al. [45] used this to demonstrate that the superposition of self-similar velocity profiles can lead to accurate predictions of depth-averaged wake velocities downstream of arrays.
The results from the fixed mesh Fluidity ADM k − ω model, with the aforementioned correction terms, Eq. (6) and Eq. (9), and the ADM k − ω SST model have been compared against this self-similar profile, Fig. 9. The Fluidity ADM wake profiles also resemble a self-similar solution with a good agreement observed between the lateral profiles at 8D and 12D. There are however discrepancies in the bypass region where the Fluidity results suggests that the flow is accelerating whereas Eq. (18) sets the bypass flow equal to the ambient. Furthermore, both ADM models predict a narrower profile compared to the experimental data and the self-similar profile. This could be due to the different turbulence length scales in the onset flow in the flume with vertical, lateral and streamwise scales typically in the ratio 1:3:5 [45]. This variation in length scales would not be captured via the RANS models implemented and therefore the model has underpredicted the amount of mixing in the lateral direction.

Three turbines
The flow past three turbines has also been modelled in Fluidity with mesh optimisation and using the inlet values stated above. In this section the results for the ADM k − ω model, with and without the correction terms, Eq. (6) and Eq. (9), as well as the results for the ADM k − ω SST model are presented. Fig. 10 displays the variation of velocity deficit, turbulence intensity and non-dimensional length scale, l/D, in the streamwise direction along the centreline for the Stallard experimental data and the various Fluidity simulations. Note that l is not available for the experimental data since ω was not experimentally measured.
It is evident that the original k−ω model underpredicts the velocity deficit significantly. The RANS correction terms, Eq. (6) and Eq. (9), help improve the velocity deficit prediction by increasing the peak deficit, Fig. 10a, and shifting it further downstream. The reason for this lies in the fact that the S ω term results in an increase in ω around the actuator disc which results in greater TKE dissipation, leading to a drop in the peak I value, Fig. 10b, which helps delay wake recovery. Examining the variation of turbulence length scale, Fig. 10c, reveals how the correction terms in the Fluidity ADM k − ω model help capture the short circuiting of the turbulence cascade. It also demonstrates how the Fluidity ADM k − ω SST model is also able to capture the short circuiting via the blending functions and without the need for additional correction terms. This reflects the popularity of the k − ω SST model in the literature [11,41,43,44,46], but to the best knowledge of the authors, it is the first time that a comparison between the k − ω and k − ω SST model has been presented to demonstrate the advantages of using the latter.
The velocity deficit and turbulence intensity comparisons in the lateral direction at various distances downstream of the turbines are also presented in Fig. 11-12. The root mean square error (RMSE) between the simulated results and the experimental data have also been computed and are presented in Fig. 13. Although the numerical results predict lower velocity deficit and turbulence intensity values than the experimental data, this match improves further downstream where both the k − ω ADM model with correction terms and the k − ω SST ADM model show good agreement with the experiments. This is reflected in the reduction of RMSE in the far wake (i.e. x > 8D). In fact, even the near wake results are better than those presented in [39], which is not surprising since the latter is a depth-averaged model which does not take into account 3D effects. Furthermore, the wake widths have increased compared to those presented in Fig. 9 and therefore follow the experimental results closer. It appears that the presence of the additional turbines has encouraged more lateral mixing which has resulted in wider wake profiles.

Mesh optimisation
In order to demonstrate the potential advantages of the mesh optimisation available in the Fluidity ADM model, the adaptive runs are compared against fixed mesh solutions. The fixed mesh adopted for the present test case is illustrated in Fig. 14. For the purpose of the adaptive simulations, this mesh is refined in regions of high curvature in the velocity, k, and ω fields, and the minimum and maximum values of the element edge length were set to l e /D = 0.037 and l e /D = 1.11 respectively. The minimum edge length specified is smaller than the finest resolution of l e /D = 0.055 used in the fixed mesh and the maximum edge length used in the adaptive simulation is slightly larger than the rotor diameter and three times greater than the coarsest resolution of l e /D = 0.370 used in the fixed mesh. Both fixed and adaptive mesh simulations ran until steady state was achieved and the resulting wake profiles were compared to ensure no loss of accuracy as a result of using mesh optimisation, Fig. 15. The RMSE between the wake profiles of Fig. 15 and the experimental results are presented in Table 1, illustrating a close agreement between the wake profiles of the fixed and adaptive mesh simulations.
The fixed mesh simulation used 3.8 × 10 5 vertices compared to a maximum of only 1.0 × 10 5 vertices used in the adaptive simulation. Despite the computational time taken up by the adaptivity algorithms, the adaptive mesh simulation still takes ∼ 32 hours compared to ∼ 71 hours needed for the fixed mesh simulation when both simulations were run in parallel on 16 cores. This adaptive runtime corresponds to ∼ 512 CPU hours which is an order of magnitude greater than the ∼ 30 CPU hours required by the fixed mesh BEM model of Malki et al. [47] to simulate the flow past a single turbine. However, the Malki et al. model is a steady state finite volume model with a carefully constructed mesh consisting of 1.5 × 10 6 elements which uses a fine resolution close to the turbine. Hence, extending the model to simulate the flow past tidal turbine arrays requires the construction of an adequate mesh prior to the simulations. On the other hand, the dynamic mesh optimisation capabilities of the Fluidity model can automatically adjust the mesh to ensure the physics is accurately captured, regardless of the number of turbines, which saves time on the pre-processing as well as the overall runtime. Furthermore, the Fluidity runtime is only a fraction of the 0.14 million CPU hours required by the moving mesh RANS model of Afgan et al. [11], although the latter is an unsteady RANS model and the substantial runtime is predominantly due to the constraint on timestep where a value of 2 × 10 −5 s had to be used [11].
An important point to note is that in this simple case, since the location of the wake can be predicted, a non-uniform mesh has been carefully created for the fixed mesh Fluidity simulation which only uses a fine resolution in locations of interest. However, the location and extent of the wake is not always known a priori in the real world, e.g. due to time dependent variations in the tidal flow, and a high resolution uniform mesh would be required potentially everywhere. In this extreme case, if the finest resolution was used throughout the domain, this would have led to 2.8 × 10 6 vertices. This is more than seven times larger and would have required significantly more computational power to achieve the same level of accuracy as the adaptive mesh simulation. This would be further exacerbated if a larger domain and a large array of turbines was considered.
In summary, the Fluidity ADM model has combined the advantages of a highly adaptive mesh CFD solver with the benefits of appropriate RANS modelling to accurately characterise the wakes formed behind the rotors. The model has successfully identified the locations within flow that play a critical role in the development of the wake profile and has optimised the mesh in these regions, maintaining a coarse resolution everywhere else, Fig. 16. This approach has ensured that the momentum and turbulence perturbations are accurately predicted, whilst keeping the computational overhead relatively low.  Table 1: RMSE values in velocity deficit and turbulence intensity comparing the fixed and adaptive mesh simulation results to the experimental data. The corresponding wake profiles are shown in Fig. 15.

Mycek et al. test case
In order to further validate the ADM model presented, a comparison between the numerical model and results from an experimental study of a pair of three-bladed 1/30th scale prototypes of horizontal axis tidal turbines has also been conducted. In this case, the turbines are axially aligned one behind another as opposed to the in-line arrangement of the previous test case. The experimental dataset chosen is taken from a study by Mycek et al. [20,21] where the effect of turbulence intensity on the wake structure of the turbines was examined. Mycek et al. conducted experiments at two different background turbulence intensities, I ∞ = 3% and I ∞ = 15%, with similar inflow conditions. The two turbines used in both cases were running at a fixed tip speed ratio, TSR = 3.67, relative to the onset flow of the front turbine [21]. Using the thrust curves provided, the corresponding thrust coefficients were determined as C T = 0.81 (I ∞ = 3%) and C T = 0.75 (I ∞ = 15%) [20].
Initially the wake behind the upstream turbine, in the absence of the downstream turbine, was recorded and then the wake behind the second turbine, positioned 6D downstream of the first turbine, was also recorded. This approach was emulated in the numerical study and Fig. 17 shows the numerical domain used in Fluidity which is based on the dimensions of the IFREMER flume where the experiments were conducted. The same boundary conditions previously described in section 4.1 have also been applied here, and Table 2 describes the inlet conditions used in the Fluidity ADM simulations. The rate of decay of turbulence in the streamwise direction was not available in this case and the ω in value chosen has been motivated by the values used in [46]. The blockage ratio in this case is 6.2% and as with the previous test case, Eq. (15) has been used to compute u 0 .  The results from the adaptive mesh simulation are shown in Fig. 18-19. The Fluidity ADM model has successfully captured the increased rate of wake recovery expected as a result of the increase in background turbulence intensity. In fact for the I ∞ = 15% case the wake of the upstream turbine is almost fully recovered 6D downstream where the second turbine is positioned and both turbines produce similar wake profiles. This is definitely not the case with I ∞ = 3% where the rate of wake recovery is much slower and the second turbine is placed in the wake of the upstream turbine which has also resulted in a greater velocity deficit.
where u x is the instantaneous turbulent velocity fluctuations in the x direction and similarly u y is the instantaneous turbulent velocity fluctuations in the y direction. Shives et al. [46] also conducted a numerical study of these experiments and describes how the experimental I 2D values can be used to derive I values normalised by the upstream velocity, Eq. (16), which take into account all three velocity components. Given that the TKE in the RANS models is based upon all three components of velocity, the experimental I values presented here are normalised following the same procedure as Shives et al.
As with the Stallard test case, mesh optimisation has been used throughout these simulations to refine the mesh in regions of high curvature in the velocity, k, and ω fields and the minimum and maximum values of the element edge length were set to l e /D = 0.063 and l e /D = 0.631 respectively. These l e /D values are not refined enough for accurate study of turbine blade turbulence. However, given that the actual turbine blades are not being modelled here, the l e /D values chosen are adequate for the ADM model adopted. Eddies are not modelled in the RANS models employed, if an LES approach had been followed, then a far finer resolution would certainly have been necessary. Also, note that the maximum edge length is only a factor of ten greater than the minimum which was found to be sufficient for accurate wake characterisation in this case. In order to ensure a smooth transition from small elements to larger elements, mesh gradation algorithms are used in Fluidity that constrain the rate of growth in desired edge lengths along an edge. Throughout this study, a conservative gradation parameter value of 1.5 has been used. The limiting factor regarding this value and the ratio of the minimum to maximum edge length depend on the numerical schemes used and the physics involved. For larger scale simulations these can be extended to higher values without jeopardising the accuracy of the simulations, although this has not been explored in this study.
As observed previously in section 4.4, the original k − ω model overpredicts the rate of wake recovery for every case considered, Fig. 20-23. In contrast, with the help of the correction terms this is delayed and the model is able to follow the experimental results much better. Once again, a good agreement with the experimental data is produced via the k − ω SST model for both low and high turbulence intensity cases without the need to include additional turbulence correction terms. Furthermore, there is very little difference between the output from the k − ω SST model and the k − ω model with correction terms, especially in the far wake. These results are consistent with [46] where an ADM model coupled with k−ω SST RANS was used. The experimental data for the I ∞ = 15% case appears to be skewed towards one side of the tank and Mycek et al. put this down to the removal of the flow conditioning screens from the flume, which was necessary to achieve the more turbulent flow. Given that this was not modelled in the numerical simulations, since sufficient data regarding the inflow was not available, the outputs from the ADM models do not agree as well as the I ∞ = 3% case. That being said, a good agreement can be seen from 6D onwards in both velocity deficit and turbulence intensity values for the upstream turbine, Fig. 21, as well as the downstream turbine, Fig. 23. In comparison to the previous test case, a closer agreement between the simulated and experimental results is observed here and the RMSE values between the simulations and the experimental data are smaller than those observed in section 4.4 with the maximum values presented in Table 3.
The test case presented here is a scenario where inaccuracies in the wake characterisation of the upstream turbine will be magnified further down-  stream. This is because the momentum sink term applied depends on the velocity at the downstream turbine location which lies within the wake of the upstream turbine. The Fluidity model has accurately modelled the wake behind the upstream turbine which has in turn led to accurate wake characterisation of the downstream turbine. Once again, this has been achieved by correctly establishing the influential regions within the domain and refining the mesh therein, Fig. 20-21. This is an important result considering turbines are likely to be placed in the wake of upstream turbines in tidal turbine arrays in the real world.
Furthermore, the RANS models deal with the different ambient conditions effectively to produce a good match with the experimental results in both low and high background turbulence cases. It is interesting to see that, in both test cases considered in this study, the results from the k − ω SST model closely follow those produced via the k − ω model with correction terms. Given that the k − ω SST model does not require correction terms, it will be the RANS model of choice. Having said that, based on the results presented, there is no disadvantage to using the k − ω RANS model, as long as appropriate correction terms are introduced to account for the short circuiting of the turbulence cascade.

Conclusions
An ADM model incorporating a momentum sink term and RANS models which can simulate flow past horizontal axis tidal turbines and account for turbulence characteristics has been developed within an adaptive mesh CFD solver. The model has been verified against the Conway semi-analytical solution and a brief comparison against an OpenFOAM model based on the same ADM methodology is also presented to portray the benefits of the P1 DG − P2 finite element discretisation scheme employed in the Fluidity ADM model. Mesh optimisation techniques have been successfully employed to obtain an accurate representation of the turbine wakes whilst maintaining a relatively low computational cost and such techniques can be used in the future to aid larger scale simulations. Furthermore, the validity and versatility of the model has been assessed by reproducing two different experimental flume tests with the numerical model producing closely matching results by capturing the important turbulence characteristics, especially in the far wake. This has provided valuable insight into the importance of accounting for the short circuiting of the turbulence cascade due to the presence of the turbine when using RANS models. It has also been demonstrated for the first time that while the original k − ω model is able to capture the short circuiting with the aid of additional correction terms, the k − ω SST model is able to capture the far wake turbulence characteristics without the need to include correction terms that will require tuning. Having said that, correction terms may be beneficial to represent processes in the near wake of a tidal turbine that are not described by an ADM approach [43].
The versatility of the model needs to be tested further via a comparison against established higher fidelity numerical models. For example, a Fluidity model capable of mesh optimisation based on ADM that uses LES for turbulence modelling has been developed by Creech et al. [13] to assess the performance of wind turbines, and Afgan et al. [11] use both RANS and LES with a 3D moving mesh model of a tidal turbine to investigate wake profiles. In order to improve the ambient turbulence computation, wall functions need to be introduced to improve the bottom boundary representation of the current Fluidity model. Furthermore, C T is assumed to be constant in this work; however thrust curves, where C T is expressed as a function of u 0 , can be used instead. This represents the performance of a real device more accurately as it enables the model to take into account the cut-in speed below which the turbine does not operate and the rated speed above which C T decreases to maintain a constant power yield [9]. This would be a worthy addition to the model when simulating realistic scenarios. The next step in this research will be to seamlessly extend this model in order to investigate flow past arrays of tidal turbines in realistic tidal flow channels, e.g. the Inner Sound of the Pentland Firth. The large scale simulations will be carried out in Fluidity using a single layered model for efficiency, where 3D dynamics are unimportant, and in regions where the turbines are located the number of layers can be increased and the ADM model will be inserted. In order to capture the transient characteristics of such realistic flows, a fixed mesh approach would be challenging since the same mesh cannot be optimally generated for both flood and ebb flows without compensating for numerical accuracy. This is a situation where the mesh optimisation techniques presented can be employed to maintain a reasonable accuracy with realistic computational expenses. Furthermore, a goal is that this 3D adaptive ADM model can be used to help validate and extend the layout optimisation approach for tidal turbine arrays, developed by Funke et al. [7], which currently uses a shallow water model.
where ν is the kinematic viscosity, y is the distance to the nearest wall and Also, the definition for the dynamic eddy viscosity is where a 1 = 0.31 and S is the strain invariant S = 2 (S · S), with the second blending function, F 2 , defined as