Judith Juvin

Victor Malherbe

Marie-Charlotte Belhomme

Stéphanie Castex

Agathe Martinez

Hassan Khartabil

Arnaud Haudrechy
email: arnaud.haudrechy@univ-reims.fr

A stereo-electronic interpretation of a challenging orthoamide Overman rearrangement rationalized by molecular modelling

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the course of our studies on the synthesis of human neuraminidase 1 specific inhibitors, 1 one key transformation has been used a modified Overman rearrangement proposed by Chida's group ten years ago. 2 This reaction was performed onto orthoamides instead of imidates and benefits from time consuming protection steps. 3 However, as observed by the authors on several examples, [3][4][5][6][7] this process has proven rather difficult. In uncertainty, we have decided in this article to detail the mechanistic aspects and to carefully study the chances of success of the different steps involved in this Overman rearrangement.

Our study has decisively been extended to include molecular modelling investigations. In fact, orthoamides are recognized as reagents included in the portfolio of Overman rearrangements, yet questions remain concerning their structures as well as their behaviors in Overman rearrangements leading to the observed reactivities.

Therefore, we have focused on understanding their structures and reactivities. The relative energies and conformer percentages have been computed to clarify the type of orthoamides that are present in the reaction media and to investigate the reaction mechanisms. From this, conclusions have been drawn concerning the mechanisms of the reactions, the role of temperature, and the origin of reactivity. This theoretical investigation is intended to rationalize, optimize, and extend our experimental observations.

A The Overman rearrangement, a sometimes demanding transformation

The Overman rearrangement 2 was described four decades ago as a [3,3]-sigmatropic rearrangement of allylic trichloroacetimidates leading to allylic trichloroacetamides (Fig. 1), and this transformation has been applied to the synthesis of elaborated nitrogenated structures, like in cases of modified -aminoacids 1, 8 a substituted pyrrolidine 2, 9 polyoxamic acid 3, [10][11][12] a tetrodotoxin intermediate 4, 13 and thymine polyoxin C 5.

Fig. Replacing trichloroacetonitrile with trifluoroacetonitrile proved to be very efficient and yields were above 90%. [9][10][11]13 The major disadvantage of this method is the gaseous character of toxic trifluoroacetonitrile, nevertheless it is possible to use it as a solution in tetrahydrofuran. 10,31,32 In addition, Metz's group 33 described the Pd(II) catalysis offering another opportunity for improvement as excellent yields were obtained at room temperature for similar starting materials and Gonda's group has observed interesting results with microwave irradiation.

34,35

However, in our hands and on large scale, the Overman rearrangement remains a challenging transformation because of unreproducible results and handling difficulties with unusual reagents.

B The orthoamide Overman rearrangement

The major difficulty, when the Overman rearrangement is chosen as a key step in a total synthesis, relies on the necessity for the isolation of a specific alcohol, with orthogonal selective protections of other functions. The orthoamide Overman rearrangement, whose proposed mechanism appears in Figure 2, initially discovered by Vyas group 36 has been nicely applied by Chida's group as a key step in the total synthesis of broussonetine F, 3,4 sphingofungin F, 5,6

and kaitocephalin (Fig. 3), 7 and in the course of the synthesis of highly nitrogenated carbohydrates, 37 under high temperature conditions.

Experimental results

Interested by the key step of the total synthesis of broussonetine F and especially by the transformation 6a / 7a, 3,4 we have synthesized the starting orthoamide 12 inspired by the method described by Chida's group, and tested the orthoamide Overman rearrangement under thermic conditions and also under micro-wave irradiation.

Despite numerous experiments, we did not succeed to repeat the reaction to obtain the benzyl substituted derivative 13 (R = Bn) as efficiently as Chida's group, and consequently, we turned our attention to the TBDPS substituted series, otherwise more adapted to the continuation of the total synthesis we have anticipated. The observed isolated yields were slightly higher under microwave conditions (36% compared to 0% with usual heating on a 100 mgscale), nevertheless disappointing in the scope of a scale-up anticipation (Fig. 4).

Please do not adjust margins

Please do not adjust margins Supported by the analysis of the samples using ERETIC (Electronic REference To access In vivo Concentrations) NMR method,

38
showing that by-products may exist significantly (estimated to approximately 40%), we have concluded that this approach is not easy at it seemed at first glance. Disappointed by the low yield observed, even after an extensive optimisation, we have decided to carefully think about the rearrangement, order to detect reasons for this rather problematic transformation.

In this article, we would to share with the readers several stereoelectronic introspections into orthoamide diastereoisomers with consequent differential behaviours, in terms of stabilities and reactivities.

A orthoamide exists as a mixture of two diastereoisomers in equilibrium

The reaction between the precursor diol 14 and trichloroacetonitrile under basic conditions gave two diastereoisomeric orthoamides 12a (R) and 12b (S), in a 1:1 ratio estimated by 1 H NMR thanks to nuclear Overhauser effects, separable by HPLC. However reequilibration after standing in solution in CDCl 3 at room temperature for 24h has been observed, as already observed by Chida's group on related derivatives (Fig. 5). Each diastereoisomer exists as a mixture of conformationally flexible structures, presumably due to the 1,3-dioxolane system, [START_REF] Coiffier | For an analogous discussion on flexibility of furanoses[END_REF] most likely controlled by the semi-equatorial positioning of the bulkier substituent (CCl 3 ), as exemplified in figure 6, and consequently, conformers and 12b-maj should be the major ones. A first supposed element in the discussion is the possible stabilizing hydrogen bond between NH 2 and OPMB in the case of 12a-maj, even if this "7-ring" is not very favorable. Consequently, we have anticipated that 12a-maj may be the most stable diastereoisomer (Fig. 7).

B Generation of the trichloroacetimidate 15 from the two orthoamide diastereoisomers

Importantly, during the reaction conditions of Overman rearrangement, each of these diastereoisomers 12a and 12b may lead to the crucial intermediate 15, bearing a trichloroacetimidate on allylic position (Fig. 8). However, this compound may be contaminated by the formation of the unproductive regioisomer 16, all these species being hypothetically and fortunately in equilibrium.

The difference in behaviours between regioisomer 15 and the "Dead end" regioisomer 16 may be explained by a careful observation of the general figure 1. Indeed, the oxygen moiety should be on allylic position, in order to perform the Overman rearrangement. If for any reason, the equilibrium is difficultly displaced back to the diastereoisomeric orthoamide mixture, and unfortunately favouring the formation of regioisomer 16, this would be a dramatic sterile dead end for the process. This fact has been previously invocated by Danishefsky, as an insurmountable issue, resulting in failure during one approach in pancratistatin synthesis. This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins We believe that all the involved species (12a, 12b, 15 and 16) deserve complete stereoelectronic and/or modelling studies to understand their behaviours towards a productive Overman rearrangement.

C Stereoelectronic considerations of the different conformers of 12a

In a first step, we only considered the conformational aspects of the 1,3-dioxolane rings and figured out that the necessary subsequent selective opening obeys to stereoelectronic considerations. [START_REF] Deslongchamps | Stereoelectronic effects in organic chemistry[END_REF] In a second step, we investigated the different orbital overlaps and for that purpose we drew the most probable staggered conformers, regarding the starting from the possibly most stable 12a (R), named 12a-i, 12a-ii and 12a-iii (Fig. 9).

Fig. 9 Most probable staggered conformers of 12a.

We postulated that the antiperiplanar orbital overlap generates the desired regioisomer 15 (indicated in green in figure 10).

Fig. 10 Only one favorable overlap for most probable staggered conformers of 12a.

We deduced from the figure 10 that 12a-iii, satisfying the antiperiplanar overlap, may be the only productive conformer.

D Stereoelectronic considerations of the different conformers of 12b

The same reasoning has been adapted to the most probable staggered conformers, always regarding the C-N bond, for the possibly less stable 12b (S), named 12b-i, 12b-ii and 12b-iii (Fig. 11).

Fig. 11 Only one favorable overlap for most probable staggered conformers of 12b.

We deduced from the figure 11 that 12b-i, satisfying the antiperiplanar overlap, may be the only productive conformer.

These first considerations have been completed by molecular modelling studies.

E Molecular modelling studies

The Overman rearrangement is a process involving several transformations (Fig. 12), with all species in equilibrium, each entity existing as many conformers. The hypothesis of our work relies on the supposition that some stable species have the disadvantage not to lead to the desired Overman transformation, which would explain the difficulty of this process. A detailed discussion of the different species, their stabilities and behaviours, underpinned with modelling studies, is described in the Please do not adjust margins Please do not adjust margins supplementary materials section. The aim of this paragraph is to emphasize the most important points supporting our discussion.

-Discussion on orthoamides The structures and reactivities of orthoamides have been investigated in order to rationalize the difficulty of this Overman rearrangement. To get a quick overview on the conformers of orthoamide diastereoisomers 12a and 12b, full geometry optimisations of the six most probable staggered conformers (Fig. 13 and Table 1) have been performed. Even if several molecular conformations have been considered in each case, only results for the most stable ones will be reported here. A first observation on this table 1 is that the diastereoisomer 12a (negative values) seems to be a little more stable than the diastereoisomer 12b. The computed difference of Gibbs free energy is small and sometimes negligible. The smallest values correspond to 12a at 298K (∆G ranges from 0.6 to 1 in absolute value) and at 453K (∆G ranges from 0.7 to 1.2 in absolute value). Nevertheless, this first theoretical consideration shows that 12a-i and 12a-ii are the most stable conformers of 12a. 12b-ii is the most stable conformer of 12b. Unfortunately only 12a-iii and 12b-i deliver productive imidates (Fig. 10 and11) whereas these conformers appear to be the most unstable according to Table 1 (even if interconversions between all conformers are clearly very rapid and solve this issue). In order to justify the stabilities of 12a-i, 12a-ii and 12b-ii, the relatively recent IGM methodology was used. [START_REF] Lefebvre | [END_REF][42][43] This tool is particularly attractive to capture and characterize noncovalent interactions and has shown a weak hydrogen bond in 12a-i (2.890 Å where usually 2.4-2.5 Å are observed, see ESI).

As an intermediate conclusion regarding the behaviours of orthoamides 12a and 12b, a priori, the different conformers have differences in energies, to a greater or lesser extent, and could react differently. Importantly, 12a-i, 12a-ii and 12b-ii are the most stable ones, while 12a-iii and 12b-i are the most favorable ones for the formation of a useful imidate in the Overman rearrangement (see Fig. 10 and11), these parameters may slow down the process.

-Discussion on prereactive complexes The possibility of weakening of C-O bonds in the prereactive complex B (Fig. 12) corresponding to the diastereoisomer 12a-i (presumably giving a "Dead end" imidate, Fig. 8 and 10) was investigated using Intrinsic Reaction Coordinate (IRC) calculations at the wB97X-D/ Def2TZVP level (Fig. 14). This study revealed that one of the C-O bond is weaker (1.406 Å compared to 1.391 Å) with a lower IGM-IBSI index (1.144 compared to 1.196). The prereactive complex B should consequently evolve to cleave this weaker C-O bond, giving the "Dead end" imidate, with the imidate function positioned on an homoallylic position (unproductive for the Overman process). The same approach allowed to study the weakening of C-O bonds in the prereactive complex A (Fig. 12) corresponding to the diastereoisomer 12a-iii (corresponding this time to the "Useful" imidate), revealing that one of the C-O bonds is weaker (1.404 Å compared to 1.388 Å, see ESI) with a lower IGM-IBSI index (1.156 compared to 1.208). The prereactive complex A should consequently evolve to cleave this weaker C-O bond, giving the "Useful" imidate, with the imidate function positioned on an allylic This journal is © The Royal Society of Chemistry 20xx

Please do not adjust margins Please do not adjust margins place compared to the double bond (productive for the Overman process).

Interestingly, the molecular modelling study indicates that the prereactive complex A resulting from the less stable conformer 12aiii is far more stable compared to the prereactive complex B coming from the more stable conformer 12a-i (Fig. 15). This fact may "confuse" and slow down the overall Overman process. ) and the addition on homoallylic position barrier lowers (54.5 -5.7 = 48.8 kcal.mol -1

) but both are similar and high (difference < 1 kcal.mol -1 = 0.8 kcal.mol -1

). This similarity in barrier indicates that for the addition reaction, there is a competition between addition on allylic position and addition on homoallylic position barrier. Although the two reactive reaction channels (allylic and homoallylic) are in competition, prereactive complexes show a clear bias for the prereactive complex A corresponding to the "Useful" imidate addition on allylic position (0.5 kcal.mol -1 vs 5.7 for the prereactive complex B). Nonetheless, the free energy of the orthoamide conformers (12a-iii with the ideal antiperiplanar orbital overlap) is not the lowest one, and the fractional population of this orthoamide conformers is small (5.37 at 453K, 4.22 at 373K, and 2.92 at 298K), which explain the experimental difficulties.

-Discussion on transition states Transition state structures (TS) have been located on the potential energy surfaces at the wB97X-D/ Def2TZVP level using the Gaussian 16 program. 44 IRC analyses have been undertaken to examine the reaction path from the transition state structures of 12a on a potential energy surface and to check the species connected by a transition state. Satisfyingly, two reactants with a possible stabilizing hydrogen bond between NH 2 and OPMB have been identified from the endpoints of IRC, which is in accordance with our preliminary hypothesis (Fig. 16). The transition state TS1, that leads to the "Useful" imidate on allylic position, has shorter distances for the C-O broken bonds (2.36 Å vs 2.57 Å) and the H-O being formed bonds (1.68 Å vs 1.99 Å) compared to the distances in the transition state TS2, that leads to the "Dead End" imidate on homoallylic position. In the second case (TS2), the C-N distance is shorter (1.29 Å vs 1.33 Å), this one looking like a semi-planar transition state. Some additional computations were also performed to accurately predict the barrier of the two reactions, leading to the imidate product substitutions on allylic and on homoallylic positions. The computed Gibbs free energies for these reactions G° ‡ are 49.2 and 54.1 kcal.mol -1 for 12a-i and 12aiii respectively. These (high) Gibbs free energy barriers may explain the requirement of elevated temperature in this Overman process.

-Overall conclusions explaining the problematic formation of imidates in our Overman process Overall, several parameters may "confuse" and slow down our Overman process. Indeed, the different conformers of orthoamides have differences in energies, to a greater or lesser extent, and can react differently. Importantly, 12a-i, 12a-ii and 12b-ii are the most stable ones. In terms of formation of a "Useful" imidate, 12a-iii and 12b-i are the most favorable ones. Additionally, the prereactive complex B resulting from the less stable conformer 12a-iii is far more stable compared to the prereactive complex A coming from the most stable conformer 12ai.

Finally, the computed Gibbs free energies for the reactions G° ‡ are equal to 49.2 and 54.1 kcal.mol -1 for 12a-iii and 12a-i respectively. These (high) Gibbs free energy barriers may explain elevated experimental temperature conditions.

Conclusion

The orthoamide Overman rearrangement is a useful transformation, affording multiply elaborated structures, however suffering from harsh reaction conditions (especially in terms of elevated temperatures). Molecular modelling studies on orthoamides stabilities and transition states evolutions buttress these observations and we believe that Please do not adjust margins Please do not adjust margins stereoelectronic effects may explain the reasons of these difficulties. The requirement of a specific steric arrangement allows to thread the way into the desired product. We hope this study will contribute to help the readers to understand their results, anticipating transformations for better processes.

Fig. 2 Fig. 3

 23 Fig. 2 Expected general mechanism for the orthoamide Overman rearrangement.

Fig. 4

 4 Fig. 4 Orthoamide Overman rearrangement under microwave irradiation.

4 5

 4 Two stable diastereoisomers due to the orthoamide formation.

Fig. 6

 6 Fig.6 Major conformers of orthoamide diastereoisomers (12a-maj and 12bmaj).

Fig. 7

 7 Fig. 7 Stabilising hydrogen bond in 12a-maj.

Fig. 8

 8 Fig.8 Equilibrium between conformers of orthoamide diastereoisomers (12a-maj and 12b-maj) and trichloroacetimidate regioisomers.

Fig. 12

 12 Fig. 12 Species and transformations in the Overman process.

Fig. 13 Tab. 1 .

 131 Fig. 13 Optimised structures (wB97X-D/ Def2TZVP level) of the conformers of 12-a and 12-b. Tab. 1. Computed energies at wB97X-D/ Def2TZVP (kcal.mol -1 ) level and dH(NH2)•••O(PMB) in Å

Fig. 14

 14 Fig. 14 Prereactive Complex B obtained from the endpoints of IRC calculations at the wB97X-D/ Def2TZVP level for 12a-i, and schematized structures with bond lengths and IGM-IBSI indexes.

Fig. 15 1 )

 151 Fig. 15 General reaction mechanism and free energy profile. Barrier heights for the two possible reactions are calculated to identify the most reactive reaction channels (allylic vs homoallylic positions). The reaction step has a high activation barrier (~50 kcal.mol -1 ) in the two cases. These (high) Gibbs free energy barriers may explain the requirement of (elevated) temperature in this Overman process and (slow) reaction time. The addition on allylic position barrier raises (50.1 -0.5 = 49.6 kcal.mol -1

Fig. 16

 16 Fig. 16 Transition states structures TS1 (left) and TS2 (right) for 12a orthoamide diastereoisomer at the wB97X-D/ Def2TZVP level.

1

  The Overman rearrangement, as reported in 1974 and examples of chiral structures elaborated with this transformation Ajouter l'abréviation de TBDMS Cyclic structures, including unusual carbohydrates are the most common targets of traditional Overman rearrangement.

15-30

Overall, highly functionalized oxygenated structures require diverse drastic temperature conditions in aprotic high boiling points solvents, such as 165-180°C in dichlorobenzene or 220°C in t-butylbenzene (in sealed tubes) or varying between 135°C under very low pressure (0.1 mm Hg) and 150°C in xylene. However, despite these harsh conditions, the product was sometimes obtained in low yields. Consequently, two strategies have been developed to overcome these drawbacks.

This journal is © The Royal Society of Chemistry 20xx J. Name., 2013, 00, 1-3 | 3

Acknowledgements

We would like to thank the Région Grand-Est for the financial support of the PhD thesis of Judith Juvin. Financial support by CNRS, "Conseil Régional Champagne Ardenne", "Conseil Général de la Marne", Ministry of Higher Education and Research (MESR) and EU-programme FEDER to the PlAneT CPER project is gratefully acknowledged. We thank the MaSCA (Maison de la Simulation de Champagne-Ardenne, France) for various supports including computing facilities (http://romeo.univ-reims.fr) and the CRIANN computational center (http://www.criann.fr) for additional support. We thank Alexis Vallée for HPLC purifications. Also, we would like to thank Prof. E. Hénon for regular and lasting scientific discussions and for his work on related systems.

Author Contributions

Judith Juvin:

Conflicts of interest

There are no conflicts to declare.


Mikaël Grondeau
email: grondeau.mikael@hotmail.fr

Sylvain Guillou
email: sylvain.guillou@unicaen.fr

Philippe Mercier
email: philippe.mercier@unicaen.fr

Emmanuel Poizot
email: emmanuel.poizot@lecnam.net

Wake of a Ducted Vertical Axis Tidal Turbine in Turbulent Flows, LBM Actuator-Line Approach

Keywords: tidal turbine, wake, turbulence, numerical simulation, Lattice Boltzmann Method, Large Eddy Simulation, Actuator Line Model

Vertical axis tidal turbines are devices that extract the kinetic energy from tidal currents. Tidal currents can be highly turbulent. Since ambient turbulence affects the turbine hydrodynamic, it is critical to understand its influence in order to optimize tidal farms. Actuator Line Model (ALM) combined with Large Eddy Simulation (LES) is a promising way to comprehend this phenomenon. In this article, an ALM was implemented into a Lattice Boltzmann Method (LBM) LES solver. This implementation gives good results for predicting the wake of a vertical axis tidal turbine placed into a turbulent boundary layer. The validated numerical configuration was then used to compute the wake of a real size ducted vertical axis tidal turbine. Several upstream turbulence rates were simulated. It was found that the shape of the wake is strongly influenced by the ambient turbulence. The cost-to-precision ratio of ALM-LBM-LES compared to fully resolved LBM-LES makes it a promising way of modeling tidal farms.

Introduction

In the context of the global energy crisis, developing renewable energies is of primal importance. Among the renewable resources available are tidal currents. Tidal currents can be exploited with tidal turbines. Tidal turbines extract the kinetic energy contained in tidal currents, mostly in areas where maximum currents are over 1.5 m.s -1 [START_REF] Multon | Energies Marines Renouvelables[END_REF]. There are numerous technologies for tidal turbines and most of them can be classified into two categories: the horizontal axis and the vertical axis. This work is realized in collaboration with Hydroquest (Figure 1), a vertical axis tidal turbines manufacturer. Thus, this paper investigates the behavior of vertical axis tidal turbines. Most tidal sites are characterized by strong turbulence intensities [START_REF] Clark | Turbulence in Marine Environments (TiME): A framework for understanding turbulence and its effects on tidal devices[END_REF]. Several test campaigns have already been performed, such as in the sound of Islay where the turbulence intensity was found to be between 9.5% and 10.3% [START_REF] Milne | Characteristics of the turbulence in the flow at a tidal stream power site[END_REF] or in the Grand Passage where the turbulence intensity was found to be between 10% and 20% [START_REF] Hay | Turbulence and Bottom Stress in Grand Passage and Minas Passage[END_REF]. The presence of turbulence in most tidal sites has led companies and academics to study the effects of ambient turbulence over tidal turbines. Mycek et al. [START_REF] Mycek | Experimental study of the turbulence intensity effects on marine current turbines behaviour. Part I: One single turbine[END_REF] carried out such study for a horizontal axis tidal turbine model. They showed that going from 3% to 15% turbulence intensity decreases the power coefficient by a maximum of 10%. It also reduces the wake recovery by more than 5D. Tidal turbines will most likely be installed into farms. The results from Mycek et al. show that interactions between turbulence and tidal turbines will directly affect the layout and performances of farms. Predicting those interactions is thus crucial for reducing the cost of tidal farms. Several studies of vertical axis tidal turbines have been conducted during the past few years. Recently, Ouro et al. [START_REF] Ouro | Three-dimensionality of the wake recovery behind a vertical axis turbine[END_REF] carried out an experimental study of a vertical axis tidal turbine model in a turbulent flow with approximately 15% turbulence intensity for the streamwise velocity. They have highlighted different regions of the wake and the persistence of the turbine footprint up to 14D. Hochman et al. [START_REF] Hochman | The effect of inflow conditions on vertical axis wind turbine wake structure and performance[END_REF] studied the impact of the presence of an atmospheric boundary layer on vertical axis wind turbine model. They showed that there is little impact of the atmospheric boundary layer on the turbine. However, little is known about the impact of several turbulent inflow conditions over a single real size tidal turbine.

Vertical axis tidal turbines are home to complex phenomena such as dynamic stall. During a revolution of the rotor, the blade of a vertical axis turbine have an oscillating angle of attack that passes through a maximum positive value and a minimum negative value. At a low enough tip-speed ratio, the angle of attack can be superior to the static stall angle. The main effect of the dynamic stall is to delay the stall, thus increasing the performance of the turbine. It is therefore important to consider the dynamic stall when studying vertical axis tidal turbine. The downstream half of the blade trajectory is located in the wake of the upstream half. Thus, there is a strong asymmetry between the upstream half and the downstream half. Brochier et al. [START_REF] Brochier | Water channel experiments of dynamic stall on Darrieus wind turbine blades[END_REF] showed how vortex shedding occurs during the rotation of the blades and how it affects the right/left symmetry. Ouro et al. [START_REF] Ouro | An immersed boundary-based large-eddy simulation approach to predict the performance of vertical axis tidal turbines[END_REF] proposed a Large Eddy Simulation (LES) geometry resolved study of a vertical axis tidal turbine. They showed the vortex shedding strong dependence on tip-speed ratio and how it affects the hydrodynamic of the turbine.

Experimental studies are limited to turbine models and thus suffer from scale effects. In situ measurement are still expensive and difficult to acquire. Blade-resolved Computational Fluid Dynamics (CFD) approaches are appropriate to study vertical axis tidal turbines [START_REF] Ouro | An immersed boundary-based large-eddy simulation approach to predict the performance of vertical axis tidal turbines[END_REF]. Because the hydrodynamic of vertical axis tidal turbines is highly unsteady, LES are well adapted here [START_REF] Ouro | An immersed boundary-based large-eddy simulation approach to predict the performance of vertical axis tidal turbines[END_REF]. LES computational cost is still prohibitive for real size tidal turbine [START_REF] Ahmed | Fluctuating loads on a tidal turbine due to velocity shear and turbulence: Comparison of CFD with fields data[END_REF] and simplified approaches are preferred solutions. Shamsoddin and Porté-Agel [START_REF] Shamsoddin | Large Eddy Simulation of Vertical Axis Wind Turbine Wakes[END_REF] showed that Actuator Line Model (ALM) LES is more appropriate than Actuator Cylinder LES for modeling vertical axis turbine. Liu et al. [START_REF] Liu | An actuator line-Immersed boundary method for simulation of multiple tidal turbines[END_REF] modeled a horizontal axis tidal turbine model and showed that ALM-LES is well suited for modeling tidal turbine in different turbulent conditions.

This paper proposes the study of a ducted three-bladed H-type Darrieus contra-rotating vertical axis tidal turbine. This study was carried out with a Lattice Boltzmann Method (LBM) ALM LES approach. The LBM is a highly parallelizable unsteady CFD method [START_REF] Guo | Lattice Boltzmann Method and Its Application in Engineering[END_REF]. Since it has a low numerical dissipation [START_REF] Marié | Comparison between lattice Boltzmann Method and Navier-Stokes high order schemes for computational aeroacoustics[END_REF], it is well suited for propagating the wake of tidal turbines. The actuator line model used considers dynamic stall effect. Four upstream turbulence intensities were simulated. The chosen turbulence intensities are 0%, 4%, 8% and 12%. The main goal of this study was to highlight the effects of several turbulence rates on a single vertical axis tidal turbine, thus providing a better understanding to tidal turbines manufacturers. This work is divided into three sections, after this Introduction. Section 2 gives elements on the theoretical background of the approach. Since this approach has not been validated yet, validation elements are also presented in this section. Section 3 presents the results of ALM-LBM simulations of the vertical axis turbine. Section 4 gives the conclusions and prospects of this work.

Materials and Methods

This section is divided into four subsections. Section 2.1 succinctly presents the LBM models used. Section 2.2 introduces the implemented ALM model. Then, the approach is validated with an experimental case in Section 2.3. Section 2.4 is dedicated to the real-size turbine simulation set-up.

The Lattice Boltzmann Method

The LBM is an unsteady CFD method based on the Lattice Gas Automata and the Boltzmann equation [START_REF] Guo | Lattice Boltzmann Method and Its Application in Engineering[END_REF]. The LBM solver used here is the open-source library Palabos [15], developed at the University of Geneva. LBM describes the fluid from a molecular point of view through distribution functions f . Distribution functions represent the density of molecules at a certain location in space with a given velocity. Thus, it is a continuous function of space, time and velocity. Space is discretized into a Cartesian grid, or lattice, with a grid spacing ∆x. Cartesian grids in Palabos can be divided into several levels, by a factor 2 in mesh size and time step between two adjacent levels. This multi-block approach was described by Lagrava et al. [START_REF] Lagrava | Advances in multi-domain lattice Boltzmann grid refinement[END_REF].

Velocity space in Palabos is discretized by defining a finite set of directions of propagation. At each lattice node lies n distribution functions f i , each representing the density of molecules going in a direction. The set used in this article is the D3Q27 set with 27 directions of propagation, where 3 stands for three dimensions. Time in LBM is discretized into constant time steps ∆t. The fluid behavior is resolved into two steps, the collision step and the streaming step. The collision step is calculated with Equation (1):

f i (x, t) = Ω( f i (x, t)), (1)
where Ω is the collision operator and f i are the post-collision distribution functions. The collision operator represents the collision between molecules and is the core of the method. There are many collision operators available in Palabos and even more in the literature. The collision operator used in this paper is a regularized BGK collision operator. It was described by Malaspinas [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: Recursivity and regularization[END_REF] and was chosen for its stability [START_REF] Malaspinas | Increasing stability and accuracy of the lattice Boltzmann scheme: Recursivity and regularization[END_REF]. The collision operator used allows fluid forcing. Indeed, an external force is needed for implementing the actuator line model. The streaming step consists in the propagation of distribution functions onto the neighboring nodes pointed by the discrete directions of propagation. It is calculated with Equation (2):

f i (x + c i ∆t, t + ∆t) = f i (x, t)), (2)
where c i are the directions of propagation and ∆t is the time step. To reduce the mesh size, a LES model is used. It is a static Smagorinsky model and was described by Malaspinas and Sagaut [START_REF] Malaspinas | Consistent subgrid scale modelling for lattice Boltzmann methods[END_REF]. Following the work of Shamsoddin and Porté-Agel [START_REF] Shamsoddin | Large Eddy Simulation of Vertical Axis Wind Turbine Wakes[END_REF], we chose C s = 0.1. All the boundary conditions used are based on a second-order accurate finite-difference scheme described by Skordos [START_REF] Skordos | Initial and boundary conditions for the lattice Boltzmann method[END_REF] and Latt et al. [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF]. It was chosen because of its stability at high Reynolds number [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF]. For generating the inflow turbulence, a Divergence Free Synthetic Eddy Method (DFSEM) [START_REF] Polleto | A new divergence free synthetic eddy method for the reproduction of inlet flow conditions for LES[END_REF] was implemented by Palabos, Grondeau et al. [START_REF] Grondeau | Numerical simulation of a tidal turbine model in a turbulent flow with the Lattice Boltzmann Method[END_REF]. To increase the numerical stability of simulations, viscous damping zones are placed on the side boundaries and on the outlet.

Actuator Line Model

Since it was not present in Palabos, an actuator line model was implemented. The chosen model is a modified MIT dynamic stall model, which was described in detail by Shamsoddin and Porté-Agel [START_REF] Shamsoddin | A large-eddy simulation of vertical axis wind turbine wakes in the atmospheric boundary layer[END_REF]. This latter approach has been slightly modified. In their article, the averaged angle of attack α f it is generated with an eighth-order polynomial. In this study, it is generated with a cubic spline in order to better capture the evolution of the angle of attack. Actuator Line Models, introduced by Shen and Sørensen [START_REF] Sørensen | Numerical Modling of Wind Turbine Wakes[END_REF], represent the blades with a force f AL computed from the airfoil tabulated C L and C D . The dynamic stall model acts on those two coefficients and has a non-negligible influence [START_REF] Shamsoddin | Large Eddy Simulation of Vertical Axis Wind Turbine Wakes[END_REF]. The blades of the turbine are simplified into lines, which are discretized into a number of elements (Figure 2). The force f AL is applied to the fluid at the location of these elements and is calculated with Equation (3).

f AL ρ = (c∆l)( U 2 rel 2 )(C L DS e L + C D DS e D ), (3)
where c is the chord of the blade, ∆l is the spacing between elements and U rel is the norm of the relative velocity of the blade. Spacing ∆l is equal to the mesh size at the turbine location. C L DS and C D DS are the lift and drag coefficients, respectively, modified by the dynamic stall algorithm described in [START_REF] Shamsoddin | A large-eddy simulation of vertical axis wind turbine wakes in the atmospheric boundary layer[END_REF]. e L is the lift unit vector and e D is the drag unit vector. The relative velocity is calculated with Equation ( 4):

U rel = (U f luid .N)N + (U f luid .T -U ω )T, (4)
where U ω = ω R R, with R the radius of the rotor and ω R its rotational speed. U f luid is the fluid velocity calculated at the element location. The vector N is normal to the blade displacement and the vector T is tangent to the blade displacement. The angle of attack α is calculated with Equation ( 5):

α = tan -1 U f luid .N U f luid .T -U ω . ( 5
)
The actuator line moves at the rotational speed of the rotor ω R at each lattice iteration. Shamsoddin et al. [START_REF] Shamsoddin | Large Eddy Simulation of Vertical Axis Wind Turbine Wakes[END_REF] applied the force on a single node of the mesh. To avoid numerical instabilities, the force is distributed around each element of the actuator line. The distribution kernel is η( , d) and is given in Equation ( 6):

η( , d) = 1 3 π 3/2 exp - d 2 , ( 6
)
where d is the distance between a node of the lattice and an actuator line element. Studies about the choice of the parameter already exist. Some studies suggest choosing = 4 ∆x [START_REF] Bachant | Actuator line modeling of vertical-axis turbines[END_REF][START_REF] Mendoza | Validation of an actuator line model coupled to a dynamic stall model for pitching motions characteristic to vertical axis turbine[END_REF], where ∆x is the mesh size at the turbine location. In Palabos, the domain is divided into cuboids of n × n × n nodes; here, n = 14. Those cuboids are then distributed as equally as possible between the processors.

To do non-local operations, the nodes surrounding a cuboid need to be loaded with it, thus degrading the efficiency. To keep the efficiency of Palabos and stay as local as possible, a compromised between the methods of Bachant et al. [START_REF] Bachant | Actuator line modeling of vertical-axis turbines[END_REF] and Shamsoddin and Porté-Agel [START_REF] Shamsoddin | Large Eddy Simulation of Vertical Axis Wind Turbine Wakes[END_REF] is made and we chose = 2 ∆x.

The same value was chosen by Liu et al. [START_REF] Liu | An actuator line-Immersed boundary method for simulation of multiple tidal turbines[END_REF] for modeling a horizontal axis tidal turbine. The force applied to lattice nodes is calculated with Equation ( 7):

f AL (x) = η( , d) f AL . (7)

Accuracy of the Model

The code was validated with the experimental results from Rolin and Porté-Agel [START_REF] Rolin | Exprimental investigation of vertical-axis wind-turbine wakes in boundary layer flow[END_REF]. The case consists in a three-bladed vertical axis wind turbine placed in a turbulent boundary layer. The rotor characteristics are given in Table 1 and the turbine is illustrated in Figure 3. The blades are straight and their profile is a NACA0018. The wind tunnel is 28.00 m long, 2.57 m wide and 2.00 m tall. The turbine is located 20.00 m downstream of the tunnel contraction such that the turbulent boundary layer is fully developed [START_REF] Rolin | Exprimental investigation of vertical-axis wind-turbine wakes in boundary layer flow[END_REF]. The turbine half-height is located 0.13 m above the ground. The reference velocity is the axial velocity at the turbine mid-height, U H/2 = 9.43 m.s -1 . The turbine was modeled with the ALM dynamic-stall model presented in Section 2.2. Tabulated C L and C D were taken from Kumar et al. [START_REF] Kumar | Low Reynolds Number Vertical Axis Wind Turbine for Mars[END_REF] for low angles of attack (inferior to 25 • ) and from Sheldahl and Klimas [START_REF] Sheldahl | Aerodynamics Characteristics of Seven Symmetrical Airfoil Sections through 180-Degree Angle of Attack for Use in Aerodynamic Analysis of Vertical Axis Wind Turbine[END_REF] for high angles of attack [START_REF] Shamsoddin | Large Eddy Simulation of Vertical Axis Wind Turbine Wakes[END_REF]. The mast of the turbine has a radius of approximately 4 × 10 -3 m and it is also modeled with a force. Its drag coefficient is 1. The H-bars are also modeled with a force. Their length is approximately 7 × 10 -3 m, their lift is neglected, and they have a drag of 0.05, estimated with a 2D kω URANS simulation. The upstream velocity profiles are generated with a DFSEM boundary condition. The DFSEM, based on the SEM of Jarrin et al. [START_REF] Jarrin | Reconstruction of turbulent fluctuations for hybrid RANS/LES simulations using a Synthetic-Eddy-Method[END_REF], was introduced by Poletto et al. [START_REF] Polleto | A new divergence free synthetic eddy method for the reproduction of inlet flow conditions for LES[END_REF]. It generates a random upstream turbulence that respects a prescribed Reynolds tensor. Input parameters of the DFSEM are the average axial velocity profile, the average Reynolds tensor profile, the size of the DFSEM eddies and the number of DFSEM eddies. Upstream average axial velocity and fluctuations were extracted from Rolin and Porté-Agel [START_REF] Rolin | Exprimental investigation of vertical-axis wind-turbine wakes in boundary layer flow[END_REF]. The Reynolds tensor is symmetric and all its diagonal components are equal to the axial velocity fluctuations. Eddies are spherical and their diameter is 0.15 m. The number of eddies is set to have an eddy density of 2.0. The eddy density is defined as the ratio between the volume of eddies and the volume of the DFSEM domain [START_REF] Polleto | A new divergence free synthetic eddy method for the reproduction of inlet flow conditions for LES[END_REF]. The top and side boundary conditions are velocity-imposed boundary conditions set to match the prescribed inflow velocity profile. The outlet is a pressure-imposed boundary condition. The bottom boundary condition is a velocity-imposed boundary condition modified with a wall law, based on the formulation of Spalding [START_REF] Spalding | A single formula for the law of the Wall[END_REF].

Computational domain dimensions were set to match the wind tunnel dimensions, except for the length. The computational domain is 6 m long and the turbine is located 12D downstream of the inlet. Because the distribution kernel parameter varies with the mesh size ∆x, three mesh sizes were tested to evaluate the influence of both mesh size and distribution kernel. The LBM is a weakly compressible method. The time step was set to keep the Mach number identical in all cases, Ma max ≈ 0.08. Table 2 summarizes the various meshes, visible in Figure 4. Simulations were carried out at CRIANN [32] (France) on Broadwell CPU cores (2.4 GHz).

U De f = 1-< U x > /U H/2 (8)
U * x =< U x > /U H/2 (9)
TKE + = TKE W -TKE 0 (10)
Three wake quantities were investigated. The average axial velocity deficit is defined in Equation ( 8), where < U x > is the time-averaged axial velocity. The dimensionless average axial velocity is calculated with Equation ( 9). The dimensionless turbulent kinetic energy (TKE) is defined in Equation [START_REF] Ahmed | Fluctuating loads on a tidal turbine due to velocity shear and turbulence: Comparison of CFD with fields data[END_REF], where TKE W is the TKE calculated in the wake and TKE 0 is the TKE calculated upstream of the turbine. Both TKE are divided by U 2 H/2 . Figure 5 shows the location of the turbine and the plane where the statistics are calculated. The convergence lasted 70 rotations and data were averaged over 60 rotations of the turbine. An additional plane was extracted 2.0D upstream of the turbine to evaluate differences between the experimental flow and the numerical flow (Figure 6). LBM average axial velocity profiles are close to experimental profiles. Differences are observed in the ground vicinity. Indeed, the wall law does not allow predicting the proper flow behavior in the viscous sub-layer. Similarly, high differences are observed on the TKE in this area. The mesh resolution does not allow for a proper turbulence production. The turbulence generated by the DFSEM close to the bottom is not propagated on a long enough distance and, in lack of production, it is dissipated. Above the flume bottom, TKE profiles from the ALM-LBM simulations are close to the experimental one.

The small TKE peak observed at z/D = 1.8 comes from the influence of the mesh interface between two levels of the mesh. Figure 7 shows the average axial velocity deficit in the wake of the turbine and Figure 8 the average axial velocity. The maximum velocity deficit at 2.0D is different for the three meshes. The velocity deficit for Mesh No. 3 is too small compared to the experimental one. This difference tends to disappear after 4.0D and the maximum velocity deficit is correctly predicted. The only difference that persist in the far wake is the asymmetry of the profile in Figure 7. It is more pronounced on the experimental profile than on the simulation profiles. Figures 9 and10 show the average TKE in the wake of the turbine. The asymmetry difference can be observed in Figure 9. The ALM model is functioning with a wide angle of attack range (between -60 • and 60 • ) because of the low tip-speed ratio. This could induce some discrepancies. Figure 11 shows the contour of TKE + in a plane 2.0D downstream of the turbine. In-plane velocity vectors are also plotted. Two contra-rotating vortices in the negative y side are identified in both the experimental and simulation results. Further away in the wake, 6.0D downstream of the turbine, the bottom vortex is no longer visible in the experimental plane (Figure 12) while it is still there in the simulation. In the experiment, the bottom vortex vanishes because of the turbulent boundary layer. In the simulations, because the turbulence is not energetic enough close to the flume bottom, the contra-rotating vortices pair is conserved. This difference in wake structure could also induce the discrepancies observed on the wake asymmetry. Vertical profiles plotted in Figures 8 and10 are relatively close to experimental profiles. Only differences close to the bottom persist. Errors of the LBM-LES simulations are given in Table 3. Mesh No. 2 gives the best results overall. It is surprising to find out that it is not Mesh No. 3. In Section 2.2, it is explained that the parameter depends on the mesh size and is equal to = 2∆x. Therefore, the distribution volume for Mesh No. 3 is rather small compared to the blade chord. This could explain the differences for Mesh No. 3. It would be interesting to increase the distribution parameter as the mesh size decreases, but for reasons of efficiency, it is not possible to do so at this time. Mesh No. 2 gives relatively good wake prediction both for shapes and quantities. The configuration with ∆x = D/30 was kept for modeling the real-size tidal turbine. Most of the differences observed between the experiment and the simulations probably come from the low mesh resolution close to the bottom of the flume. Resolving the flow in this area would cost a lot of computational resources and is not adapted to simplified approaches such as ALM. Since the real-size turbine will not be as close to the viscous sublayer as the turbine from Rolin et al. [START_REF] Rolin | Exprimental investigation of vertical-axis wind-turbine wakes in boundary layer flow[END_REF], its influence should not be as significant.

U * x =< U x > /U H/2 z * = z/D U * x 0 2 4 •10 -2 T KE 0 = T KE/U 2 H/2 2
where U x δ is the axial velocity at height δ and z is the distance from the sea-bed. For the simulations presented afterwards, δ is equal to 10.0 m and U x δ = 1.475 m.s -1 . The parameter n depends on the flow. Here, n = 6.5 was chosen, relatively to the velocity profile measured in the Alderney Race from Thiébaut et al. [START_REF] Thiébaut | Merging velocity measurments and modeling to improve understanding of tidal stream resource in Alderney race[END_REF] and shown in Figure 14. The tip-speed ratio is calculated with λ = ωR/U x δ = 1.5. The Reynolds tensor was chosen to have a constant turbulence intensity from the sea-bed to the surface. The prescribed Reynolds tensor is diagonal and its components are equal. It is constructed to have I x = I y = I z = [0, 4, 8, 12] %. The eddies are spherical and their radius is equal to L = 8.0 m. The bottom boundary condition is a velocity-imposed boundary condition modified with the law of the wall of Spalding [START_REF] Spalding | A single formula for the law of the Wall[END_REF]. The top and side boundary conditions are velocity-imposed boundary conditions that match the velocity profile given by Equation [START_REF] Shamsoddin | Large Eddy Simulation of Vertical Axis Wind Turbine Wakes[END_REF]. The outlet is a pressure-imposed boundary condition. The depth was chosen close to the depth of the Alderney race tidal site [START_REF] Thiébaut | Merging velocity measurments and modeling to improve understanding of tidal stream resource in Alderney race[END_REF]: z = 43.12 m. The domain is 460 m long and 200 m wide. The turbine is located 4D eq downstream of the inlet. The numerical domain is divided into three sub-domains. The finest domain, where the turbine is located, goes up to 3.2D eq downstream of the turbine. The intermediate domain goes from 3.2D eq up to 10.5D eq and the coarsest domain is for the remaining domain. Figure 15 is a cross-sectional view of the mesh and Figure 16 is a three-dimensional view of the computational domain. The finest mesh was chosen according to Section 2.3: ∆x = D/30 = 0.257 m. There are 12.6 × 10 6 LBM nodes in the mesh and it takes 48.5 h.CPU to compute a single revolution.

Results

This section gives the results of the simulations of the turbine presented in Section 2.4. It is divided into two subsections. Section 3.1 is about the influence of inflow conditions over averaged wake quantities. Section 3.2 shows how inflow turbulence influences the shape of the wake.

Wake Recovery

Averages were computed over 60 rotations, after a convergence of 40 rotations. The results were made dimensionless for confidentiality reasons. The average axial velocity < U x > was divided by U re f . Figure 17 shows the horizontal average axial velocity profiles in the wake of the turbine. Profiles from the bottom rotors and from the top rotors are different, independently of the turbulence rate. The asymmetry according to the y direction is stronger for the bottom profiles than for the top profiles. Because of the upstream velocity profile, the fluid axial velocity at the bottom rotors height is lower than the velocity at the top rotors height. Thus, the tip-speed ratio is higher for the bottom rotors. The asymmetry is counter-intuitive since it should weaken as the tip-speed ratio increases. Differences between the bottom and top profiles tend to be reduced by the ambient turbulence. The higher the turbulent rate is, the lower the asymmetry is. The asymmetry diminishes significantly after 1.6D eq and even disappears for the higher turbulence rate. The same can be observed on the vertical axial velocity profiles, as visible in Figure 18. In the close wake, the velocity deficit is greatly influenced by the ambient turbulence. In the far wake, the upstream turbulence intensity seems to have little influence over the horizontal average axial velocity profiles. The turbine footprint is still present 5.6D eq downstream of the turbine and the maximum velocity deficit is lower than 20%. However, on the vertical axial velocity profiles, the turbulence seems to affect the velocity deficit in the far wake. The higher the upstream turbulence is, the lower the velocity deficit is.

Figure 19 shows the average TKE profiles in the wake of the turbine. In the close wake, 0.8D eq downstream of the turbine, four peaks can be observed on the I = 0% and I = 4% profiles. The two peaks from the positive y side come from the positive y rotors and the two negative y peaks come from the negative y rotors. On the I = 8% and I = 12% turbulence intensity profiles, only three peaks are visible: two peripheral peaks coming from the negative y rotors and positive y rotors and one big central peak. This highlight the influence of ambient turbulence over the merging of the wake of the positive y rotors and the negative y rotors. The TKE levels are greatly influenced by the upstream turbulence. A higher upstream turbulence increases the TKE in the wake of the turbine. After 4.0D eq downstream of the turbine, the TKE levels in the wake are close to the TKE of the ambient flow. However, 5.6D eq downstream of the turbine, the footprint is still visible for the I = 0% profile and the TKE is not zero. x/D eq = 5.6 Profiles of average wake quantities have been described. The upstream turbulence influence have been highlighted. However, a more comprehensive approach is needed to better understand the influence of ambient turbulence over the wake of the turbine. Such approach is carried out in the next section.

Wake Structures

Contours of the velocity deficit in the wake of the turbine are visible 1.2D eq downstream of the turbine in Figure 20, 3.2D eq downstream in Figure 21 and5.2D eq downstream in Figure 22. Velocity vectors are also projected onto the extracted planes. The influence of ambient turbulence on the shape of the wake is visible at each downstream distance. In the close wake, 1.2D eq downstream of the turbine, two distinct wakes are observed on the I = 0% figure while a single wake is observed on the other profiles. This confirms the influence of turbulence over wakes merging, a phenomenon already highlighted in Section 3.1. Two vortices can be observed at the location of the top corners of the turbine. In Figure 20, those vortices are observed independently of the turbulence rate. However, their intensity seems to be diminished by the ambient turbulence. They appear to be giving the wake a T shape.

I = 0 % I = 4 % I = 8 % I = 12 %
Figures 23 and24 show instantaneous iso-contours of the Lambda2 criterion in the wake of the turbine for, respectively, I = 0% and I = 12%. Planes of the average < U y > velocity are also represented. In both cases, the vortices are visible on the < U y > planes at x/D eq = 1.2. At x/D eq = 4.2, the vortices are visible for I = 0% but much less obvious for I = 12%. On the instantaneous iso-contours, those vortices are visible as sets of smaller vortices. For I = 0%, two sets are clearly observable up to x/D eq = 4.2. For I = 12%, the two sets are only observable up to x/D eq = 1.2. The ambient turbulence seems to affect the coherence of the vortices sets. This was seen with the velocity contours and vectors.

At the moment, it is not possible to know whether they are generated by the stator or by the rotor. Further in the wake, the ambient turbulence affects these vortices and causes them to lose their consistency. The shape of the wake is thus influenced and the T shape is smoothed into a square shape. In Figure 22, 5.2D eq downstream of the turbine, the T shape is observable only on the I = 0% figure.

As shown in Section 3.1, a strong vertical asymmetry was observed in the close wake. It is possible that this asymmetry is caused by these two vortices. Indeed, at the location of the bottom rotors, the radial velocity is directed toward the turbine center. This phenomenon crushes the wake of the bottom rotors toward the center of the turbine. Since the vortices are weaker for high turbulence rates, the asymmetry is greatly reduced. This is observed in Figure 17.

Conclusions

The implementation of an ALM in a LBM solver was realized. The ALM-LBM-LES approach was validated with the experimental case of Rolin and Porté-Agel [START_REF] Rolin | Exprimental investigation of vertical-axis wind-turbine wakes in boundary layer flow[END_REF]. A mesh sensitivity study was performed and the mesh ∆x = D/30 was found to have the highest accuracy. Some discrepancies with the experiment were observed and are most likely due to the low mesh resolution in the viscous sublayer. It is to be noted that there is no tip-loss correction in the ALM-LBM-LES model. Yu et al. [START_REF] Yu | Study on actuator line modeling of two NREL 5-MW wind turbine wakes[END_REF] showed that adding a tip-loss correction when modeling a horizontal axis wind turbine can reduce the error of the predicted power. It could be an interesting addition to the ALM-LBM-LES model. Nevertheless, the TKE and velocity deficit are relatively well predicted. The shape of the wake is also correctly computed and the two contra-rotating vortices pairs are observed in the close wake, 2.0D downstream of the turbine.

The ALM-LBM-LES computational cost is a lot lower than LBM-LES blade-resolved approaches. A LBM-LES blade-resolved study of a horizontal axis tidal turbine model was performed by Grondeau et al. [START_REF] Grondeau | Modelling the wake of a tidal turbine with upstream turbulence: LBM-LES versus Navier-Stokes LES[END_REF]. The computational time required to compute a single revolution of the turbine is 450 h.CPU, whereas it is 14 h.CPU for the tidal turbine model simulated here. The same horizontal axis tidal turbine model was studied by Elie et al. [START_REF] Elie | Simulation of horizontal axis tidal turbine wake using a Weakly-Compressible Cartesian Hydrodynamic solver with local mesh refinement[END_REF] with an actuator-disc LES approach. Their code uses an explicit weakly-compressible Navier-Stokes solver. For a mesh size of ∆x = D/30, it takes 7.2 h.CPU to compute a revolution of the turbine. It is half the time needed here. However, they only meshed the turbine area with a mesh size of ∆x = D/30, thus the number of elements is 7.6 × 10 5 compared to 5.2 × 10 6 here where the wake is also refined at ∆x = D/30.

The validated mesh configuration and ALM-LBM-LES model were then used to compute the wake of a real-size tidal turbine. The influence of the turbulence over the wake was assessed. It was observed that the close wake is strongly modified by the upstream turbulence. The influence of ambient turbulence on the far wake is weaker but can still be observed on velocity contours. The presence of two vortices at the top corner of the turbine was highlighted. It was showed that they may have an influence on the vertical asymmetry of the wake. Their intensity and consistency is weakened by the ambient turbulence.

The influence of ambient turbulence on the far wake is limited. Even if the top-corner vortices lose their intensity, they are still observable with I = 12%, 140.0 m downstream of the turbine. Given the way they affect the shape of the wake, it is possible that they contribute to maintaining the velocity deficit between them. It would then be interesting to use the same numerical configuration to investigate the behavior of the turbine with other turbine configurations, for example with an opposite direction of rotation and without the stator.

After 2.0D eq , the velocity deficit is relatively low. From the horizontal profiles of the average axial velocity, it is estimated to be around 10% at x/D eq = 2.4 for I = 12%. Since it has been shown that the ALM-LBM-LES gives relatively good wake prediction in the close wake, it would be interesting to study a farm with longitudinal spacing starting from [2.0, 2.5]D eq . Indeed, one way to optimize tidal farms is to find the smallest distance possible between turbines. A four-turbine farm modeled with the ALM-LBM-LES is estimated to use around [200,300] h.CPU for each revolution of the turbines. Considering the computational resources we have, this would allow us to model several farm configurations.

Figure 1 .

 1 Figure 1. Hydroquest/CMN tidal turbine to be installed at Paimpol-Bréhat test site (France).

Figure 2 .

 2 Figure 2. Actuator line scheme with: U ∞ the upstream velocity, Ω the rotational speed, θ the angular position, C the center of rotation, U ω the tangential velocity, U rel the relative velocity, α the angle of attack, L the lift force, and D the drag force.

Figure 3 .

 3 Figure 3. Vertical axis wind turbine tested in Rolin and Porté-Agel [27].

Table 2 .

 2 Configurations for the ALM-LBM mesh sensitivity study. h.CPU is wall-clock-time multiplied by the number of CPU used. Mesh D/∆x ∆x (m) ∆t (s) Number of Nodes h.CPU per rev. No. 1 15 11.16 × 10 -3 3.04 × 10 -5 1.5 × 10

Figure 4 .

 4 Figure 4. Meshes used for the ALM-LBM mesh sensitivity study. Meshes are represented in cubes of 14 × 14 × 14 nodes. Turbine is represented with the blue cylinder. Plane xy is at turbine mid-height and plane xz is at canal mid-width.

Figure 5 .

 5 Figure 5. Average axial velocity deficit at different locations in the wake of the turbine. Axis origin is the center of the turbine. Dimensions are scaled by the turbine diameter. The results are from ALM-LBM simulation with Mesh No. 2. U De f = 1-< U x > /U H/2 .

Figure 6 .

 6 Figure 6. Average axial velocity and TKE upstream of the tidal turbine[START_REF] Rolin | Exprimental investigation of vertical-axis wind-turbine wakes in boundary layer flow[END_REF]: experimental data and ALM-LBM results for the three meshes. Turbine footprint is drawn in horizontal black dashed-lines.

Figure 7 .Figure 8 .• 10 - 2 Tx = 4 Figure 9 .• 10 - 2 Tx = 4 Figure 10 .

 7810249102410 Figure 7. Average axial velocity deficit y profiles in the wake of the tidal turbine[START_REF] Rolin | Exprimental investigation of vertical-axis wind-turbine wakes in boundary layer flow[END_REF] at TSR = 1.1: experimental data and ALM-LBM results for the three meshes. Turbine footprint is drawn in horizontal black dashed-lines.

Figure 11 .

 11 Figure 11. Contour of average TKE normalized by the upstream velocity U H/2 , 2.0D downstream of the turbine: experimental data from [27] on the left and from the ALM-LBM simulation with Mesh No. 2 on the right. Average velocity vectors are projected onto the plane. Turbine footprint is delimited in black.

Figure 12 .

 12 Figure 12. Contour of average TKE normalized by the upstream velocity U H/2 , 6.0D downstream of the turbine: experimental data from [27] on the left and from the ALM-LBM simulation with Mesh No. 2 on the right. Average velocity vectors are projected onto the plane. Turbine footprint is delimited in black.

Figure 14 .

 14 Figure 14.Axial velocity profiles obtained from the power law with n = 6.5 and from ADCP measurements in the Alderney Race during ebb[START_REF] Thiébaut | Merging velocity measurments and modeling to improve understanding of tidal stream resource in Alderney race[END_REF].

Figure 15 .

 15 Figure 15. Cross-sectional views of the mesh. The mesh is represented in cubes of 14 × 14 × 14 nodes. Dimensions are scaled by the equivalent diameter D eq = 25.0 m. Axis origin is the center of the turbine. (a): side view at the turbine mid-width. (b): top view at the turbine mid-height.

Figure 16 .

 16 Figure 16. Illustration of the computational domain. Dimensions are scaled by the equivalent diameter D eq = 25.0 m. The mesh on the domain boundaries is represented in cubes of 14 × 14 × 14 nodes. Stator is in red and rotors are in blue.

6 Figure 17 . 6 Figure 18 .

 617618 Figure 17. Average axial velocity y profiles in the wake of the tidal turbine at TSR = 1.5. ALM-LBM results for the four turbulence rates. Black curves are extracted at the bottom rotors mid-height. Red curves are extracted at the top rotors mid-height. Turbine footprint is drawn in horizontal black dashed lines. Dimensions are scaled by the equivalent diameter D eq = 25.0 m.

• 10 - 2 T 2 ref

 1022 KE + = T KE/U

Figure 19 .

 19 Figure 19. Average TKE + y profiles in the wake of the tidal turbine at TSR = 1.5. Horizontal profiles extracted at the turbine mid-height. ALM-LBM results for the four turbulence rates. Turbine footprint is drawn in horizontal black dashed lines. Dimensions are scaled by the equivalent diameter D eq = 25.0 m.

Figure 20 .I = 0 % I = 4 %I = 8 % I = 12 %Figure 21 . 4 %I = 8 % I = 12 %Figure 22 .

 20481221481222 Figure 20. Contour of axial velocity deficit, 1.2D eq downstream of the turbine. Velocity vectors are projected onto the plane. Turbine footprint is delimited in black. ALM-LBM simulations with four turbulence rates. U De f = 1-< U x > /U re f .

Figure 23 .

 23 Figure 23. Average < U y > velocity planes at x/D eq = 1.2 and x/D eq = 4.2 in the wake of the turbine for an upstream turbulence rate of I = 0%. Dimensions are scaled by the equivalent diameter D eq = 25.0 m. Vortices are represented with instantaneous iso-contour of the Lambda2 criterion, Lambda2 = -1.0 × 10 -5 . x and y axis are centered at the turbine location.

Figure 24 .

 24 Figure 24. Average < U y > velocity planes at x/D eq = 1.2 and x/D eq = 4.2 in the wake of the turbine for an upstream turbulence rate of I = 12%. Dimensions are scaled by the equivalent diameter D eq = 25.0 m. Vortices are represented with instantaneous iso-contour of the Lambda2 criterion, Lambda2 = -1.0 × 10 -5 . x and y axis are centered at the turbine location.

Table 1 .

 1 Rotor configuration from Rolin and Porté-Agel[START_REF] Rolin | Exprimental investigation of vertical-axis wind-turbine wakes in boundary layer flow[END_REF].

	Quantity	Value	Unit
	Number of blades n	3	Rotor
	Radius R	0.083	m
	Blades height H	0.155	m
	Chord length c	0.03	m
	Solidity Sd	0.54	Tip
	Speed Ratio λ	1.1	Reynolds
	number (radius) 5.0 × 10 4

Table 3 .

 3 Root-mean-square of the error on the U * x and TKE + z profiles. Error is averaged between z/D = 0.316 and z/D = 1.25.

	Mesh U *
	No. 1	22.0%	48.3%	5.3%	37.6%	5.2%	24.4%
	No. 2	11.0%	18.4%	5.8%	14.4%	6.9%	16.0%
	No. 3	6.8%	19.6%	5.1%	37.1%	6.0%	43.2%

x at 4D TKE + at 4D U * x at 8D TKE + at 8D U *

x at 10D TKE + at 10D

Acknowledgments: This work was funded by l'Agence de l'Environnement et de la Maîtrise de l'Energie (ADEME) [38] in the frame of the OCEANQUEST project and the Région Normandie in the frame of the project SEMARIN. Computational resources were provided by the Centre Régional Informatique et d'Applications Numériques de Normandie (CRIANN) [32] and by the Manche county council. The authors would like to thank Hydroquest [39] for the data provided concerning the Hydroquest/CMN tidal turbine. The authors would like to thank the Palabos project [15], without whom this work would not exist. The authors would also like to thank the ALGLIB project [40] for the open-source interpolation algorithm.

Funding: This research was funded by Agence de l'Environnement et de la Maîtrise de l'Energie grant number 1682C0085 and Conseil Régional de Normandie grant number 17E01557.

Turbine Set-Up

The turbine studied here has the same design as the Hydroquest turbine in Figure 1. The turbine has four rotors. The top rotors are shifted by 60 degrees compared to the bottom rotors. The left and right rotors are contra-rotating. Rotors are H-type Darrieus with straight blades. Blades chord and profile are confidential. Tabulated C L and C D are obtained with XFLR5 software with N crit = 9 and 300 panels. Table 4 summarizes the rotors characteristics. The tip-speed ratio λ is chosen close to the operating point of the turbine, λ = 1.5. axis and H-bars are modeled with actuator lines. As for the blades, forces from axis and H-bars are distributed around the actuator lines elements. The stator is modeled with an off-lattice approach, described in Guo et al. [START_REF] Guo | Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method[END_REF]. A no-slip boundary condition or a wall model would require an unaffordable fine mesh considering the Reynolds numbers at stake. A free-slip condition is used to model the interactions between rotors and stator. The off-lattice algorithm is set to impose a free-slip boundary condition on the stator walls. The bottom of the stator is located 5 m above the sea-bed. The inflow velocity profiles are generated with a DFSEM boundary condition. The average axial velocity profile is constructed with the power law given in Equation ( 11):

Author Contributions: Software, M.G., P.M. and E.P.; validation, M.G.; writing-original draft preparation, M.G.; writing-review and editing, M.G., S.G. and E.P.; supervision, S.G.; project administration, S.G.; funding acquisition, S.G.

Conflicts of Interest:

The authors declare no conflict of interest.

Abbreviations

The