Modelling Deep Green tidal power plant using large eddy simulations and the actuator line method

The Deep Green technique for tidal power generation is suitable for moderate ﬂ ows which is attractive since larger areas for tidal energy generation hereby can be used. It operates typically at mid-depth and can be seen as a “ ﬂ ying ” kite with a turbine and generator attached underneath. It moves in a lying ﬁ gure-eight path almost perpendicular to the tidal ﬂ ow. Large eddy simulations and an adaption of the actuator line method (in order to describe arbitrary paths) are used to study the turbulent ﬂ ow with and without Deep Green for a speci ﬁ c site. This methodology can in later studies be used for e.g. array analysis that include Deep Green interaction. It is seen that Deep Green creates a unique wake composed of two velocity de ﬁ cit zones with increased velocity in each wake core. The ﬂ ow has a tendency to be directed downwards which results in locally increased bottom shear. The persistence of ﬂ ow disturbances of Deep Green can be scaled with its horizontal path width, D y , with a velocity de ﬁ cit of 5% at approximately 8 e 10 D y downstream of the power plant. The turbulence intensity and power de ﬁ cit are approximately two times the undisturbed value and 10%, respectively, at 10 D y .


Introduction
Tidal energy is emerging as a potential provider of sustainable electric power generation. The used technologies span from devices that operate at the surface or at mid-depth to turbines mounted on the bottom [1]. In this work the tidal power plant Deep Green developed by Minesto operating at mid-depth is studied, see Fig. 1. Tidal power plants operate in turbulent currents which makes a robust design of the equipment and control system important. The surrounding tidal mean flow and its large-scale motion is fairly well-known and can be observed using, e.g., Acoustic Doppler Current Profiler (ADCP) [2,3], but the small-scale motion (both spatially and temporally) is more demanding to observe and less known [4,5]. Here computational fluid dynamics simulations can be used to study the turbulence characteristics and its impact on power plants. (e.g., dynamical loads on the structure [6,7] and the corresponding requirements on the control system algorithms). Turbulence characteristics can beside design consideration also be important for the power quality from a tidal power plant [8]. It is a parameter that is increasingly important for power grids with an increasing share of renewable energy.
Some studies that have paid particular attention to the turbulence characteristics of tidal currents [9,10], and bottom-mounted turbines [7,11,12], have mainly been focusing on the turbulence intensity, the length scale of turbulent eddies, and to some degree anisotropic feature of turbulence. Other ways to characterize turbulence by quantities such as probability density functions, and intermittency have been suggested as well [e.g. Ref. [13]]. Previous Large Eddy Simulations (LES) results of the undisturbed tidal flow for the Deep Green test site at Holyhead in Wales show that both the turbulence characteristics and horizontally averaged boundary layer profiles differ between the accelerating and decelerating phases, for a given instantaneous mean flow [14].
Tidal and wind power plants with horizontally mounted turbines generate a wake of reduced flow and changed turbulence characteristics downstream of the plants. The flow around, and the wake behind, wind power plants and arrays of power plants have been studied extensively, [e.g. Refs.
[15e19]]. It is seen that the wake extension depends on the roughness of the ground, and that the return to undisturbed flow conditions is faster over rough land than over smooth seas. Furthermore, it has been shown that landbased wind farms increase the mixing in the lower atmosphere [20,21]. From geometrical reasons it can be argued that shallow water conditions increase the importance of horizontal mixing compared to vertical mixing, especially for arrays. This difference for tidal power plants (in shallow waters where the flow is bounded by the seabed and sea surface) compared to wind power plants (in atmospheric conditions where the flow is only bounded by the ground) stresses the importance to study tidal power plant array design and configuration, and how the local surrounded area is affected. Aspects of ocean energy extraction have been the subject in a few studies [e.g. Refs.
[22e28]], but not for the Deep Green technology.
The working principle of Deep Green is similar to a wind kite. Deep Green is pre-programmed to fly in an ∞-shaped trajectory, referred to as lemniscate, see Fig. 1. The actual path is controlled online via a control system that takes temporal and spatial flow variability into account. The hydrodynamics of the wing enables it to move almost perpendicular to the flow at a speed several times the flow speed, while the water is pressed through the turbine. The high speed-ratio between the kite and the flow enables the power plant to operate in currents with lower speed and to be smaller, using less material compared to many other tidal power technologies. Multiple Deep Greens will eventually form plant arrays. Its main components are shown in Fig. 1: the wing, the axial turbine, the nacelle that comprises the generator and power electronics, the rudders and struts, and the tether that attach it to the foundation on the seabed. The electricity from the generator in the nacelle is transferred ashore to the grid via the tether and the local power grid. The design that is studied here and tested outside of Wales has a wing span of 12 m and a generator with a rated power of 500 kW. The test site for the Deep Green is situated outside Holyhead at the east coast of Wales. The depth and maximum tidal current are approximately 80 m and 2 ms À1 , respectively.
In the present study, the Actuator Line Method (ALM) is used to model Deep Green in a tidal flow predicted by large eddy simulations. ALM has been in use for many years in e.g. wind power industry in order to estimate power output and to optimize plant farm design [e.g. Refs. [17,18]]. This is however the first time it is adapted and used to model a moving tidal kite in a LES. ALM uses a blade-element approach to model for example a complete rotating turbine with its individual blades. The main advantage with the method is that it is much less resource demanding than a set-up that resolves the geometry, since a coarser mesh can be used which in turn results in both less computational cells and enables longer time steps. This method can then be used to study tidal power plant array configurations and environmental effect. The adaptation and eventually the simulations are here performed in the OpenFOAM open source software. The wake characterization is part of present work whereas environmental and array analysis are outside the scope. The results are presented for a typical tidal power generation situation along with a resolution discussion. This is to the best of our knowledge the first time the wake characterization of a moving power plant using large eddy simulations has been presented.
The aim of this paper is to present our initial work using an adaption of ALM to predict the scale of the wake and downstream flow disturbance of a moving tidal kite using LES. The analysis is, however, not meant to be used for estimating the power output or efficiency of the power plant. The balance of resources (domain extent and mesh resolution) has therefore been focused on having a large enough domain and fine enough resolution to capture and resolve the large scale turbulence that is supposed to affect any downstream power plants rather than resolve the detailed flow in the vicinity of the power plant and its wing. Furthermore, the presented simulations are not meant to capture all aspects of tidal flow and the hydrodynamics of the power plant. The physics have been simplified neglecting any effects of a free-surface, tidal channel topography, density stratification, and tidal-cycle variations [29,30]. The power plant is simplified modelling its main components which is the wing and turbine, but not the nacelle, actuators and the tether.

Methodology
Large eddy simulations and the actuator line method are used to analyze tidally oscillating turbulent flow with the tidal power plant Deep Green. This section starts by presenting the numerical set-up including governing equations, turbulence model, boundary conditions, and computational domain. It is concluded by a description of the modified ALM used to model the Deep Green.

. Governing equations
The filtered Navier-Stokes equations are Here uðu; v; wÞ denotes the resolved velocity vector and its components in the x-, y-, and z-direction, respectively. p is the pressure, r is the density, n is the kinematic viscosity, and t is time. The sub-grid stress tensor Τ is modeled via a one-equation eddyviscosity concept, using a transport equation for the sub-grid-scale kinetic energy and a local length scale [31,32]. The body forces F T and F DG are used to represent the tidal and Deep Green forcing, respectively. The body force is activated when developing the fully developed tidal flow without the Deep Green, using cyclic conditions in the flow direction (referred to as Precursor LES). A T is the force amplitude, u ¼ 2p= T is the frequency where the period T is set to 12 h, and e x denotes the unit vector in the x-direction (general flow direction). When introducing and studying a Deep Green unit, cyclic conditions are not used in the flow direction, since it would represent an infinite array of Deep Green units once the flow disturbances reach the downstream boundary and in turn the upstream boundary. Thus, the flow with Deep Green is given by a time-resolved inlet boundary condition obtained from the Precursor LES. The tidal body force term is deactivated and the Deep Green force term F DG (described in Section 2.2) is activated. This is referred to as Deep Green LES.

Computational domain and boundary conditions
Two computational domains (rectangular cuboids) with different sizes are used in this study. A smaller domain (4 Â 2 Â 1H in the x-, y-, and z-directions), where H ¼ 80 m is the depth, is used in a sensitivity study whereas a longer domain (16 Â 2 Â 1H) is used for the detailed studies of how the Deep Green affects the flow further downstream. The sensitivity cases, see Table 1, are referred to with a "S" whereas the large domain cases are referred to as "DG" (with Deep Green) and "0" (for undisturbed flow without Deep Green). The computational domains are discretized using equidistant meshes with cubic cells. The top boundary, modeled as a flat surface, has a slip boundary condition whereas the bottom boundary has a rough wall boundary condition. Cyclic boundary conditions are applied in both horizontal directions for the Precursor LES but only in the y-direction (cross-flow) for the Deep Green LES. In the Deep Green LES the pressure boundary conditions are given by a homogeneous Neumann boundary condition at the upstream (x ¼ 0Þ boundary and Dirichlet boundary condition at the downstream boundary. For the velocity and sub-grid scale properties, a temporally and spatially varying boundary condition obtained from the precursor simulation (with tidal forcing) is applied at the upstream boundary, and a homogeneous Neumann boundary condition at the downstream boundary.
Initial conditions for the Precursor LES are given by mapping a developed velocity field from a tidal cycle large eddy simulation. That simulation was performed during several tidal cycles in order to reach fully developed turbulent and oscillating tidal flow conditions. The mapped velocity field represents the conditions at a time t 0 where the volume averaged flow speed over the vertical extension of the lemniscate is approximately 1:6 ms À1 , which is the design flow for the Deep Green power plant under study (Fig. 2). The computational domain in the tidal cycle LES used to create the initial field was (16 Â 4 Â 1H) using (2096 Â 512 Â 128) grid points in (x; y; z) directions. Dgrid is the side length of the cubic cells of the mesh, N is the number of elements in the ALM to discretize the foil, ε is the width of the Gaussian function in the ALM, Db is the ALM foil element width. a The cells are cubic and for case S4eS6 a refinement box extends (0 : 250 m; À 68 : 68 m; 12 : 68 m). b S1 and S4 use all three Gaussian width constraints whereas S2-3, S5-6 and DG only use the drag and lift constraints. c As S3 but with a differencing scheme using a blend of 98% second-order linear and 2% first-order upwind differencing for the convection terms of the momentum and turbulent kinetic energy equations.

Solver and discretization
The large eddy simulations are performed using OpenFOAM, employing a collocated finite volume approach. The PIMPLE algorithm is used with one outer corrector step and two corrector steps for each time step. The time derivatives are discretized using an implicit second-order scheme that uses the present and two previous time-step data, referred to as the backward scheme in OpenFOAM. The advection and diffusion terms are in general discretized using a second-order central differencing scheme [33]. However, for certain cases with small Gaussian width of the ALM, oscillations are observed in the velocity and pressure field as the flow is affected by the ALM source term. To remove these oscillations a differencing scheme using a blend of 98% second-order linear and 2% first-order upwind differencing is introduced in the sensitivity study (Section 2.3) and used in large domain cases 0 and DG for the convection terms of the momentum and turbulent kinetic energy equations. For ALM the time-step should for stability reasons preferably be in the order of the time it takes for a point in the ALM discretization to pass through one cell. This sets the time

Post processing and sampling of mean quantities
A large portion of the results in Section 2.3 and 3 will be presented as spatially and temporally averaged values. Here a temporal and a spatial mean of a quantity f is denoted as f and hfi x , respectively, where the subscript x describes the direction/volume of the spatial average. This means for example for temporal averages of the resolved velocity in the x-direction that: CuD is the spatial average over the complete domain (volume average), CuD x;y is the spatial average over the horizontal plane z in the xand y-directions (area average), and CuD x;y;z l is the spatial average over the depths covered by the lemniscate in the x-and y-directions (volume average).
The sampling of quantities in the Deep Green LES is performed when the wake of the Deep Green is fully developed. For the small domain cases (S1eS6 in Table 1 Table 1) the sampling starts after 400 s which equals approximately 20 full lemniscate revolutions and that the disturbances from Deep Green have been advected a distance of at least 10 widths of the lemniscate, D y . The sampling and averaging is thereby done during 600 s (between t 0 þ 400 s and t 0 þ 1000 s). The sampling time is a compromise of two conflicting constraints; a) long sampling time in order to get representative time averaged values given the temporal variability, b) short sampling time since it is an accelerating flow with changing forcing and velocity profile. The chosen sampling time length is considered to be reasonable keeping in mind that the main goal is to find the approximate rather than precise extent of the flow disturbances and wake downstream of the power plant.

Actuator line method
The actuator line method [e.g. Ref. [17]] is used to model the wing and axial turbine of the Deep Green. The deepGreenFoam solver that is implemented in this work is based on the libraries of turbinesFoam [34], which is a stand-alone user-contributed module for OpenFOAM. ALM describes the body forces that arise when Deep Green sweeps through the domain, using a blade-element approach.
The body force used in Equation (2) is obtained as the summation of each of the blade elements i, as where F foil is the actuator line force and r is the vector between the actuator point and the cell center where the source is to be applied, see Fig. 3. ε is the width of the spherical Gaussian function. The function is used in order to avoid too high spatial force gradients (which may cause numerical oscillations) between adjacent cells. The actuator line force is determined as per span-wise unit length for foil element i where "foil" is used as a general term for any component modeled with ALM. Here U rel;i represents the magnitude of the relative velocity vector, and C i is the chord length, of each foil element. The unit vectors e L and e D are in the lift and drag directions. The relative velocity is determined as where v f is the foil spatial velocity. C L;i ¼ C L ða i ; Re i Þ and C D;i ¼ C D ða i ; Re i Þ are the lift and drag coefficients given in lookup tables, a i is the angle of attack, and Re i ¼ c i U rel;i =n is the Reynolds number. The coefficients are given for 0 < a < 31 . The center foil element typically operates in 9 < a < 13 , whereas the tip elements typically operates in 5 < a < 20 . The Gaussian width ε is in turbi-nesFoam dynamically and individually determined for each time step and foil element as the maximum value of the three constraints based on the 1) mesh size as ε mesh ¼ c mesh 2DL cell where is related to the cell volume, 2) drag as ε drag ¼ c drag C d C=2, or 3) lift as a function of the chord length as ε chord ¼ c chord C where c chord ¼ 1=4. The Gaussian width and its determination is further discussed in the sensitivity analysis in Section 2.3 since it affects both the sampled fluid velocity that represents the free stream velocity and in turn the tip vorticity of the foil.
The lift and drag coefficients for the wing of Deep Green have been determined by steady-state analysis in OpenFOAM using the k À u turbulence model. The wing profile is close to a NACA6315 profile with similar lift and drag coefficients. They were determined for ten sections, see Fig. 4 and then linearly interpolated to the N elements in Deep Green LES. The coefficients were determined using a free-stream velocity of 12 ms À1 (similar to the speed of Deep Green through the water) while varying the angle of attack. C L and C D are known to be fairly insensitive to the Reynolds number, Re ¼ U rel c=n, as long as it is high enough [35]. Hence, the lift and drag coefficients are here considered as constant with respect to the Reynolds number. The lift and drag coefficients are set identical for all foil elements adjusting the chord lengths (3.3e0.5 m) accordingly. In ALM, u is typically sampled at the point along the actuator line where the force is to be applied (usually C i =4 downstream of the leading edge along the chord). This is reasonable since then everything concerning each actuator line element is taken into account in the same point and the local fluid velocity is used. It is, however, problematic in the sense that F foil;i is determined using C L;i and C D;i under the assumption that U rel;i is obtained using a free-stream velocity, while the sampled fluid velocity here in most situations is influenced by the body force in the preceding time steps. u in this position can therefore not be considered as a true free-stream velocity. How much it is affected depends on the element velocity v f , actuator line force, the Gaussian width which in turn is affected on its three constraints, and by the LES set-up via the time step and local velocity. A longer time step, and higher local flow and element velocity make the free stream velocity less affected of the body force from previous time steps and vice versa.
The temporally and spatially varying position, orientation (transformation matrix M), and velocity of the foil are given as input to deepGreenFoam. They have been determined in a prior control system and force simulation using a dynamic system model. During these simulations the model of Deep Green consisted of the kite itself (including wing, nacelle with turbine, and actuators), the tether, the control system and the tidal current environment. The flow was set to a constant plug flow of 1.6 ms À1 . The obtained input data is determined in the center of the foil 1=4 chord length downstream of its leading edge. The position is in deepGreenFoam transformed using M to any actuator line element along the actuator line. The Deep Green wing (wing span of 12 m) is discretized with N elements ( Table 1). The axial turbine with its relatively small diameter of 1:5 m is considered to be a point source and its thrust is therefore described with one element. The remaining components of the Deep Green (tether, nacelle and actuators) are not included. The influence from the nacelle and the actuators are believed to be quite small in comparison to the wing and the turbine. The tether was not included since it was still under development and its crosssectional profile affects how to set it up using ALM to a large degree. There is also an elongation of the tether that might make the modelling challenging and it was left out for time being. The center of the lemniscate are x c ¼ 20:2 m; 128:8 m for the short (4H) and long domain (16H), respectively, y c ¼ 0 and z c ¼ 47:3 m. The lemniscate horizontal width D y and vertical extension D z are 64:4 m and 22:4 m.

Sensitivity study
This work is performed during the design and development phase of the Deep Green and there is therefore up to now no measurements available of the wake of the kite. Accordingly, the performance of the modified ALM is evaluated based on a sensitivity study rather than comparing against observations. The main aim of this work is to estimate flow disturbances at some distance downstream rather than in the vicinity of the Deep Green. The results are therefore mainly both evaluated and shown at some   4. Snapshots of the iso-surface of vorticity magnitude in black: U ¼ 1 s À1 and the actuator force in green, and the velocity magnitude U at the boundaries for t ¼ 400 s for simulations (S1, S2, etc) described in Table 1. distance of Deep Green rather than in the close vicinity of Deep Green. The model results are here not used to design the Deep Green as such but how it in a large scale affects the flow. The sensitivity study is performed by varying model settings, that are known to affect the solution in order to find a set-up, where the velocity deficit in the downstream wake is fairly independent of these settings. Here, the Gaussian width, mesh resolution, position of local fluid velocity sampling, time step, and number of elements to describe the wing are altered (Table 1).
Figs. 4e7 present the results of the sensitivity study. The figures are to some extent overlapping but are though given in order to facilitate the understanding of the three-dimensionality of the flow. Fig. 4 shows snapshots of the velocity magnitude U at the side and downstream boundaries, together with iso-surfaces of vorticity magnitude U ¼ 1 s À1 . Fig. 5 shows the time averaged velocity magnitude at z ¼ z c and snapshots of the vorticity magnitude U ¼ 1 s À1 and 0:2 s À1 . Fig. 6 shows plots of the time-averaged velocity normalized with the inflow velocity, in four cross-wise transects downstream of Deep Green. Fig. 7 compares all the sensitivity cases, showing the time-averaged normalized velocity for one cross-wise transect at x ¼ x c þ 4D y and z ¼ z c .
In S2 and S5 (see Table 1) the mesh-size Gaussian-width constraint is relaxed compared to S1 and S4. The relaxation leads to a decreased Gaussian-width ε and in turn ε=Dx ratio (except for the most central elements of the foil). In Figs. 4 and 5 it can be seen that this change leads to a substantial increase in tip-vorticity, comparing S1 and S2. There is no tip-vorticity of the magnitude U ¼ 1 s À1 in S1, whereas in S2 the length of the iso-surface is approximately the same as the lemniscate length. This length is in turn of the same magnitude as the downstream extension of the isosurface of tip-vorticity in the resolved steady-state case, or even more so for S3. The same trend with increasing vorticity is seen comparing S4 and S5. In Fig. 6 it can be seen that there are large differences in the velocity deficit close to the power plant (x x c þ 2D y ), comparing S1 with S2 and S4 with S5. In both Figs. 6 and 7 it can, however, also be seen that the velocity deficit is fairly similar (somewhat smaller for S1 and S5) for all these cases at x ¼ x c þ 4D y .
A small Gaussian width ratio ε=Dx (S2, S3, S5 and S6) results in increased force gradients in the projected ALM body force and in turn in local oscillations in the velocity fields. These oscillations are visible in Fig. 4 as fringes at the vorticity iso-surfaces and in Fig. 5 as "shadows" of vorticity in front of the power plant. To counteract these oscillations a differencing scheme using a blend of secondorder linear and first-order upwind differencing is used for the convection terms of the momentum and turbulent kinetic energy equations in S3*. This scheme limits the linear scheme towards an upwind scheme using a coefficient that sets the amount of limiting. A coefficient of 1 gives the strongest limitation and 0 tends to give a second-order linear scheme. A similar limitation is used by [36] using a coefficient of 0.1 in the vicinity of the foil and 0.02 in the rest of the computational domain. Here the coefficient is set to 0.02 (small limitation) which removes the oscillations while still not significantly altering the wake results comparing S3 with S3* in Figs. 4e7.
Another approach to decrease the oscillations is to increase the mesh resolution. This is done in S4 via a mesh refinement box (cell side length Dx ¼ 0:3125 m instead of Dx ¼ 0:625 m) in the area of the Deep Green. This refinement box is extended further downstream to also allow a more resolved flow in the area downstream of Deep Green. By comparing S2 with S4, it can be seen in Figs. 4 and 5 that the oscillations disappear. Furthermore, the iso-surface of the tip-vorticity has a smaller core diameter, and it persists much longer. The velocity deficit seen in Fig. 6 is, however, fairly similar with a slightly larger deficit for the coarser mesh case S2 for x x c þ 2D y .
The influence of the position of the sampling point of u is studied by comparing S2 with S3 and S5 with S6. It is either sampled, as usually is done in ALM, at a position one fourth of the chord length downstream of the leading edge (case S2, S5), or at the leading edge (not presented in this paper), or not sampled at all but given as a fixed velocity equal to the mean stream-wise flow speed u ¼ ð1:6; 0; 0Þ in case S3 and S6. The velocity in the first two positions are affected by the source term from previous time steps, although to a different degree, whereas it is not in the third alternative. Comparing S2 with S3, and S5 with S6, it can be seen in Figs. 4e7 that the effect on the vorticity and the velocity deficit close to the Deep Green is large but that the velocity deficit, once again, is fairly similar further downstream (x ¼ x c þ 4D y ).
The time step dependence is studied by decreasing Dt ¼ 0:1 s in S3* to 0:05 s in S3**. By comparing S3* and S3** in Figs. 4e7, it is seen that the decreased time step did not significantly alter the results.
Martinez et al. [37] show that the power output using ALM for wind power begins to converge for an actuator point resolution Db=Dx 0:75, where Db is the width of the discrete blade element.
In the present sensitivity analysis, the wing is divided in either 10 or 40 elements and the minimum mesh size is 0.625 m (S1eS3) and 0.3125 m (S4eS6). This results in that the sensitivity analysis spans an actuator point resolution of 0:48 < Db=Dx < 1:92. Although not studied in a specific case comparison, it can be seen that all cases except S1 is close or within the proposed actuator point resolution of 0.75 (Db=Dx ¼ 0:48 for the cases using N ¼ 40 and a mesh size of 0.625 m which is used for the case DG presented below). Fig. 7 shows the velocity deficit at x ¼ x c þ 4D y for all cases in the sensitivity study. It is seen that the deficit is fairly similar for all set-ups (some differences is seen along the center line y ¼ 0 and for S1 and S5 in the peaks at y=D y ¼ ±1=4). It is therefore concluded that results on velocity deficit further downstream of the Deep Green do not depend critically on finer mesh resolution or shorter time step. For computational demand reasons it is acceptable to use the unrefined mesh and a time step Dt ¼ 0:1s. Furthermore, it is seen that the location of the point where u is determined in Equation (7) does not significantly affect the velocity deficit at x ¼ x c þ 4D y .

Results and discussion
Case DG and 0 are used to study a longer portion of the disturbed flow downstream of Deep Green ( Table 1). The set-up follows the conclusions from the sensitivity study in Section 2.3. Some of the details of the numerical settings are repeated for clarity (last bullet point valid for case DG only): The convection terms of the momentum and turbulent kinetic energy equations are discretized using a scheme that blends 98% second-order linear and 2% first-order upwind differencing. No local mesh refinement. A long domain of 16H is used, to allow studies of a longer wake. The sampling time for mean velocities is increased to 600 s, starting at t ¼ 400 s, since the analysis is performed to study an extended wake in a longer domain.
The Deep Green wing is described with 40 elements and the Gaussian width is determined using the drag and lift constraints. The fluid velocity is sampled at a position one fourth of the chord length downstream of its leading edge. The distance from the inlet to the center of the lemniscate is x c ¼ 2D y .
The instantaneous velocities at t 0 (initial conditions) at depth z c for case DG and 0 are shown in Fig. 8 to give a visual impression of the undisturbed flow field. Fig. 8 shows that the velocity fluctuations are larger in the mean flow direction (x-direction) than in the crosswise and vertical directions. It is seen that the largest flow structures are found in the x-direction, followed by the y-and z-direction. Furthermore, there seems to be no locking of turbulent eddies due to the domain size. The integral length scales in the horizontal directions are estimated to be approximately 0:4 < L x < 0:7H and 0:1 < L y < 0:15H over the lemniscate vertical extension for the x-and y-direction, respectively, using the autocorrelation function of u and v velocity similar to [9,38]. It is challenging to estimate the vertical integral length-scale for w using the same methodology since there is a vertical gradient in the mean flow. Due to this and since the power plant is affected by the horizontal variability in the vertical velocity as well, we therefore here estimate the horizontal length scale of the vertical velocity instead of a true vertical integral length-scale. This is found to be in the same magnitude as L y . Although difficult to directly match the sizes of the structures in Fig. 8 with the integral length-scales, it can still be seen that the overall trend with the largest scale in the x-direction is seen. Furthermore, it can be concluded that the domain extension in both horizontal directions seems to be large enough for estimating the length scales since the domain size divided by the integral length scale is approximately 20e30 compared to advisable ratios indicating 2e6 or more [38]. At the same time it can be concluded that the mesh resolution is fine enough ðL x =Dx z60; L y =Dx z15 Þ to resolve these turbulent structures at the size of the horizontal length-scales. Experiences from field tests (pers. comm.) have shown that it's the turbulent structures of the  same scale as the wingspan that affects the control system the most. These scales are here typically of the same magnitude as L y and therefore well resolved. An overview of the flow for case DG is presented in Figs. 9 and 10. In Fig. 9a) it is seen that a substantial amount of vorticity (compared to bottom-generated vorticity) is generated by Deep Green. The tip vortices last some distance downstream, where they diffuse and break up in a more irregular vorticity pattern. The forces from Deep Green acting on the flow have components in the z-directions which in turn yield a vorticity plume that is directed downwards. The flow returns to a more undisturbed flow further downstream, at approximately x > 10D y . There is, however, still additional vorticity compared to the inlet conditions in the upper parts of the domain all the way to the end of the domain. Fig. 9b) shows that the downward redirection of the flow leads to an increased shear stress at the bottom. It can also be seen in the xy-planes that the velocity disturbance persists all the way to x ¼ 10D y , however, much more pronounced at x ¼ D y and 4D y . Fig. 10 presents a close-up of the flow in the vicinity of Deep Green. All sub-figures of Figs. 9 and 10 are from the same time step and as can be seen in Fig. 9, the Deep Green is moving downwards and towards the reader (out of the paper) in Fig. 10. The center of the green ring, which is a slice of the iso-surface of the magnitude of the AFM force F DG at this particular xz-plane, can be seen as the center of the Deep Green at this time. Besides the present position, three "previous" passages of the Deep Green can be seen. In Fig. 10a) the streamlines and the increase in velocity magnitude indicate the expected velocity increase at the low-pressure side of the wing. The velocity component plots in Fig. 10bed) shows how this is manifested in each component as the Deep Green wing moves along its lemniscate and orientation. The velocity gradients in the advected disturbances are of the same magnitude as in the area of the present position of Deep Green which indicates that the mesh resolution is fine enough to maintain the flow structures as they advect downstream. It can also be seen that there are no visible oscillations in the pressure or the velocity fields as discussed in Section 2.2. It can furthermore be noted that the tests with two times higher resolution (case S4eS6) typically gives higher local vorticity but that the mean velocity deficit further downstream is similar.
In sections 3.1 to 3.4 Figs. 11e16 are used to discuss specific aspects of the results. Fig. 11 show the three-dimensionality of the flow using the time averaged velocity deficit. Figs. 12, 13 and 15 show the velocity, power deficit, and turbulence intensity, as averaged profiles at different x-positions. The results are here normalized with the results from the undisturbed case 0. The extension of the lemniscate, z c À D z =2 < z < z c þ D z =2 and À D y =2 < y < D y =2, is shaded. Panel a) in Figs. 12, 13 and 15 show the results along ÀH < y < H in order to study the cross-wise spatial variation, and panels b,c) show the vertical variation. In the latter panels, the variable has been horizontally averaged over the lemniscate width À D y = 2 < y < D y = 2. The variables are normalized with the horizontally averaged variable at the center of the lemniscate in panel b), and as a function of depth in panel c). Fig. 14 shows the time averaged velocity and turbulence intensity at z ¼ z c . Fig. 11 shows that the velocity decreases most directly downstream of the lemniscate, but only to a limited degree (or even increases depending on height) downstream of the inner parts of the rings of the lemniscate (y=D y z±1=4). As can be seen in the z-plane in Fig. 11a) and in the x-planes in Fig. 11c) the maximum deficit is found at the centerline at approximately x c þ 0:5D y and zzH=2. The position (along the centerline) with the largest deficit is found where the Deep Green passes the center of the lemniscate twice during each revolution. It can further be seen that the position of the maximum deficit for x > x c þ 0:5D y is shifted from the centerline to y=D y z À 1=4 and zzz c . This shift is also seen in Fig. 11a) as a transition where eventually the deficit is composed of two deficit zones at z ¼ z c for x ! x c þ D y . The strength of these deficit zones is then further decreased by entrainment of surrounding flow with increasing downstream distance x. This can also be seen in Fig. 12a) where there are large cross-wise gradients in the velocity deficit in the vicinity of the Deep Green, but that these gradients decrease with increasing downstream distance. The deficits are smaller than for horizontal-axis tidal turbine using similar analysis methodology but using D y for Deep Green in this study and the turbine dimeter for the latter as length scales [23]. It should, however, be noted that D y is much larger than the turbine diameter for corresponding rated turbine power. Furthermore it can be seen in Fig. 11 that while the position of the maximum velocity deficit is positioned at approximately z c the zone with increased velocity is directed downwards. A surface (used for visualization of the results) with a variable height z, that approximately follows this zone, is shown upper left in Fig. 11b). This height is indicated as a red line in the y-and x-plane plots. It can be seen that the increased velocity zone persists to approximately x c þ 2D y where it reaches the bottom flow area. At similar downstream distance and height, it can be seen in the surface plot that the velocity deficit zones end. Fig. 12a) shows that the velocity deficit varies significantly in both the x-and y-directions. The vertical variation of the velocity deficit is also given in Fig. 12b and c). Here it is seen that there is a pronounced deficit peak at z=Hz0:55 for x ¼ x c þ 2D y , whereas the deficit is more evenly distributed over the lemniscate height for x ! x c þ 4D y . The increased velocity close to the bottom can be seen in Fig. 12c). It results in increasing bottom shear stress as seen in Fig. 9b). The largest increase (approximately 15%) compared to the undisturbed case is seen at z=Hz0:025 and between x c þ 4D y < x < x c þ 6D y . It should here be noted that the prescribed inlet velocity boundary condition may lead to an overestimation of this effect compared to a pressure driven simulation with a free surface.

Power deficit
The power deficit is presented in Fig. 13, under the presumption that the available power PfU 3 . It shows (Fig. 13a) a large cross-wise variation close to the power plant and remains a substantial length downstream thereof, however, more evenly distributed. Fig. 13b) shows that the horizontally averaged available power is more evenly distributed over the lemniscate vertical extent. It is also seen that the available power is higher in the upper and lower portions of the water column. Fig. 13c) shows that there is a deficit of almost 20% in the lemniscate area for x x c þ 6D y and that it is fairly evenly distributed over the trajectory operating depth for x ! x c þ 4D y . The power deficit is decreased to approximately 10% at x ¼ x c þ 10D y . Fig. 13a) shows, however, also that the available power actually increases beside the lemniscate of the power plant. A similar feature can be seen in Fig. 13c) where the available power increases below and above the power plant. These zones with increased available power arise due to the blockage of Deep Green, causing some of the flow to take a path around the power plant. Here it should be noted that the cyclic boundary conditions at the sides of the domain can be seen as a case where there is an infinite number of power plants in a single-row array. For a single power plant, the flow velocity increase besides the lemniscate would most likely be smaller but covering a larger cross-wise distance from the actual power plant. That scenario needs however other boundary conditions and is not studied here. Although not tested here, it can be worth-while to take advantage of the power increase outside the lemniscate area using a staggered arrangement if designing a power plant array with multiple rows. Similar findings have been seen for wind power arrays [e.g. Ref. [39]]. Another way to take advantage of the increased velocity outside the lemniscate can be to alter the operating depth in consecutive rows. This is, however, for geometry reasons not believed to be the first option to choose for tidal power due to navigation constraints etc.

Turbulence intensity
The turbulence intensity is often used to characterize turbulence while evaluating tidal power sites. It is typically defined as the square root of the variance of the velocity divided with its time Fig. 11. Time averaged velocity deficit u DG =Cu 0 D x;y . y ¼ ÀDy=4 and z ¼ zc are indicted with black lines. a) z-plane at z ¼ zc and y-plane at y ¼ À Dy=4. b) Upper panel shows a surface at variable height z, indicated as a red line in lower panel showing an y-plane at y ¼ ÀDy=4 and in the panels in c). c) x-planes at xc, xc þ Dy=2; xc þ Dy; xc þ 3Dy= 2; xc þ 2Dy, and xc þ 5Dy=2 at right. u DG =u 0 x;y ¼ 1 is indicated with a dashed line.  ðCu DG D yl =Cu 0 D x;y Þ 3 , where Cu 0 D x;y;zc and Cu 0 D x;y are the horizontally and time averaged velocities during undisturbed conditions, for the zc-and each z-planes, respectively. The velocity is sampled at the xcyc-plane and five xy-planes downstream thereof (x ¼ xc; x ¼ xc þ 2Dy; x ¼ xc þ 4Dy; x ¼ xc þ 6Dy; x ¼ xc þ 8Dy; andx ¼ xc þ 10DyÞ. a) Sampled at transects at zc. b-c) Spatially averaged over À0:5 < y=Dy < 0:5 (shaded area in a)). averaged value, TI ¼ s=U, where the variance is defined as s 2 ¼ u0 2 þ v0 2 þ w0 2 , s is the standard deviation, and u 0 ¼ uÀ u etc., are the fluctuating parts of the velocity components. The turbulence intensity estimation is fairly straightforward for flows that vary from a steady mean value. The flow cases in this study, however, are spatially inhomogeneous with a velocity gradient from bottom to surface, and temporally varying. For the cases with Deep Green there is also a moving power plant that affects the flow making it further inhomogeneous. Here we will discuss how the definition of mean velocity influences the turbulence intensity estimation. In Fig. 14a) the time averaged velocity U DG is presented. The low U DG values found in the center of the trajectory (y ¼ 0) and approximately 0:5D y downstream thereof result in high local turbulence intensity when used for normalization as seen in Fig. 14b). In Fig. 14c), instead of U DG the mean velocity CU 0 D x;y from case 0 without Deep Green is used for normalization of the standard deviation. It can be seen that the turbulence intensity is still high downstream of Deep Green in Fig. 14c) but much smaller than in Fig. 14b). The use of the undisturbed flow in case 0, as in Fig. 13c) can be seen as using undisturbed flow conditions upstream of a power plant during a measurement situation and will be used hereafter. In Fig. 14d) the turbulence intensity s DG =CU 0 D x;y is normalized by the turbulence intensity Cs 0 D x;y =CU 0 D x;y for the undisturbed flow. This can be simplified as s DG =Cs 0 D x;y which equals the standard deviation at each cell in a specific layer for case DG normalized with the horizontally averaged standard deviation for all cells in the same layer for case 0. It can in Fig. 14d) be seen that the turbulence intensity is substantially increased downstream of the power plant (s DG =Cs 0 D x;y ! 3 in the center of the domain for x x c þ 6D y ).
The turbulence intensity s DG =CU 0 D x;y;zc is also presented as transects in Fig. 15ab). It can be seen in Fig. 15a) that the turbulence intensity has a large cross-wise variation close to the power plant. The turbulence intensity and its cross-wise variation decrease downstream and the turbulence intensity is fairly homogenous in the wake for x ! x c þ 4D y . Fig. 15b) shows the spatial average over the lemniscate width in order to encompass the large crosswise variation of the turbulence intensity. It shows that the turbulence intensity in the lemniscate area is larger in the upper parts of the vertical extension of the lemniscate at x ¼ x c and in the lower part at x ¼ x c þ 2D y . The turbulence intensity is then eventually more evenly distributed further downstream. In Fig. 15c), the turbulence intensity s DG =CU DG D y l is shown. It can be seen that the turbulence intensity at x ¼ x c þ 2D y peaks at z=Hz0:38 which coincides well with the position of the transition from increased to decreased velocity in the wake seen in Fig. 11b). At x ¼ x c þ 10D y the turbulence intensity is approximately a factor 2.5 larger for case DG than for the undisturbed case 0 in the upper part, and 1.5 in the lower part of the lemniscate vertical extension.

Reference length scale of disturbances
The persistence of a flow disturbance is often estimated as a ratio of reference length scales, where this ratio typically is in the order of 10. Typical reference length scales for Deep Green could be the turbine diameter, D Turbine , the span of the wing, D wingspan , or the vertical or horizontal extension of the lemniscate D z or D y , respectively. It is seen Fig. 15a) that the lemniscate area-averaged velocity and power deficits are approximately 5% and 10% at 8D y and 10D y , respectively. 10D y corresponds to approximately 30D z , 50D wingspan , and 400D turbine . The lemniscate width therefore seems to be a representative reference length scale for Deep Green. Another option is to use an equivalent diameter of the lemniscate. Using the swept area converting it to an equivalent circular area gives D ¼ ð4D y D z =pÞ 0:5 , which for the present set-up equals 0:7D y , which in turn also seems to be a reasonable reference length scale. This equivalent diameter is here based on that the swept area is calculated as a rectangle with the sides D y and D z . It is, however, not defined how to characterize the swept area for Deep Green, and for time being, D y can be used as the reference length scale. Note that other parameters, such as flow velocity and depth, have not been altered in this study. These parameters will most likely influence Fig. 15. Turbulence intensity depth profile (z/H) presented as the normalized velocity standard deviations s DG =CU 0 D x;y;zc , s DG =CU 0 D x;y , and Cs DG =U DG D yl , where CU 0 D x;y;zc and CU 0 D x;y are the horizontally and time averaged velocities during undisturbed conditions, for the zcand each z-planes, respectively. For the undisturbed case the full domain length is used for averaging. The velocity is sampled at the xcyc-plane and five xy-planes downstream thereof (x ¼ xc; x ¼ xc þ 2Dy; x ¼ xc þ 4Dy; x ¼ xc þ 6Dy; x ¼ xc þ 8Dy;and x ¼ xc þ 10DyÞ. a) Sampled at transects at zc. b-c) Spatially averaged over À0:5 < y=Dy < 0:5 (shaded area in a)). the results in general but the choice of reference length scale only to limited degree. Besides comparing velocity and available power, a turbulence intensity comparison between case DG and 0 can be seen in Fig. 16b). It is seen that the increase in turbulence intensity persists longer than for the velocity and power deficit. It is still approximately 1.5 times or more larger for x x c þ 10D y .

Summary and conclusions
The Deep Green technology, here studied, is suitable to use in moderate tidal flows which makes it promising since larger areas of tidal flow generation then can be foreseen. It is a novel technique that presently is undergoing full scale testing. In the development phase, the ability to study turbulent characteristics of the flow is important in order to find a good ratio between cost and robust design and plan for eventual power plant arrays. Turbulent characteristics are here studied by using large eddy simulations for a specific site outside Wales, with Deep Green modeled using the actuator line method. The actuator line method is here further developed in order to be able to model the arbitrary trajectories of the power plant. The main emphasis during this study has been to discuss the flow features at some distance downstream of Deep Green rather than in the close vicinity thereof.
Factors affecting the performance of the actuator line method is discussed in a sensitivity study. The relative flow velocity can be acquired in different ways. It is, however, seen that using a sampling point of the local flow velocity at either the leading edge of the foil or at the 1/4 chord length downstream thereof does not influence the extension of the wake significantly. In a mesh sensitivity study it is found that the mesh size (used here) influences the flow characteristics in the vicinity of the Deep Green, but not significantly for the wake analysis further downstream. The Gaussian length, ε, in the Gaussian function used to distribute the body force in the actuator line method influences the tip vorticity to a large degree. Decreasing its length by decoupling the Gaussian length and the mesh size may, however, come with a price of numerical oscillations in the velocity field in the vicinity of the actuator line. This is due to more locally applied body forces along the actuator line resulting in larger force gradients. These oscillations are suppressed by using a blend of a second-order central differencing scheme (98%) and a first-order upwind scheme (2%) for convection terms of the momentum and turbulent kinetic energy equations.
The flow disturbance length scale is determined using the deficit in the mean flow velocity. It is found to be 8 times the horizontal lemniscate width, D y , for a velocity deficit of approximately 5%. It can furthermore be seen that the Deep Green has a tendency to direct the flow downwards which results in increased bottom shear. It is also shown that the turbulence intensity is substantially increased downstream of Deep Green. It is almost two times the undisturbed values up to 10D y downstream of the power plant. The power deficit is approximately 10% at the same distance. The power deficit and turbulence intensity increase can be used to discuss how to balance costs for installation and maintenance against energy output for different power plant mutual distances. The velocity deficits are significantly smaller for Deep Green than horizontalaxis tidal turbines at corresponding number of length scales downstream of the power plant. It must, however, be noted that the length scale D y used here is much larger than the turbine diameter, often used as a length scale for horizontal-axis tidal plants at equal rated power. This difference in length scale must be taken into account while doing comparisons between different tidal power plant designs and array distributions.
Future development of the power plant model may include the inclusion of the tether and anchoring depending on its design and size. As eventually more components of Deep Green, not considered here, may be included an enhanced mesh resolution analysis would be fruitful. That is especially true if design optimization studies are to be performed. The overall set-up may be further developed by adding surface effects for shallow sites and topographical feature or even large boulders if large enough to influence the flow field. Present model can be used for array design and arrangement analysis and for environment studies regarding turbulence effects and flow changes along the sea bottom. It would in that respect also be interesting to perform further research regarding how Deep Green is affected by different levels of turbulence.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.