Pascale Minet

Saoucene Mahfoudh

SERENA: an energy-efficient strategy to schedule nodes activity in wireless ad hoc and sensor networks

Keywords: energy efficiency, node activity scheduling, network lifetime, sleeping node, spatial reuse, coloring algorithm, slot assignment, wireless ad hoc networks, sensor networks, energy consumption

In wireless ad hoc and sensor networks, an analysis of the node energy consumption distribution shows that the largest part is due to the time spent in the idle state. This result is at the origin of SERENA, an algorithm to SchEdule RoutEr Nodes Activity. SERENA allows router nodes to sleep, while ensuring end-to-end communication in the wireless network. It is a localized and decentralized algorithm assigning time slots to nodes. Any node stays awake only during its slots and the slots assigned to its neighbors, it sleeps the remaining time. SERENA is based on distributed and localized two-hop coloring. The node's color is then mapped in time slot. Thus, each node is ensured to get at least one time slot, it also gets additional time slots proportionally to its traffic rate. Such a solution adapts to varying traffic rates and supports late node arrivals. A performance evaluation allows us to compare SERENA coloring algorithm with existing ones such as DLF, both in terms of number of colors and complexity. Simulation results show that SERENA enables us to maximize network lifetime while increasing the number of user messages delivered. We quantify the slot reuse and evaluate the impact of the frame size on network performance. We then study how to dimension buffers at the router nodes. Finally, we show how SERENA improves the node energy consumption distribution and maximizes the energy efficiency of wireless ad hoc and sensor networks.

SERENA: une stratégie efficace en énergie pour ordonnancer l'activité des noeuds dans les réseaux sans fil ad hoc et les réseaux de capteurs Résumé : Dans les réseaux sans fil ad hoc et les réseaux de capteurs, une analyse de la répartition de la consommation énergétique des noeuds montre que la plus grande part d'énergie est consommée dans l'état oisif. Ce résultat est à l'origine de SERENA, un algorithme pour ordonnancer l'activité des noeuds. SERENA permet aux noeuds routeurs de dormir tout en assurant la connectivité du réseau. C'est un algorithme localisé et décentralisé s'appuyant sur un algorithme d'assignation de slots temporels aux noeuds du réseau. Un noeud quelconque ne reste éveillé que durant ses slots et les slots attribués à ses voisins à un saut; il dort le reste du temps. SERENA utilise un algorithme distribué et localisé de coloriage à deux sauts. Un slot est ensuite associé à la couleur du noeud. Ainsi, chaque noeud est assuré de disposer d'au moins un slot. Il dispose de slots supplémentaires proportionnellement à son trafic. Cette solution s'adapte aux variations de trafic et supporte les arrivées tardives de noeuds. Une évaluation de performances permet de comparer SERENA à d'autres algorithmes de coloriage comme DLF, à la fois en termes de nombre de couleurs et complexité. Les résultats de simulation montrent que SERENA maximise la durée de vie du réseau tout en accroissant la quantité de données remises. Nous quantifions la réutilisation spatiale des slots et évaluons l'impact de la taille de la trame sur les performances du réseau. Finalement, nous montrons comment SERENA améliore la répartition de la consommation énergétique des noeuds et maximise l'efficacité énergétique des réseaux sans fil ad hoc et des réseaux de capteurs.

Mots-clés : efficacité énergétique, ordonnancement de l'activité des noeuds, durée de vie du réseau, noeud endormi, réutilisation spatiale, algorithme de coloriage, assignation de slots, réseaux mobiles ad hoc, réseaux de capteurs, consommation énergétique. With the increasing number of applications in many domains (such as detection of forest fire or seismic event, wild life protection, building and bridge monitoring, emergency rescue, target tracking, exploration mission in hostile environments and home monitoring), wireless ad hoc and sensor networks have a promising future. However, nodes in such networks can have a limited amount of energy. Moreover, this energy can be very expensive, difficult or even impossible to renew. That is why, energy efficient strategies are needed in order to maximize both network lifetime and the amount of data delivered.

The originality of this paper consists in identifying the highest energy consumption of a node and proposing a solution that improves network energy efficiency by scheduling node activity. This paper is organized as follows. In Section 2, we first analyze the node energy consumption distribution in wireless ad hoc and sensor networks and identify the main energy costs. Section 3 is a brief overview of the state of the art related to the scheduling of wireless nodes activity. Section 4 presents a new algorithm SERENA, SchEduling RoutEr Nodes Activity, whose originality consists in allowing router nodes to sleep. Time slots are assigned to nodes in a decentralized and localized way. A node is awake only during its slots and the slots granted to its neighbors, it sleeps the remaining time. In Section 5, we evaluate the impact of SERENA on the network lifetime and the amount of user data delivered. Section 6 deals with slot reuse and the impact of the size of the periodic frame on network performance. Section 7 shows how to dimension buffers at the router nodes. Finally, Section 8 quantifies the ability of SERENA to improve energy efficiency in wireless ad hoc and sensor networks and we conclude in Section 9.

State

Power value Transmitting P transmit = 1.3W Receiving P receive = 0.9W Idle P idle = 0.74W Sleeping P sleep = 0.047W Table 1: Power value in each radio state.

In Table 1, we report the reference values of P transmit , P receive , P idle and P sleep taken from a Lucent silver wavelan PC card [START_REF] Feeney | Investigating the energy consumption of a wireless network interface in an ad hoc networking environment[END_REF] for an IEEE 802.11b network. These values are used in the performance evaluation reported in this paper.

We now evaluate the energy consumption of wireless nodes, when a routing efficient strategy is used. This routing efficient strategy chosen consists in minimizing the energy consumed by the flow transmission in selecting the path with minimum energy dissipated, and avoiding depleted nodes.

For that, we will compute the energy dissipated by the end-to-end transmission of a flow packet, denoted cost(f low). Instead of using the number of hops between source and destination to select the best path (i.e.; every link has a cost of one), as done in OLSR [START_REF] Adjih | Optimized Link State Routing Protocol[END_REF], we will use cost(f low) as the criterion to choose the best path.

The energy cost of a flow on its path P is equal to:

cost(f low) = i∈sender(f low) cost transmission (i),
where i is a sender of f low on its path P .

When a transmitter transmits one packet to next hop, because of the shared nature of wireless medium, all neighbors of the source receive this packet even if it is intended to only one of them. Moreover, each node situated between transmitter range and interference range receives this packet but it cannot decode it. These two problems generate loss of energy. So to compute the energy dissipated by one transmission, we must take into account these losses as follows [START_REF] Allard | Evaluation of the energy consumption in MANET Adhoc-Now 2006[END_REF]:

cost transmission (i) = E transmit + n * E receive ,
where:

n is the number of non-sleeping nodes belonging to the interference area of the transmitter i, The Routing Efficient strategy, called RE is based on the OLSR routing protocol [START_REF] Adjih | Optimized Link State Routing Protocol[END_REF]. In order to avoid depleted nodes, RE modifies the MPR (Mul-tiPoint Relay) selection of OLSR to take into account the residual energy of nodes. These new MPRs are used to build the energy efficient routes, whereas the classical MPRs are used to optimize network flooding. For more details, the reader can refer to [START_REF] Mahfoudh | A comparative study of energy efficient routing strategies based on OLSR[END_REF].

Simulations have been performed for different wireless networks, where the network density (the average number of neighbors per node) is fixed to 10. Nodes whose number varies from 50 to 200, are randomly distributed in the network area. Network throughput is set to 2 Mbps. The initial energy of nodes is equal to 100 Joules. User traffic consists of 30 flows, with randomly chosen sources and destinations, and a throughput of 16 Kbps. The size of a message is 512 bytes. Messages of the routing protocol are not taken into account. Each result is the average of 5 simulation runs. Figure 1 illustrates the average on all network nodes of the energy dissipated in the different states defined previously. Notice that the Receiving state has been splitted into three substates:

Receive: when this node is the message destination, Overhearing: when this node is a one-hop neighbor of the transmitter but not the message destination, Interference: when this node is a two-hop neighbor of the transmitter. Furthermore, it clearly appears that the highest part of energy (about 50%) is surprinsingly dissipated in the idle state. The second highest energy cost is due to interferences (about 30%), even if the RE routing actively contributes to energy saving. The third cost is due to overhearing (about 10%). Finally, the energy dissipated in the Transmit and Receive states are small (about 3%). From these results, we can conclude that the most efficient strategy consists in allowing nodes to sleep. Solutions already exist for non router nodes: see for INRIA instance ZigBee [START_REF][END_REF]. Our purpose is to propose a new solution allowing router nodes to sleep and to evaluate its benefit on the node enenergy consumption distribution.

3 State of the art

Energy efficient strategies

The energy constrained nature of wireless nodes requires the use of energy efficient strategies to maximize network lifetime. We can classify these strategies in four categories: Some strategies find the optimum node transmission power that minimizes energy consumption, while keeping network connectivity, like for instance [START_REF] Ingelrest | Optimal Transmission Radius for Energy Efficient Broadcasting Protocols in Ad Hoc Networks[END_REF][START_REF] Cardei | Topology Control in Ad hoc Wireless Networks with Hitch-hiking[END_REF];

Other strategies reduce the volume of information transferred by:

aggregating information with the use of clusters, like [START_REF] Heinzelman | Energy-efficient communication protocol for wireless microsensor networks[END_REF][START_REF] Xia | Near-optimal node clustering in wireless sensor networks for environment monitoring[END_REF] or without, like [START_REF] Kalpakis | Maximum Lifetime Data Gathering and Aggregation in Wireless Sensor Networks[END_REF], [START_REF] Tang | Extending Network Lifetime for Precision-Constrained Data Aggregation in Wireless Sensor Networks[END_REF], [START_REF] Kulik | Adaptive Protocols for Information Dissemination in Wireless Sensor Networks[END_REF];

decreasing the frequency of information refreshment with distance, like [START_REF] Nguyen | Scalabilty of the OLSR protocol with the Fish Eye extension[END_REF];

avoiding to transfer information to uninterested nodes;

Other strategies focus on energy efficient routing in order to minimize the energy consumed by the end-to-end transmission of a packet, to avoid depleted nodes and reduce the number of unsuccessful transmissions, like [START_REF] Shah | Energy Aware Routing for Low Energy Ad Hoc Sensor Networks[END_REF], [START_REF] Ganesan | Highly-resilient,energyefficient multipath routing in wireless sensor networks[END_REF], [START_REF] Srinivas | Minimum Energy Disjoint Path Routing in Wireless Ad-Hoc Networks[END_REF], [START_REF] Senouci | Energy efficient routing in wireless ad hoc networks[END_REF], [START_REF] Kwon | Energy-Efficient Interference-Based Routing for Multi-hop Wireless Networks[END_REF], [START_REF] Hassanein | Reliable Energy Aware Routing In Wireless Sensor networks[END_REF], [START_REF] Shresta | Reception Awarness for Energy Conservation in Ad Hoc Networks[END_REF];

Finally, the last category allow nodes to sleep in order to spare energy, provided that the network and application functionnalities are still ensured, like [START_REF] Cardei | Energy-efficient target coverage in wireless sensor networks[END_REF], [START_REF] Carle | Energy-Efficient Area Monitoring for Sensor Networks[END_REF], [START_REF] Cardei | Improving wireless sensor network lifetime through power aware organization[END_REF].

In this paper, we focus on the last category, nodes activity scheduling, for the reasons given in the previous section.

Node activity scheduling

All solutions scheduling node activity determine time intervals during which a node must be awake and those during which it can sleep, knowing that the application and network functions must always be satisfied. In IEEE 802.15.4 [START_REF][END_REF], the MAC layer can operate in two modes: beacon-enabled mode in which all nodes can sleep. A superframe is used to indicate to each router when it can send data and when it can sleep. This mode is used in tree and star topologies. The second mode, non-beacon mode, is used in mesh topology. In this mode only non-routing nodes can sleep. Data intended to these nodes must be kept by their parent. In this paper, we propose a solution allowing any node, even router, to sleep. Some authors propose in [START_REF] Cardei | Energy-efficient target coverage in wireless sensor networks[END_REF] to extend network lifetime by dividing the network nodes in disjoint sets, such that each node set meets the network and application functions. These sets are activated successively, and at any time only the nodes of one set are active. All others nodes are in the sleep state. The problem consists in maximizing the number of disjoint sets. It has been shown NP-complete. The solution proposed is centralized. These authors have shown in [START_REF] Cardei | Improving wireless sensor network lifetime through power aware organization[END_REF] that network lifetime can be improved by allowing non-disjoint sets.

In [START_REF] Carle | Energy-Efficient Area Monitoring for Sensor Networks[END_REF], a distributed and localized solution is proposed. It consists in selecting a connected dominating set of sensor nodes (i.e. each node is either in this subset or is a neighbour of a node in this subset). Only the nodes of this set are active. All other nodes can change their state to sleep mode. Some solutions take advantage of spatial reuse to determine the time intervals dedicated to node activity. Indeed, during the same time intervall two transmitters can transmit simultaneously and succesfully if they do not interfer. Spatial reuse can be obtained by means of a coloring algorithm. We will now focus on such solutions.

Centralized and distributed one-hop coloring

One-hop graph coloring has received a lot of attention from researchers (see [START_REF] Brelaz | New methods to color the vertices of a graph[END_REF][START_REF] Peemoller | A correction to Brelaz's modification of Brown's coloring algorithm[END_REF][START_REF] Kubale | A better practical algorithm for distributed graph coloring[END_REF][START_REF] Hansen | Distributed largest-first algorithm for graph coloring[END_REF][START_REF] Kosowski | On greedy graph coloring in the distributed model[END_REF][START_REF] Kuhn | On the complexity of distributed graph coloring[END_REF]). One-hop graph coloring consists in coloring each vertex of the graph such that two adjacent vertices have not the same color and the number of colors used is minimum. This problem has been shown NP-complete in [START_REF] Garey | Computers and intractability: a guide to theory of NP-completeness[END_REF] for the general case, whereas graphs with maximum vertex degree less than four, and bipartite graphs can be colored in polynomial time.

The first one-hop graph coloring algorithms proposed were centralized (see [START_REF] Brelaz | New methods to color the vertices of a graph[END_REF][START_REF] Peemoller | A correction to Brelaz's modification of Brown's coloring algorithm[END_REF]). Among the greedy algorithms (i.e. no color backtracking), Dsatur, presented in [START_REF] Brelaz | New methods to color the vertices of a graph[END_REF], where the vertex with the highest number of already colored neighbor vertices is colored first, exhibits very good performances, even if it is not optimal. It is then followed by Largest First, where the node with the highest degree is colored first.

Distributed one-hop graph coloring algorithms also exist. Some of them resort to randomization to select the color as for instance [START_REF] Johansson | Simple distributed (∆ + 1) -coloring of graphs[END_REF][START_REF] Busch | Contention-free MAC protocols for wireless sensor networks[END_REF] and [START_REF] Finoccho | Experimental analysis of simple distributed vertex coloring algorithms[END_REF]. The color selected by a node can be used only if it does not conflict with the colors chosen by its neighbors. Other algorithms are strictly deterministic. The authors of [START_REF] Kosowski | On greedy graph coloring in the distributed model[END_REF] require that the results of the distributed algorithm and its centralized version be identical for any graph. This constraint is not required in the case of wireless ad hoc and sensor nodes. In [START_REF] Hansen | Distributed largest-first algorithm for graph coloring[END_REF], it is shown that in some network configurations, Distributed Largest First uses more colors than Largest First and the reverse is true in some other network configurations. This is because with Distributed Largest First, a node is allowed to keep its selected color, not only if it has the highest degree among its neighbors (as Largest First does), but also if there is no color conflict among its neighbors. Another approach consists in finding maximum independent sets and then coloring these sets independently, as in [START_REF] Herman | A distributed TDMA slot assiognment algorithm for wireless sensor networks[END_REF] and [START_REF] Kosowski | On greedy graph coloring in the distributed model[END_REF], because both problems are related [START_REF] Kuhn | On the complexity of distributed graph coloring[END_REF].

INRIA

Complexity of a distributed coloring algorithm

The efficiency of a distributed coloring algorithm, [START_REF] Hansen | Distributed largest-first algorithm for graph coloring[END_REF][START_REF] Kuhn | On the complexity of distributed graph coloring[END_REF], can be evaluated by: the number of colors needed to color a graph G: closer this number to the chromatic number of G, more efficient the algorithm. its time complexity, expressed in the case of a distributed algorithm, by the maximum number of rounds needed to color each node. A round is defined such that every node can:

send a message to all its one-hop neighbors, -receive the messages sent by them, -perform some local computation based on the information contained in the received messages.

Let n be the number of vertices and ∆ the largest vertex degree. For onehop graph coloring, the algorithm proposed in [START_REF] Johansson | Simple distributed (∆ + 1) -coloring of graphs[END_REF] runs in O(log n) rounds, but uses a number of colors close to ∆, whereas Distributed Largest First runs in O(∆ 2 log n) rounds [START_REF] Hansen | Distributed largest-first algorithm for graph coloring[END_REF].

Distributed deterministic two-hop coloring

In this paper, we are interested in two-hop coloring. Indeed, in wireless ad hoc and sensor networks, interferences are generally assumed to be limited to two-hops. Hence, two transmitters at a distance strictly higher than two transmission range can simultaneously transmit. The coloring algorithm used must be deterministic and distributed.

To extend a deterministic algorithm of one-hop graph coloring to two-hop graph coloring can rise some difficulties. Indeed a node can communicate directly only with its one-hop neighbors. The information coming from its two-hop neighbors is received two rounds later. Let N 2 (N ) denote the set of nodes at a distance up to two-hop from N , we can distinguish two classes of algorithms: the simplest ones, such as the extension of Largest First, are based on identical rounds, called decision round. In a round, a node sends a message to its one-hop neighbors, this message contains its color and the colors of its one-hop neighbors. It receives the messages from its one-hop neighbors and takes a decision if it has the highest priority. Its decision is based on decisions already taken by nodes in N 2 (N ) having a higher priority than N . The priority of a node is fixed and does not depend on the round. the more complex ones, such as Distributed Largest First and Dsatur, alternate proposal rounds and decision rounds. To propose a color, a node N must know all the decisions taken in the previous decision rounds, by nodes in N 2 (N ). To decide, a node N must know all the decisions taken in the previous decision rounds, and all the proposals made in the previous proposal round, by nodes in N 2 (N ).

The solution we propose belongs to the first class of algorithms that is simpler to implement and requires less messages, as illustrated by the comparative performance evaluation in Section 4.2. According to the state of the art, our algorithm will be based on the highest cardinality of the set of nodes at a distance up to two-hop. In case of equal cardinality, the node identifier will be used as a secondary criterion.

Slot assignment algorithms

Slot assignment in wireless ad hoc and sensor networks consists in assigning slots to wireless nodes in such a way that the same slot is never used by two transmitters at a distance less than or equal to two hops. In such a context, slot assignment and two-hop coloring are very close. Probabilistic algorithms exist, such as [START_REF] Bao | A new approach to channel access scheduling for ad hoc networks[END_REF] and [START_REF] Busch | Contention-free MAC protocols for wireless sensor networks[END_REF], where a node randomly chooses a slot. It can use this slot, only if no conflict is detected on this slot up to two-hop. In [START_REF] Herman | A distributed TDMA slot assiognment algorithm for wireless sensor networks[END_REF], the number of slots allocated to a node depends on the number of colors seen by this node. All these solutions are inadequate in case of non-uniform traffic distribution, because the number of slots allocated to a node does not depend on its traffic rate.

The first deterministic solution based on slot assignment is TRAMA [START_REF] Rajendran | Energy-efficient, collision-free medium access control for wireless sensor networks[END_REF]. It consists in a neighborhood discovery protocol, a schedule exchange protocol and an adaptive election algorithm that selects the transmitter and receiver(s) for each time slot. Only nodes having data to send contend for a slot; notice however, that a node does not know which of its 1-hop and 2-hop neighbors have data to send. The node with the highest priority in its two-hop neighborhood wins the right to transmit in the slot considered. Each node declares in advance its next schedule containing the list of winning slots and for each winning slot its receiver(s). This new schedule is declared in the last slot of the current schedule. Hence, each neighbor node has to listen to this slot. In order to tolerate packet loss, a schedule summary is sent in each data packet. The adaptivity of TRAMA comes at a price: its complexity.

FLAMA [START_REF] Rajendran | Energy-efficient, application-aware medium access for sensor networks[END_REF] optimizes TRAMA for data gathering applications in sensor networks. In such an application, each node (except the sink) has one outgoing flow in the data gathering tree rooted at the sink and one incoming flow per child. The protocol is simplified both in terms of 1) message exchange (the schedule is no longer sent, the receiver is implicitly the parent of the transmitting node) and 2) processing complexity (the priority of a node is the weighted sum of its traffic rate and a pseudo-random function of the node identifier and slot number). The number of slots allocated by FLAMA to a node with a given traffic rate highly depends on node priority computation.

With regard to slot assignment, we focus on a deterministic algorithm, less complex than TRAMA and more generic than FLAMA, ensuring that: each node is guaranteed to receive at least one slot, the number of slots granted to a node depends on its traffic rate.

Moreover, late node arrivals can create conflicts. The impact of a late arrival should be limited to the neighborhood of the joining node. SERENA allows any router node to stay awake only during its slots and those assigned to its one-hop neighbors and to sleep all the remaining time. The number of slots assigned to a router node is proportional to its traffic rate. In any case, any router node is ensured to have at least one slot per frame period. SERENA is decentralized and localized. It adapts to traffic and topology changes. In SERENA, any node N has a priority equal to the cardinality of the set of nodes up to two-hop, denoted N 2 (N ). In case of equal cardinalities, the node with the smallest identifier wins. SERENA consists of a two-hop coloring algorithm and a slot assignment algorithm.

The aim of the two-hop graph coloring is to color all nodes with the smallest number of colors, in such a way that two different nodes at a distance less than or equal to two-hop have not the same color. A node N can select its color if and only if all nodes in N 2 (N ) with a higher priority have already selected their color. It then selects the smallest color unused in N 2 (N ).

As soon as all nodes in N 2 (N ) are colored, the slot assignment algorithm starts. Its purpose is to assign time slots to each node. Any node N computes its number of slots depending on its traffic. More precisely, node N receives the slot reserved to its color and it computes its number k of additional slots according to the following formula:

k = traf f ic(N ) i∈V isibleColor(N ) traf f ic(i) * (Size -|V isibleColor(N )|) .
Among these k slots, kk are requisitionable, where

k = traf f ic(N ) M ∈N 2 (N ) traf f ic(M ) * (Size -|V isibleColor(N )|) .
where:

Size is the size of the frame;

|V isibleColor(N )| denotes the cardinal of the set of colors visible by N up to two-hop;

traf f ic(N ) is the bandwidth request of node N ; it is computed from the traffic submitted by N on the last period and the traffic pending on N ; traf f ic(i) denotes the highest bandwidth request of nodes having color i up to two-hop from N . Notice that several nodes in N 2 (N ) can have the same color: this is perfectly acceptable insofar as these nodes are not at a distance less than or equal to two hops. That is why the maximum of the bandwidth request must be taken in the computation.

A node N selects its k additional slots among the available ones if and only if all nodes in N 2 (N ) with a higher priority have already selected their additional slots. If all nodes succeed in assigning their additional slots, the algorithm is over. Otherwise, the node N that is unable to get at least its k additional slots requisitions them among the requisitionable slots used in N 2 (N ).

SERENA slot assignment algorithm is run periodically to adapt to traffic changes. The SERENA message sent by a node N to its one-hop neighbors contains the following fields: the node identifier, the node priority, the node color, the node traffic indication, the list of allocated slots and the indication whether they are requisitionable or not, and for each of its one-hop neighbors the node identifier, -the node priority, -the node color, -the node traffic indication, -the list of allocated slots and the indication whether they are requisitionable or not.

Depending on the algorithm progress, some fields can be missing.

We can notice three important features of this algorithm:

no node starvation: each node is ensured to get at least one slot per frame. In other words, the minimum throughput guaranteed to a node is equal to Bandwidth/Size, where Bandwidth denotes the network bandwidth and Size the size of the MAC slotted frame; node fairness: if traffic is uniformly distributed and all nodes have the same degree, they all receive the same amount of time slots; adaptivity to varying traffic rates and non-uniform traffic patterns, where a few number of nodes submit a high part of network traffic. Some optimizations of SERENA are possible such as ordering the messages to transmit in a slot, in such a way that a node can detect the soonest possible that no message is destinated to it. For instance, we can: put the broadcast messages at the beginning of a slot, order the point-to-point messages by increasing destination identifier.

Performance evaluation of SERENA two-hop coloring

In this section, we compare the performance of SERENA with two other coloring algorithms: Distributed Largest First, DLF, and an algorithm close to SERENA, called the smallest identifier algorithm. In this algorithm, the priority of a node is equal to its identifier and the uncolored node N with the smallest identifier in N 2 (N ) is colored next. For each algorithm, we evaluate the number of colors used and the time required to color all nodes in the network.

Simulations have been performed for different wireless networks, where the network density (the average number of neighbors per node) is fixed to 10. Nodes whose number varies from 50 to 200, are randomly distributed in the network area. Each result is the average of 5 simulation runs. Simulation results show that SERENA (the algorithm using the maximum degree) outperforms the algorithm with the smallest identifier, both in terms of number of colors used (see Figure 2) and time complexity expressed in number of rounds (see Figure 3). Intuitively, with SERENA (the maximum degree algorithm), the information (i.e. the selected colors) is propagated more quickly in the network: more nodes know the selected colors and can decide. If we focus on the first node selecting its color, let N be this node. Two rounds later, the maximum number of nodes up to 2-hop from N knows that this color is chosen. At least one of them, with the maximum degree, selects its color, and so on. With the identifier variant, the node N with the smallest identifier has a degree generally less than the degree of N .

These simulation results validate the design choice of SERENA to define the priority of a node N as the cardinality of N 2 (N ), denoted |N 2 (N )|.

Compared with Distributed Largest First, DLF, SERENA shows very good performance. Both algorithms use a similar number of colors, whereas the time complexity of our algorithm is significantly lower. For 200 nodes with a density of 10, the number of rounds in DLF is 228, more than twice the number of rounds with SERENA, 105.

We now study the impact of network density, (i.e. the average number of neighbors per node) on the numbers of colors and rounds used by DLF and our algorithm. We consider a network of 100 nodes, with a node density varying from 5 to 20. Simulation results are averaged over 5 simulations and illustrated in Figures 4 and5 for the number of colors and the number of rounds, respectively. rounds: the difference between SERENA and the smallest identifier algorithm tends to vanish, when the density increases. For a density of 20, both algorithms have the same number of rounds. For all the densities studied, DLF exhibits the highest number of rounds, whereas SERNA provides the smallest one. Hence, it provides a shorter convergence time, a very interesting property in wireless ad hoc and sensor networks where energy matters.

INRIA

Let c denote the number of colors assigned by SERENA, we have:

H1 + 1 ≤ c ≤ H2 + 1,
with H1 the maximum number of one-hop neighbors, H2 the maximum number of neighbors up to two-hop.

With SERENA coloring algorithm, a fully connected graph is colored with its chromatic number, that is c = H1+1 colors. A bipartite graph is also colored with its chromatic number, but with c = H2 + 1 colors. In a one dimension network, the maximum number of rounds of SERENA coloring algorithm is reached when the node identifiers are attributed as illustrated in Figure 6. It is equal to 2 * diameter -5, where diameter denotes the network diameter expressed in number of hops. In the best case, the node with the smallest identifier is located in the center, the number of rounds is then equal to diameter/2 + 2. Compared with Figure 2, Figure 4 shows that the number of colors used by our coloring algorithm depends stongly on the network density and weakly on the number of nodes. This observation is also true to a lesser extent for the number of rounds.

Performance evaluation of SERENA slot assignment

We first justify, why in SERENA some slots are requisitionable during the slot assignment algorithm. If the slot assignment algorithm were centralized, slots would have been assigned per color. Hence, the number k of slots would have been guaranteed to node N . Let us consider the following situation, illustrated in Figure 7.a, where color 1, used at node N 1 , is reused three hops away, at node N 4 . Let us assume a frame size of six slots. Node N 2 has color 2, whereas node N 3 has color 3. Three slots are allowed to node N 1 , two slots to nodes N 2 and N 4 , and one slot to node N 4 . Colors receive their slots in increasing color order. Moreover, the number of slots allocated to color i is equal to the highest number of slots granted to a node with color i. We would have:

1. color 1 would receive slots 1, 2 and 3, node N 1 too, whereas node N 4 would only receive slots 1 and 2;

2. color 2 would receive slots 4 and 5, node N 2 too; INRIA 3. color 3 would receive slot 6 as well as node N 3 .

. However, to implement such an algorithm in a wireless network would be expensive in terms of message exchanges and computing power required by each node. That is why, we prefer a distributed and localized slot assignment algorithm, where each node locally assigns its slots. In order to limit the size of the messages exchanged, the priority used in slot assignment is this used in coloring. If the localized slot assignment assigns k slots to node N , we can face the following situation, illustrated in Figure 7.b:

..---N1-------N2-------N3-------N4---... col1 col2 col3 col1 3 slots
1. node N 1 selects slots 1, 2 and 3; meanwhile, node N 4 selects slots 5 and 6, because slots 1, 2 and 3 have already been attributed to neighbors with a higher priority than N 4 ;

2. node N 3 selects slot 4;

3. it is then impossible to node N 2 to select its two slots: there is no more available slot.

More generally, let us consider the following case where two nodes with a high priority (nodes N 1 and N 4 in the example), at a distance up to 4 hops reuse the same color and the union of the slots selected by these nodes is not equal to the slots selected by one of them (the set of slots {1, 2, 3, 5, 6} is neither equal to {1, 2, 3} nor to {5, 6}. In such a case, nodes at a distance up to two hops from these two nodes (node N 2 in the example) will be unable to find k available slots.

We now illustrate the behavior of our algorithm by a short example on a small network of 10 nodes. This network is illustrated in Figure 8, where the number besides the node identifier denotes the color assigned to this node by SERENA. The frame size is set to 2 * |N 2 (N )| slots. The traffic for each node being given, the slot assignment algorithm computes the values of k and k , and assigns to each node a number of slots given in the last line of

k 0 1 1 0 2 0 1 3 2 0 k' 0 2 2 0 3 0 2 3 2 0 Slot 1 3 3 1 4 1 3 4 3 1
Table 2: Number of slots assigned per node According to our algorithm, each node is awake during its slots and the slots attributed to its neighbors. The percentage of activity time is provided for each network node in Figure 9. The most active node is the node E that generates the highest traffic rate. Furthermore, its one-hop neighbors, C, G and H largely contribute to the network traffic, explaining why E has the highest activity. Meanwhile, nodes A and F are active at only 30%. We focus on wireless and ad hoc sensor networks where: a high majority of nodes are present initially, topology changes are limited, mobility is also limited: only a few nodes can move and their speed is slow, the maximum number of nodes up to two hops is known.

It follows that the coloring algorithm is run once, at network initialization. However, late node arrivals or node mobility can cause color conflicts: two nodes at a distance less than or equal to two hops have the same color. Late node arrivals and node mobility are detected by topology changes. More precisely, some nodes detect changes in N 2 (N ). This set is refreshed by periodic exchanges of messages with the one-hop neighborhood (see rule RS3). A new joining node has no color assigned, whereas a moving node generally has one. Two types of conflict can be created: first conflict type: the moving node becomes the one-hop neighbor or twohop neighbor of a node with the same color. Clearly, this conflict type is possible only in case of mobility.

second conflict type: two nodes that are at a distance between one and two transmission ranges become two-hop neighbors due to the late arrival of a new node or node mobility.

General case

In the general case, the algorithm proceeds in three steps to deal with late node arrival or node mobility: step 1: Conflicts are first detected. The detection is based on the exchange of the colors granted to the node and its one-hop neighbors. If a node N detects that one of the node in N 2 (N ) has the same color as itself, a conflict exists. step 2: Detected conflicts are solved: the node with the highest priority in the conflict keeps its color, the other conflicting nodes apply rule RC1 to get a new color. In this step, the new node in case of late arrival will also receive its color. step 3: This will cause each node seeing a change in colors and bandwidth requests to reassign the slots according to SERENA slot assignment algorithm.

It follows that the impact of a new joining node N is limited to N 2 (N ). This general algorithm can be simplified in some specific cases. Indeed, if there is no node mobility and any node knows all the nodes at a distance higher than one transmission range and less than two, then no conflict is created.

Specific case

In this specific case, we assume that there is no node mobility and any node knows all the nodes at a distance higher than one transmission range and less than two. This latter assumption can be met because: any node knows its GPS coordinates as well as the coordinates of any node at a distance up to two transmission ranges. This can be achieved by the routing protocol; or the network configuration is such that for any node N , there is no node M at a distance higher than one transmission range and less than two having no common one-hop neighbor with N .

In this case, the late arrival of a node N does not cause a color conflict and two steps are sufficient: step 1: A color is granted to this new node, all the other nodes keep their color. According to SERENA coloring algorithm, the new node N selects the smallest color unused in N 2 (N ).

step 2: This will cause each node seeing a change in colors and bandwidth requests to reassign the slots according to SERENA slot assignment algorithm.

The impact of the changes are limited to the nodes in N 2 (N ).

Network lifetime and user data delivered

We now evaluate the performance of SERENA in various wireless ad hoc and sensor networks. More precisely, we quantify the network lifetime obtained with and without SERENA. In both cases, the routing strategy RE is used. The network lifetime is defined as the time at which a flow destination becomes unreachable. The simulation parameters are those given in Section 2. Furthermore, the slot size is fixed to 12ms and the frame consists of 80 slots. The frame size is equal to 2 * 4 * density, where 4 * density is the average number of nodes up to two hops. Simulation results are averaged on 5 simulation runs.

It appears that with SERENA, the network lifetime can be increased by 100% (see for instance a network of 100 nodes in Figure 10). Notice that this excellent improvement would be useless, if during this extended network lifetime, the amount of information delivered to the user was not increased. Figure 11 confirms that with SERENA, the increase in network lifetime is followed by an increase in the delivered data.

Slot reuse and impact of frame size

We now focus on the slot assignment algorithm and evaluate its performance by means of simulation. We consider a network of 50 nodes. The size of the node waiting queue is set to 50. The other simulation parameters are unchanged.

Traffic is not uniformly distributed over the nodes. We evaluate the number of slots obtained by a node with regard to its traffic rate. Figure 12 depicts the number of slots assigned to ten nodes. We notice that the nodes with the highest traffic rates (see nodes 3 and 6) receive the highest number of slots (7 and 11 slots respectively), whereas nodes with a small traffic rate (see nodes 1 and 9) receive few slots (3 slots for both of them). Simulation results show that as expected, SERENA assigns a slot number proportional to the traffic rate of the node.

Figure 13 depicts the slot reuse. It provides the number of slots shared by 5 nodes, down to 0 node. Two slots among the 80 are empty and 33 slots are reused by three nodes. Three quarter of slots are reused by several nodes. A slot is in average used by 2.43 nodes. This spatial reuse is obtained thanks to the SERENA two-hop coloring and slot assignment algorithms. It follows a better network efficiency.

INRIA

We now study the impact of frame size on both network lifetime and amount of messages delivered. The number of slots in the periodic frame takes the values of 40, 80, 160 and 240. As depicted on Figure 14, the highest lifetime is obtained with a frame size of 40 slots. Network lifetime then decreases and finally stabilizes for 160 slots. However, the size of 40 slots is unsatisfactory because of an unacceptable message loss rate, as shown in Figure 15. Indeed, nodes with a high traffic rate have not enough slots to send their messages and their waiting queue overflows. That is why, we recommend the use of a frame of 80 slots, corresponding to the highest lifetime obtained with no message loss due to queue overflow. From now on, the frame size is set to 80 slots. The impact of slot size is analyzed. We evaluate network lifetime and delivery rate, when the slot size ranges from 6ms to 48ms. We observe in Figure 16 that the slot size has no impact on network lifetime, as soon as it is higher than or equal to 12 ms. A slot of 6ms provides the best network lifetime. However, as in the previous case, it leads to an unacceptable loss rate due to queue overflow (see Figure 17). Indeed, this slot size allows the transmission of only 3 messages. We also see that the size of 12 and 24ms give similar results with regard to the amount of user information delivered. With a size of 48ms, some messages are lost. Hence, the best value for the slot size is 12ms.

Buffer dimensionning

We now evaluate the average and maximum number of messages in the waiting queue of a wireless node. The flow throughput ranges from 8kbps to 32kbps. Figure 18 depicts the results obtained, when the size of the node waiting queue INRIA is set to 50. As long as the flow throughput is lower than or equal to 16kbps, no message is lost due to queue overflow and the maximum number of waiting messages per node is less than 10 for a throughput of 8kbps and 20 for 16kbps. The average value is less than the half. With higher throughputs, the network is saturated.

Distribution of node energy consumption

We now perform the same simulations as those described in Section 2, but now, using SERENA. The new distribution of node energy consumption is illustrated in Figure 19. The bar diagrams represent the energy dissipated in the Transmit, Receive, Idle, Overhearing and Interference states succesively for different size of networks: 50, 100, 150 and 200 nodes. The network density is set to 10.

Figure 19 shows that with SERENA, the energy dissipated in the idle state decreases to about 30% instead of 50%. This is the first benefit brought by SERENA. The second energy cost comes from overhearing: about 20% with SERENA instead of 10% without. We expect to improve the energy dissipated in the idle and overhearing states by optimizing SERENA as explained at the end of Section 4.1. The third enegy cost is due to interferences: about 1%. It was the second one without SERENA, with about 20%. This is explained by the fact that with SERENA, if neither the node nor its one-hop neighbors are transmitting, the node is sleeping. Hence, SERENA contributes to significantly reduce the interference phenomenon. This is the second benefit of SERENA.

Conclusion

In this paper, we have highlighted that in wireless and ad hoc sensor networks, most of the node energy is dissipated in the idle state. The second energy cost comes from interferences, even if a routing efficient strategy is used. SERENA in allowing a node to sleep while neither it nor one of its one-hop neighbor is transmitting, contributes to considerably improve the energy efficiency. First, the energy dissipated in the idle state is reduced up to 40%. Second, the energy loss due to interferences becomes negligible (about 1%). In parallel, the node spends more of its time in the Transmit and Receive states: the only useful states from the application point-of-view.

E

  transmit = P transmit * Duration, E receive = P receive * Duration, Duration is the transmission duration of a packet.

Figure 1 :

 1 Figure 1: Distribution of node energy consumption without sleeping state.

INRIA 4

 4 Presentation of SERENA 4.1 Principles of SERENA

INRIA

Figure 2 :

 2 Figure 2: Number of colors used

Figure 3 :

 3 Figure 3: Number of rounds

Figure 4 :

 4 Figure 4: Number of colors used

Figure 5 :

 5 Figure 5: Number of rounds With regard to these figures, SERENA provides an excellent performance, both in terms of: colors: it uses the smallest number of colors for densities 5 and 10 and is very close to the smallest number for density 20.

Figure 6 :

 6 Figure 6: Number of rounds in the worst case

Figure 7 :

 7 Figure 7: Distributed vs centralized slot assignment algorithm

Figure 8 :

 8 Figure 8: A wireless network and its coloring

Figure 9 :

 9 Figure 9: Node activity time

4. 4

 4 Late node arrivals and node mobility 4.4.1 Assumptions

INRIA

Figure 10 :

 10 Figure 10: Network lifetime with and without SERENA.

Figure 11 :

 11 Figure 11: Amount of user data delivered with and without SERENA.

Figure 12 :

 12 Figure 12: Slot assignment and traffic rate.

Figure 13 :

 13 Figure 13: Slot reuse.

Figure 14 :

 14 Figure 14: Impact of frame size on the network lifetime.

Figure 15 :

 15 Figure 15: Impact of frame size on the delivered data.

Figure 16 :

 16 Figure 16: Impact of slot size on the network lifetime.

Figure 17 :

 17 Figure 17: Impact of slot size on the delivered data.

Figure 18 :

 18 Figure 18: Average and maximum number of waiting messages per node versus flow throughput.

Figure 19 :

 19 Figure 19: Distribution of energy consumption with SERENA.

  Energy efficient strategies . . . . . . . . . . . . . . . . . . . . . . 3.2 Node activity scheduling . . . . . . . . . . . . . . . . . . . . . . . 3.3 Centralized and distributed one-hop coloring . . . . . . . . . . . 3.4 Complexity of a distributed coloring algorithm . . . . . . . . . . 3.5 Distributed deterministic two-hop coloring . . . . . . . . . . . . . 3.6 Slot assignment algorithms . . . . . . . . . . . . . . . . . . . . .

	Contents 1 Introduction
	1 Introduction
	2 Analysis of the node energy consumption distribution
	3 State of the art
	3.1 5 Network lifetime and user data delivered
	6 Slot reuse and impact of frame size
	7 Buffer dimensionning
	8 Distribution of node energy consumption
	9 Conclusion
	RR n 6388

4 Presentation of SERENA 4.1 Principles of SERENA . . . . . . . . . . . . . . . . . . . . . . . 4.2 Performance evaluation of SERENA two-hop coloring . . . . . . 4.3 Performance evaluation of SERENA slot assignment . . . . . . . 4.4 Late node arrivals and node mobility . . . . . . . . . . . . . . . . 4.4.1 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 General case . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.3 Specific case . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 2 .

 2

	Node	A	B	C	D	E	F	G	H	I	J
	Degree	4	5	7	6	8	7	7	5	5	4
	Color	1	4	2	5	1	3	4	3	2	5
	Traffic	10 20 30 10 50 10 30 50 30 10

RR n 6388

We have shown how to tune various parameters to maximize both the network lifetime and the amount of user data delivered. In conclusion, SERENA achieves a more efficient usage of node energy.


Cyril Mokrani
email: cyril.mokrani@meric.cl

Mireille Bossy
email: mireille.bossy@inria.fr

Marcos Di Iorio
email: marcos.diiorio@meric.cl

Antoine Rousseau
email: antoine.rousseau@inria.fr

Numerical Modelling of Hydrokinetic Turbines Immersed in Complex Topography using Non-Rotative Actuator Discs

Keywords: Tidal and oceanic turbulent flow, Stochastic numerical model, Hydrodynamics, Turbulent flows, Stochastic Lagrangian models, Marine current energy converters

Technologies for extracting energy from tidal currents using hydrokinetic devices have not yet achieved a sufficient degree of maturity to enable their widespread integration into power grid interconnected systems [START_REF] Khan | Ocean wave and tidal current conversion technologies and their interaction with electrical networks[END_REF], [START_REF] Güney | Hydrokinetic energy conversion systems: A technology status review[END_REF]. Recent studies have identified potential in many coastal areas which could provide significant renewable energy resources in the near future, although their implementation would require the overcoming of challenging technical, environmental and economical aspects, all of them site-specific.

Along the 6,400 km coast of Chile, in the South Pacific Ocean, the Chacao Channel has, according to initial estimates [START_REF] Cruz | [END_REF], been identified as an important energy source that could potentially contribute to a total power capacity of between 600 and 800 MW. Motivated by the need to learn more about the physical and environmental processes of highly energetic areas similar to the Chacao Channel, we are developing a numerical model that aims to comprehensively describe flow past turbines for a wide range of scales, to optimise turbine arrays and evaluate their environmental effects.

To improve device installation and performance through numerical predictions, it is necessary to model turbulence flows. Indeed, the Chacao Channel presents current velocities of up to 10 m/s at high tide and a complex bathymetry with high depth variations, which naturally generates turbulent flows. In addition, this ambient turbulence may be significantly increased by turbine blade interaction with the incident flow, mainly due to high shear stress generated downstream from the turbine.

Among the existing numerical methods for modelling turbulent flows, one can separate the Eulerian Averaged Models (EAMs) and Lagrangian Stochastic Models (LSMs). Those two approaches are fundamentally distinct for various reasons: first, Eulerian methods compute solutions (velocity and pressure fields) on a fixed mesh grid (see Figure 1-(a)) and thereby provide a snapshot of the flow at every instant. Conversely, Lagrangian methods follow water particles during the simulation and compute their velocity at each time step. Contrary to Eulerian models, the Lagrangian approach provides the trajectory (positions) of fluid particles (see Figure 1-(b)).

In addition, EAMs provide pressure and velocity fields averaged over a characteristic duration or length (chosen a priori i.e., as input parameter ). This is equivalent to saying that only large-scale turbulent processes are solved. Small-scale effects are included in the EAM's description, however they are "added" to the numerical solutions through an additional input variable (see Reynolds tensor or sub-grid tensor). Conversely, the LSM computes instantaneous fields of each water particle, and provides solutions which include both large and small-scale turbulent effects. An important point is that it is possible to deduce averaged fields from LSM solutions but not vice-versa.

Among the EAMs, the Reynolds-averaged Navier-Stokes (RANS) [START_REF] Zahle | Characterization of the unsteady flow in the nacelle region of a modern wind turbine[END_REF] and Large-Eddy Simulation (LES) methods [START_REF] Sezer-Uzol | 3-d time-accurate cfd simulations of wind turbine rotor flow fields[END_REF] have been widely developed to model turbulent flows induced by tidal turbines. LES have been used to model blade motions, predict forces and torque induced on the stator (Sezer- [START_REF] Sezer-Uzol | 3-d time-accurate cfd simulations of wind turbine rotor flow fields[END_REF][START_REF] Kang | Numerical simulation of 3d flow past a real-life marine hydrokinetic turbine[END_REF], and improve the prediction of the 2 December, 2018 power generated by the turbine [START_REF] Kang | Numerical simulation of 3d flow past a real-life marine hydrokinetic turbine[END_REF]. RANS models have also provided excellent predictions, especially regarding the local flow/structure interactions occurring during the extraction process.

However, these methods have proved to become very intensive and demanding in terms of computational resources. Accurate predictions of RANS require at least 6 million grid nodes, with a large number concentrated in the wake region [START_REF] Turnock | Modelling tidal current turbine wakes using a coupled rans-bemt approach as a tool for analysing power capture of arrays of turbines[END_REF]. This requirement can be even more important for LES models, for which a total of 185 million grid nodes are required to resolve the flow field around one single turbine device [START_REF] Kang | Numerical simulation of 3d flow past a real-life marine hydrokinetic turbine[END_REF]. These computational needs clearly lead to prohibitive computation processing times, which becomes unfeasible if multiple interacting turbines are considered in the simulation. So far, the need to reduce the demand on computer resources for EAMs has motivated a number of researchers to propose innovative methods such as coupling approaches or hybrid methods with the objective of saving computational resources associated to EAM.

Contrary to EAMs, LSMs are more recent and fewer studies have been published so far. LSMs have been originally developed by [START_REF] Pope | Turbulent flows[END_REF], based on stochastic differential equations whose solutions are the instantaneous positions and velocities of fluid particles. The main usage of LSM models in engineering is to provide information for the instantaneous velocity of the flow seen by dispersed physical particles. Such approaches have been extensively used for turbulent flows in combustion engineering, in association with an EAM for the computation of the mean flow information [START_REF] Minier | Statistical descriptions of polydisperse turbulent two-phase flows[END_REF]. Recently a stand alone version of LSM has been developed as a simulation solution for turbulent flows and mainly applied to atmospheric flows: [START_REF] Polagye | Environmental effects of tidal energy development[END_REF] used a LSM for meteorological frameworks and predicted turbulent flows fields in a realistic case of wind refinement. Originally, standalone LSM for wind computation is a downscaling method, that provides a local model limited to a given region, forced by meteorological large scale information. In order to reconstruct the mean flow field from the fluid-particle motion, a Particle-in-cell (PIC) technique is retained (see further details in [START_REF] Bossy | Modeling the wind circulation around mills with a Lagrangian stochastic approach[END_REF]) that introduces a fixed mesh according to some characteristic flow length. [START_REF] Zahle | Characterization of the unsteady flow in the nacelle region of a modern wind turbine[END_REF]; [START_REF] Bossy | Stochastic Lagrangian approach for wind farm simulation[END_REF] computed the wind circulation around mills and obtained very promising estimations of the wind variability in the wake of a wind turbine. However, even now, no studies based on tidal or ocean flow are available and the question of LSM ability to describe turbulent flows generated by tidal turbines remains unresolved. More important is the assessment of their performance. Given the high computational times associated to EAMs, reducing computational times should be very welcome.

The present study aims to evaluate the performance of the LSM "Stochastic Downscaling Method -Ocean Power Simulation (SDM-OceaPoS)" with the purpose of developing an inexpensive computation tool capable of simulating interaction between high Reynolds flows with both multiple turbine devices and complex bathymetry.

3

December, 2018

Section 2 provides details of the numerical methods employed in this study. Section 3 presents one verification case based on the boundary layer generation. Sections 4 and 5 present two benchmarks, based respectively on turbulent flows generated by a bathymetry (from [START_REF] Almeida | Wake flows behind two-dimensional model hills[END_REF]; [START_REF] Loureiro | Water tank and numerical model studies of flow over steep smooth two-dimensional hills[END_REF]) and instantaneous wake dynamics generated by an non rotative turbine (from Myers andBahaj (2010, 2012)). Since Section II is a preliminary study, conversely Sections IV and V intend to establish a framework for future optimisation of turbine arrays on specific sites.

NUMERICAL METHODS

Numerical approach

The development of the SDM-OceaPoS model is based on an existing code that models air flows named SDM-Windpos. The "air" version has been validated on atmospheric flows [START_REF] Bossy | Modeling the wind circulation around mills with a Lagrangian stochastic approach[END_REF], although its translation to the water version is still under development.

In the current version, the water is described as an incompressible flow by a Lagrangian approach. The Lagragian description of the flow uses stochastic diffusion processes given as the solution of a stochastic differential system (we refer the interested reader to [START_REF] Bernardin | Stochastic Lagrangian Method for Downscaling Problems in Computational Fluid Dynamics[END_REF] for a mathematical description of the method, and to Pope (2000, Ch. 12) to a physical introduction to this methodology). The time evolution of particle velocities is described and computed from three acceleration terms, including (i) the pressure gradients contributions, (ii) the advective and convective effects of the fluid displacement, and (iii) a diffusive term governed by a Brownian motion. These three terms allow to take both small and large scale turbulent processes into account in the simulations.

One of the main advantages of the SDM-OceaPoS approach is that in addition to providing the instantaneous fields particles, it also provides averaged fields, similar to a classical EAM. This Eulerian description of the flow is possible thanks to a fixed regular mesh grid superposed to the numerical domain. During the computation, spatial averages of Lagrangian velocities of particles located in the same cell are assessed, together with the variance and covariances. These Eulerian fields are computed at each time step after particle advections. In practice the number of particle in each cell keeps constant in order to guarantee the incompressibility constraint. This step is performed using particle optimal transport algorithm (see [START_REF] Chauvin | Stochastic Downscaling Method: Application to Wind Refinement[END_REF] for a description of proposed algorithms and [START_REF] Bossy | Local Existence of Analytical Solutions to an Incompressible Lagrangian Stochastic Model in a Periodic Domain[END_REF] for a mathematical description of incompressible Lagrangian stochastic model). Using this approach, the user has to set two input parameters which are (i) the mesh size used to assess the Eulerian fields, and (ii) the number of particles per cell.

Method for including bathymetry effects on the flow

SDM-OceaPos uses a wall treatment method to reproduce the momentum exchange present in the interaction between fluid particles and arbitrary bathymetries. The method implemented consists of reflecting any particle that enters into a region close to the ground (namely the "mirror region"). Outside this region, particles are in the fluid domain and describe the fluid dynamics. Inside the mirror region, the flow not is solved, however its effects are included in the fluid domain simulations. To make it possible, a 3-dimensional topology file is first generated prior to the computation with the main ground characteristics (that is to say, the ground position, the ground local inclinations, and the frictions coefficients). During the computations, SDM-OceaPos assesses each reflected velocity of particle that cross the mirror region from the input ground data. This step is performed by assuming that the covariances before and after the reflection keep equal. This method is applied in three steps : (i) SDM-OceaPos first detects particles that enter the sub-viscous layer at a given iterative step. During this step, the position the particle would reach without terrain (the virtual position) is computed. (ii) The position of the reflected particle is calculated according to the local ground inclination, and hit virtual position. The particle is then re-introduced into the fluid domain. (iii) The reflected velocity vector is computed according to the friction and the covariances assessed before the 4 December, 2018 reflection.

This original method is applied in Section III on a case involving a 2D hill introduced in a turbulent flow.

Model of hydrokinetic turbines

The implemented model of a non-rotative turbine is based on the theory of porous discs (Myers andBahaj, 2010, 2012). Such discs, introduced into a turbulent flow, may disturb the flow field by generating a turbulent wake as illustrated in Figure 2-(left).

Turbine models were initially implemented in the air version of SDM [START_REF] Bossy | Modeling the wind circulation around mills with a Lagrangian stochastic approach[END_REF]. They consisted in applying a correcting velocity to particles located at the disc location. This correcting velocity was assessed from an theoretical expression of the force generated by the disc on the fluid. However several problems related to this approach have been pointed out: first, this corrected velocity was added after the advection step, and therefore the turbulence generated upstream from the disc was generated indirectly by the velocity deficit, but not by the disc himself. This may lead to a misrepresentation of the wake recovery and impacts in realistic scenarios. Other limitations were that the corrected velocity was applied in a fixed control volume -typically the disc volume -and therefore any particle crossing the disc without staying in this control volume was not corrected. This not only provided a non realistic behaviour of the fluid particle, but also resulted in a flow description of flow very sensitive to the space step and the time step.

A new approach has been developed to solve these problems. It involves adding an acceleration source term to the momentum equation. In this manner, the turbulent effect induced by the turbine itself can be included in the simulation. Additionally, this source term is applied to a new control volume, defined as the volume of all particles that cross the disc between two iterative steps. With this variable control volume any particle that crosses the disc "feels" its effect.

From a mathematical point of view, the source term acceleration is assessed from the ratio between the total force applied on the control volume and its corresponding total mass, as illustrated in Figure 2-(right).

Figure 2: Illustration of the wake generated by the porous disc (left), and (right) the approach of the implemented model: particle accelerations are estimated from the total force applied on the volume V T . The direction of the particle's force is aligned with the disc-crossing particle velocity.

5

December, 2018

• The expression of the total force is given in [START_REF] Simisiroglou | The actuator disc concept in phoenics[END_REF]. For the present study, this force is assumed to be exerted in the direction of the disc-crossing particle velocity, as indicated in Figure 2. Usually the force is given according to the porosity coefficient, the area of the disc and the fluid velocity vector in the disc's vicinity. However, for the present model we used the induction factor (denoted a, see [START_REF] Bossy | Modeling the wind circulation around mills with a Lagrangian stochastic approach[END_REF]) instead of the porosity coefficient because this coefficient allows for the linking of the instantaneous particle velocities at the disc location (see U D in Figure 2-(left)) to the particle velocity away from the disc (see

U ∞ in Figure 2-(left)).
• The total mass of the control volume is estimated by the product between the total number of particles that cross the disc and the mass of one single particle. The latter is calculated by assuming that the volume of each grid cell equals the sum of the volumes of all particles in this cell. From this assumption, the particle mass can be easily deduced.

Finally, the corrected acceleration term is re-assessed at each time step for each particle that crosses the disc. This approach allows for obtaining a flow description independent from the numerical input parameters: when changing the mesh grid size (respectively the time step), the control volume changes and the subsequent acceleration keeps the same order of magnitude so that the effect of the disc on the flow remains similar.

In the following, one verification test (Section 3) and two benchmarks (Sections 4 and 5) are presented. All numerical computations have been performed with the multi-threaded open Multi-processing (openMP) version of SDM-OceaPoS, and using 23 cores, having a frequency of between 2 and 2.4 GHz. The hydrokinetic turbine model is applied in Section 5 for cases involving turbine diameter D of about one third of the depth h (D h/3).

The next section is dedicated to benchmarks carried out on cases of flows without turbines. Emphasis is placed on the generation of boundary layers and the effects of a given bathymetry on a turbulent flow.

GENERATION OF THE BOUNDARY LAYERS

A first study has involved considering verification test cases in order to assess the ability of SDM to re-produce the basic physics of a turbulent flow without turbine. The objective is to re-create experiment in channel and assess the numerical accuracy of the computed velocity fields by comparing numerical and theoretical results.

Experiments

The following results are based on the experimental study from [START_REF] Nezu | Open-channel flow measurements with a laser doppler anemometer[END_REF] Results from [START_REF] Nezu | Open-channel flow measurements with a laser doppler anemometer[END_REF] show the generation of mainly four boundary layers in a steady state (see Figure 3-b) whose structures depend both on the ground roughness and interface motion. The ground roughness generates the viscous sub-layer in which the velocity profile follows a linear law. Above this first boundary layer, a transition region appears in which the velocity profile tends towards a classical "log" law. Above the log-law region, the free-surface motion disrupts the velocity profile resulting in a new boundary layer (see [START_REF] Nezu | Open-channel flow measurements with a laser doppler anemometer[END_REF]). For the present study, one focuses only on the intermediate, and log-law regions and neglects the effect of the free-surface and the sub-viscous layer. The latter assumption is justified 6 December, 2018 by the fact that the sub-viscous layer has very small dimensions in relation to those of the turbine deployed in this study (D h/3).

Results

Numerical tests have been made using the parametrisation of the turbulent viscosity proposed in the study from [START_REF] Nezu | Open-channel flow measurements with a laser doppler anemometer[END_REF] derived from laboratory measurements. Computations have been performed using 128, 000 cells in total, and 80 particles per cell in a numerical domain of 1 m 3 . Figure 4 presents the mean velocity profiles and variances of the vertical and horizontal velocities. Theoretical solutions are also represented in the same figure in order to assess the numerical accuracy to model the boundary layers. As observed, results from SDM-OceaPos provide excellent trends of velocities in the intermediate and outer log-law regions of the flow. Mean velocity predictions are slightly under-predicted in the interface vicinity, since the effect of interface on the flow is neglected. Also, the sub-viscous layer is poorly modelled due to the assumption previously mentioned. Maximal relative error according to theoretical results reached less than 10%.

Variances of horizontal and vertical velocity are compared to the mean velocities fluctuations in Figure 4 (b) and (c). Both global trends and order of magnitudes are in good agreement with theoretical values. Maximal error reach about 13% compared to theoretical results.

Regarding the numerical performances, all computations have been performed over 200 iterations using a constant time step of 0.2 s. The corresponding time calculation for one iteration is about 0.05 s, which means that one second in the physical experiment could be modelled in about 0.25 s.

Finally, these predictions, together with the corresponding performances observed lead to the conclusion that the main structure of a high Reynolds number flow can be accurately reproduced at very reasonable computation times.

INFLUENCE OF BATHYMETRY ON THE FLOWS

The second tests are focused on the influence of the bathymetry on the flow, modelled by using the numerical method presented in sub-section 2.2. It involves reproducing the experiments performed in the laboratory and comparing them to the numerical predictions. The objective is to assess the ability of SDM-OceaPos to describe the turbulent flow dynamics involved in controlled configurations.

7

December, 2018

Experiments

All numerical results presented hereafter have been compared to experimental measurements performed by [START_REF] Almeida | Wake flows behind two-dimensional model hills[END_REF]; [START_REF] Loureiro | Water tank and numerical model studies of flow over steep smooth two-dimensional hills[END_REF]. These studies consider a 2D hill placed in a water tunnel to measure its effects on turbulent flow (see Figure 5-a). [START_REF] Loureiro | Water tank and numerical model studies of flow over steep smooth two-dimensional hills[END_REF].

The water tunnel is made of perspex and is 0.170 m deep, 0.2 m wide, and 7 m long. During the experiments, the water flow passes through straighteners located upstream of the duct inlet in order to ensure hit uniformity. Experiments correspond to Reynolds number flows of about 6.10 4 . The model hills are made of plastic moulded with a wooden former. Their shape corresponds to the inverse of a fourth-order polynomial function, with a maximal height and length in the streamwise direction of 28 mm and 108 mm respectively.

During the experiment, the curvature of the flow affects the velocity field, which results in the establishment of a zone of high-intensity velocity fluctuations along the shear layer downstream of the top of each hill. The hydrodynamical effects of the hill on local turbulence can be observed in Figure 5, especially the emergent structures upstream and the recirculation zone downstream of the hill. Two specific points, 8

December, 2018 named detachment and reattachment points, are defined as the positions where the shear stress changes sign. These points determine the main flow structure generated by the hill. According to [START_REF] Loureiro | Water tank and numerical model studies of flow over steep smooth two-dimensional hills[END_REF], the detachment point extends up to x/h = 4.8, as can be observed in Figure 5 (b).

Results

The following results correspond to computations performed using 96, 000 cells in total and 50 particles per cell. To reduce convergence times, a first computation was done without a hill and stopped when all the mean velocity, first and second order moments fields reached stabilised trends. These fields were imposed both as a boundary at the left side of the domain, and initial state in the whole domain.

Simulations including the hills have been performed using smaller time steps and more iterations to obtain converged results. In total, we ran 2, 000 iterations which lasted about 5 hours. This corresponds to a CPU rate of 9 iterations per second.

Figure 6 shows the stream lines obtained after 1, 000 iteration steps. It can be observed that the recirculation zone is well described by the code and that its dimensions are very close to those observed experimentally. The detachment point location has been assessed at 4.28 h against 4.8 h from [START_REF] Loureiro | Water tank and numerical model studies of flow over steep smooth two-dimensional hills[END_REF]. The generation and stabilisation of this recirculation zone proved to be sensitive to the time step parameter, especially when low values of this parameter were set. Using the value dt = 0.0005 s provided a stable computation and resulted in positive alignment with laboratory data. 9

December, 2018

Figure 7 presents the vertical profiles of the mean horizontal velocity predicted at different locations upstream and downstream of the hill. As observed, numerical predictions are very close to measurements for all the considered distances to the hill. Maximal error on the mean velocity is about 10% and is reached at the hill top.

The interested reader can find more information about this benchmark case in Di Iorio et al..

CASES INVOLVING NON-ROTATIVE TURBINE(S)

The following validation tests are based on experiments from Myers andBahaj (2010, 2012) involving porous discs introduced in a turbulent 2D flow. In order to estimate the ability of SDM-OceaPos to describe both the velocity deficit and the turbulent intensity generated by one or various porous disc(s), we reproduced similar experiments using the method presented in Sub-section 2.2.

Experiments

Myers and Bahaj (2010, 2012) carried out flow measurements in the wake of the discs to determine the flow recovery as a function of the thrust produced by the devices and analyse flow interaction in turbine arrays. The experiments were performed in a 21 m long by 1.35 m wide flume, with a water depth of H = 0.3 m. The mean velocity of the flow equals 0.25 m/s, which yields a Reynolds number of Re = 75, 000 and a Froude number of Fr = 0.14 based on the flow depth H. Measurements were obtained by varying the number of discs, the disc immersion depth and the disc porosity. Figure 8 shows the experimental set-up with one and three discs in the laboratory channel. Diameter of the actuator discs (denoted AD in the following) remained constant, equal to D = 0.1 m for all tests that were performed. Porous discs were located in a section of the flume which had already reached a fully-developed turbulent channel flow. They were mounted on a rig with a pivot to amplify the forces and measure the total thrust with good accuracy. To vary the disc thrust, different porosities from 0.36 to 0.91 were tested using the same upstream flow conditions. The data provided also includes tests performed with different disc immersion depths. In any case, the measurement area extends from 3D to 20D downstream where instantaneous velocities were measured using a vertical resolution of 0.1H.

10

December, 2018

Results

The following results correspond to computations performed on a 275, 000 points grid and 100 particles per cell. For all cases, simulations were performed over 500 iterations with a time step of 0.1 s, which corresponds to a CPU rate of about one iteration per 10 s. The following figures correspond to results obtained at the last iteration of the computation.

Validation tests have been conducted by comparison analysis based on turbulence intensity and velocity deficits. The velocity deficit quantifies the velocity reduction induced by the porous disc. The results presented are also compared to numerical predictions from [START_REF] Johnson | Computational actuator disc models for wind and tidal applications[END_REF]. The latter used a RANS numerical code with a k -ω turbulent model and a source term in the momentum equation to include the presence of the porous disc.

Figure 9 presents the vertical profiles of turbulent intensity computed and measured at different positions upstream from the turbine. For this case, the porous disc immersion depth was half the total domain height. In general, numerical tests performed by SDM-OceaPoS provide very good predictions for all available porosity values. The turbulent intensity is slightly under-estimated for low values of porosity coefficient. Conversely, they are overestimated for high values of C T , especially in regions close to the disc position. Relative errors decrease away from the disc, reaching about 10% for C T = 0.62 and 8% for C T = 0.99. Numerical predictions from [START_REF] Johnson | Computational actuator disc models for wind and tidal applications[END_REF] provide a maximum relative error of 40% for turbulence intensity and 30% for speed deficit. In addition, they do not correctly reproduce turbulence intensity in the disc vicinity. Results obtained with SDM-OceaPoS are significantly better at this location.

Other results obtained on velocity deficit (not presented here) revealed corresponding maximum relative errors of the order of 7% for any porosity coefficient.

Figure 10 shows predictions of speed deficit where C T = 0.86 and the disc is immersed at mid-depth. Comparisons with experimental data reveal that SDM-OceaPoS provides a good prediction of the spatial distribution of the speed deficit. The wake shape and dimensions are close to the experimental data. Low errors are observed, especially far from the disc or in the disc vicinity, where the deficit speed is slightly overestimated. Overall, the numerical results are fairly satisfactory.

Similar conclusion can be drawn for cases involving various discs. Figure 11 presents results obtained for one case involving three discs disposed in a two-row array. The first row is two discs at 1.5D lateral separation and the second one is the third disc disposed downstream at a distance of 3D. As observed in Figure 11, for a downstream distance higher than 7D, numerical prediction of the velocity is in good agreement with experimental measurements. However, for a shorter distance, the numerical prediction is significantly over-estimated. For this case, further studies are needed to understand and improve results.

CONCLUSION

The present research aims to reproduce basic turbulent flows generated in rivers and seas both upstream and downstream of an AD device, by using a Lagrangian stochastic model. From the results obtained from one verification test and two benchmarks, it can be concluded that the SDM-OceaPoS model is able to accurately represent the flow structure generated at high Reynolds numbers, the effects of bathymetry on the mean velocity field in 3D cases and the effects of AD on the local turbulence generated past the turbine device. More specifically:

• The flow structure generated in a domain volume of about ∼ 1 m 3 at a Reynolds flow of Re 10 4 could be described with a CPU rate of about one iteration per 0.05 seconds, with all computed mean horizontal and vertical velocities and variances properly predicted in intermediate and external boundary layers. These results were obtained by implementing the parametrisation of turbulence viscosity 11 December, 2018 of [START_REF] Nezu | Open-channel flow measurements with a laser doppler anemometer[END_REF], which does not include the presence of the interface, even if this does not seem essential for current industrial applications.

• A new method based on the reflection of particles has been validated using the experimental cases of [START_REF] Almeida | Wake flows behind two-dimensional model hills[END_REF]; [START_REF] Loureiro | Water tank and numerical model studies of flow over steep smooth two-dimensional hills[END_REF], involving a 2D hill immersed in a turbulent flow. Present benchmarks showed that this method provides accurate velocity profiles both upstream and downstream from the AD device, with a maximum error of horizontal velocity reaching less than 10% at the top of the hill. The numerical simulations accurately describe the generation and stabilisation of the recirculating area, which shows that the present numerical model is able to predict the main turbulent pattern associated with this benchmark case. Together with previous results, this highlights the fact that SDM-OceaPoS can be used to predict the behaviour of most of the turbulent main flows observed in rivers and seas, without a turbine device.

12 December, 2018 • A new turbine model has been implemented based on the porous disc theory. This model can reproduce with quantitative concordance all the features measured for channel flows with actuator discs used in the recent experiments of Myers andBahaj (2010, 2012). The numerical solutions and comparisons with experimental data confirm the reliability of the PDF model in simulating turbulent flows generated by AD devices. This model can resolve the mean flow velocity deficit and also the turbulence 13 December, 2018 statistics of the wake at high Reynolds numbers for various discs, while employing relatively low grid node numbers. The 3D dynamics of the flow field characterised by the asymmetric formation of arch vortices that emerge from the disc edges and interact downstream are properly reproduced only by adding an additional acceleration term in the Lagrangian equation of particle motion. This is very promising since these mechanisms are essential for predicting the interaction between the AD device and the environment. Another important concluding remark is that this approach, first developed for application in the atmospheric frame, can be easily adapted to ocean flows, especially by changing the parametrisation of the turbulent model. Moreover,predictions obtained by SDM-OceaPoS can be used to reconstruct water particle trajectories, which can yield an accurate description of physical processes inherent to turbulent flows generated from both an arbitrary bathymetry and also by a single or several AD devices (see Di Iorio et al.).

Future work will focus on implementing an improved turbine model which includes the presence of rotating blades on a free-surface. In future, we expect this model to become a powerful tool for designing array turbines in the Chacao Channel and other similar sites, as it is ideally suited to better understanding the complex flows in natural marine environments, estimating available power and evaluating the impact of new technologies in local ecosystems.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 5 CASES INVOLVING NON-ROTATIVE TURBINE(S) 10 5.1 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: Two different approaches used to model turbulent flows -(a) Eulerian approach, (b) : Lagrangian approach.

  (see Figure 3-(left)). It involves 2D flow evolving in a 20 m long, 60 cm wide and 65 cm deep open channel with smooth ground. During experimental tests, the water is stably recirculated with maximum discharge reaching 80 l/s, which generates a turbulent flow with a Reynolds number of about Re 10 4 . The friction velocity is assessed at 0.44 cm/s.

Figure 3 :

 3 Figure 3: Sketch of the experiment from Nezu and Rodi (1986) (left). Theoretical average velocity profiles in each sub-layer of a turbulent 2D flow (left) -y + = yU * /ν, U + = U/U * , with U * = 0.44 cm/s.

Figure 4 :

 4 Figure 4: Comparisons between numerical predictions (•), theoretical predictions (-) and measured (++) mean velocities (a), averaged horizontal (b) and vertical (c) fluctuations velocities in intermediate I and II outer sub-layers flow -Y + = yU * /ν, U + = U/U * , with U * = 0.44 cm/s, and h = 10 cm.

Figure 5 :

 5 Figure 5: (a) Visualisation of the considered flow with hill (Almeida et al., 1993)-(b) Zoom of the downhill zone and illustration of the recirculation zone generated by the presence of the hill[START_REF] Loureiro | Water tank and numerical model studies of flow over steep smooth two-dimensional hills[END_REF].

Figure 6 :

 6 Figure 6: Stream lines estimated from experimental measurements (Almeida et al., 1993; Loureiro et al., 2007) (a) and calculated by SDM-OceaPoS (b) -The maximum input speed is U 0 = 2.147 m/s, and the height of the hill is h = 28 mm.

Figure 7 :

 7 Figure 7: Vertical profiles of horizontal velocities at different x positions of the domain (x = 0 m -corresponds to the position of the center of the hill) -(-): Numerical predictions, (•): Measurements from Almeida et al. (1993) (h = 28 mm).

Figure 8 :

 8 Figure 8: Experimental set-up of the research carried out by Myers and Bahaj (2010) -(a) Single disc configuration, (b) Three-discs configuration.

Figure 9 :

 9 Figure 9: Vertical profiles of turbulence intensity at different downstream positions of the turbine -(-): Numerical predictions from SDM-OceaPoS, (•): Measurements, (XX) : Numerical prediction from Johnson et al. (2014) -(a) : C T = 0.62, (b) : C T = 0.91, (c) : C T = 0.99 -Position x = 0 corresponds to the position of the disc center.

Figure 10 :

 10 Figure 10: Comparison of speed deficits of a porous disc of Ct = 0.86 predicted by SDM-OceaPoS (up) and measured in channel (Myers and Bahaj, 2010) -Downstream distance is defined as (X -X D )/D where X D is the center disc position and D the diameter of the porous disc.

Figure 11 :

 11 Figure 11: Comparison of deficits in the case of three porous discs spaced by 1.5D -(top left): measured in channel (Myers and Bahaj, 2012), (bottom left): Numerical predictions from SDM-OceaPoS, and (bottom right) : vertical profiles of velocity, (•) : Experimental, (-) : SDM-OceaPoS.

 December, 2018

December, 2018

$ This work was supported by Marine Energy Research and Innovation Center (MERIC) project L1.6 Advanced modelling for marine energies. This paper is a part of the MERIC -2018 Report.