E Dombre
email: edombre@gmail.com

J C Harris

M Benoit
email: benoit@irphe.univ-mrs.fr

D Violeau
email: damien.violeau@edf.fr

C Peyrard

A 3D parallel boundary element method on unstructured triangular grids for fully nonlinear wave-body interactions

Keywords: Nonlinear wave-structure interaction, Oshore structures, Ocean engineering, Boundary element method

This paper presents the development and validation of a three-dimensional numerical wave tank devoted to studying wave-structure interaction problems. It is based on the fully nonlinear potential ow theory, here solved by a boundary element approach and using unstructured triangular meshes of the domain's boundaries. Time updating is based on a second-order explicit Taylor series expansion. The method is parallelized using the Message Passing Interface (MPI) in order to take advantage of multi-processor systems. For radiation problems, with cylindrical bodies moving in prescribed motion, the free-surface is updated with a fully Lagrangian scheme, and is able to reproduce reference results for nonlinear forces exerted on the moving body. For diraction problems, semi-Lagrangian time-updating is used, and reproduces nonlinear eects for diraction on monopiles. Finally, we study the nonlinear wave loads on a xed semi-submersible structure, thereby illustrating the possibility to apply the proposed numerical model for the design of oshore structures and oaters.

Introduction

The numerical modeling of fully nonlinear interactions between oating structures and waves in three dimensions (3D) is of great importance for the design of ocean engineering structures such as oshore wind turbines or wave energy converters, as realistic sea states may cause nonlinear motions of the structure. The problem is often addressed by means of the fully nonlinear potential ow (FNPF) approach, and has had broad success for both radiation and diraction problems. This paper presents a variation of the 3D model of Grilli et al. [START_REF] Grilli | A fully non-linear model for three-dimensional overturning waves over an arbitrary bottom[END_REF], focusing on working with surface-piercing bodies with arbitrary geometry. Additionally, the code is parallelized to work on modern computer clusters. Notably, we continue to use Taylor series expansions for the time-stepping, shown by Machane and Canot [START_REF] Machane | High-order schemes in boundary element methods for transient non-linear free surface problems[END_REF] to be faster than Runge-Kutta for the same accuracy. In the approach of Grilli et al., however, they used higher-order elements on a structured grid, which they refer to as mid-interval interpolation (MII). Unfortunately, this is best suited for simple wave propagation problems, whereas complex surface piercing objects may not always be well-suited for structured grids. As a result, this is reformulated for unstructured 3D meshes; although often structured high-order grids will be more accurate, we believe that it will be important to have the capability to consider any mesh, as increasingly complex oshore structures may not be easily adapted to a structured or block structured mesh.

For the present study, the general theory is presented in Section 2, and the discrete equations are detailed in Section 3, including the derivation of the time-stepping scheme, the assembly of the BEM system matrix, and the representation of rigid body surfaces in our numerical model. Section 4 shows applications including the verication and the validation of the numerical model with a bottom-mounted cylinder, a truncated cylinder and nally a complex structure representing a simplied semi-submersible oater.

Concluding remarks follow in Section 5.

Mathematical model

We assume the uid to be incompressible and inviscid and the ow to be irrotational.

We thus introduce a velocity potential φ which satises mass conservation, expressed as Laplace's equation within the entire uid domain Ω f (t). We assume that the boundary of the uid domain ∂Ω f (t) is divided into four parts, on which dierent types of boundary conditions can be applied, ∂Ω f (t) = Γ f (t) ∪ Γ c (t) ∪ Γ b ∪ Γ l (t) described later, including the free-surface, Γ f , the bottom boundary, Γ b , the far-eld edges (i.e., wavemaker or sidewall boundary), Γ l , and the surface of a xed or oating body under consideration, Γ c . The bottom boundary Γ b is assumed to be time independent.

Denoting the Green's function, G(x, y) = 1 4π x-y , the fundamental solution of Laplace's equation in 3D (i.e., R 3 ), the velocity potential obeys the following boundary integral equation (BIE), for every point, x, on the boundary: c(x, t)φ(x, t) = ∂Ω f (t) ∂φ ∂n (y, t)G(x, y)φ(y, t) ∂G ∂n (x, y) dS y [START_REF] Ferrant | Fully nonlinear interactions of long-crested wave packets with a three dimensional body[END_REF] where the function c(x, t) denotes the inner solid angle seen from the boundary (eld) point x, and y is taken to be the source point on the boundary, ∂Ω f (t).

On the free surface Γ f (t), the kinematic and dynamic boundary conditions state that:

     φ t (x, t) = -gz - 1 2 ∇φ(x, t) • ∇φ(x, t) for x ∈ Γ f (t) dx dt = ∇φ(x, t) for x ∈ Γ f (t) (2)
The time-integration of these equations is described in the next section.

On the solid boundary of the oating body Γ c (t), we specify a free-slip condition, which expresses the normal derivative of the potential equal to the normal component of the body velocity on that boundary:

φ n (x, t) ≡ ∂φ ∂n (x, t) = v b (x, t) • n b (x, t), ∀x ∈ Γ f (t) (3)
where n b (x, t) denotes the unit normal vector pointing inward to the solid surface Γ c (t), at point x, and v b (x, t) is the body velocity. This condition remains valid on the xed bottom and lateral boundaries, Γ b ∪ Γ l , using a zero velocity, i.e., ∂φ ∂n = 0.

To avoid the evaluation of the time derivative of the potential by use of a nite dierence scheme, we apply the same BIE technique for computing φ t . Indeed φ t satises the same eld equation, and requires the associated boundary conditions. Following Dombre et al. [START_REF] Dombre | Simulation of oating structure dynamics in waves by implicit coupling of a fully non-linear potential ow model and a rigid body motion approach[END_REF], the Neumann boundary condition satised by φ t on Γ c (t) is expressed as:

φ tn (x, t) ≡ ∂φ t ∂n (x, t) = dn dt • (v b (x, t) -∇φ) + a b (x, t)•n -v b (x, t)• (∇∇φ•n) (4)
with a b (x, t) the solid acceleration vector at the position x and time t.

Numerical scheme

Boundary Element Discretization

At each time-step, we solve the BIE problems associated to φ and to φ t by using an isoparametric BEM with at triangles. The whole set of boundaries of the domain is meshed with non-overlapping triangles. On each triangular element Γ k , we assume the eld variables and the geometry to have linear variations, described as:

         φ k (ξ 1 , ξ 2 ) = {j;xj ∈V k } φ j N j (ξ 1 , ξ 2 ) x k (ξ 1 , ξ 2 ) = {j;xj ∈V k } x j N j (ξ 1 , ξ 2 ) (5)
where V k is the set of the vertices of Γ k and (ξ 1 , ξ 2 ) denotes the co-ordinates in the reference element Γ ξ . The functions N j are the so-called shape functions, i.e., N

1 (ξ 1 , ξ 2 ) = ξ 1 , N 2 (ξ 1 , ξ 2 ) = ξ 2 , and N 3 (ξ 1 , ξ 2 ) = 1 -ξ 1 -ξ 2 .
Using a collocation method, we write that for any x i belonging to the discrete boundary of the uid domain Γ j = ∂Ω j at time t j , we have:

c(x i , t j )φ(x i , t j ) = Γj φ n (y k , t j ) G(y k , x i )dS k -φ(y k , t j ) ∂G ∂n (y k , x i ) dS k (6)
which can be, upon replacing the integral by a discrete sum over the vertices of the mesh, rewritten as:

  δ ij c i + N dof j=1 K n ij   φ j = N dof j=1 K d ij φ j n ( 7
)
with N dof the number of vertices of the mesh. Adopting notations used in [START_REF] Grilli | A computer program for transient wave run-up[END_REF], we can show that the sub-matrices of this system are dened as:

K n ij = k∈Sj Γ ξ N l k (j) (ξ) ∂G ∂n (x k (ξ), x i )J k (ξ) dξ K d ij = k∈Sj Γ ξ N l k (j) (ξ)G(x k (ξ), x i )J k (ξ) dξ (8)
where l k (j) is a local index varying from 1 to 3, S j denotes the set of elements containing the node of global index j and J k (ξ) is the Jacobian of the transformation which maps Γ ξ to Γ k . l k (j) thus associates to each Γ k in S j the local index of the node x j in Γ k . We also dene the notation S k as the set of vertices belonging to the element Γ k .

When assembling the matrix of the system [START_REF] Ma | Finite element simulation of fully non-linear interaction between vertical cylinders and steep waves. Part 1: methodology and numerical procedure[END_REF], two situations may arise. If the source point x p is not belonging to the element Γ k , both K d ij and K n ij are regular and we use a quadrature formula dened by a Dunavant rule [START_REF] Dunavant | High degree ecient symmetrical Gaussian quadrature rules for the triangle[END_REF] in order to perform numerical integration. Otherwise, the source point is belonging to Γ k , i.e. there exists x j ∈ S k such that x i = x j , in which case we treat the singularity by applying a Duy transformation, similar to the polar change of coordinates presented in Grilli et al. [START_REF] Grilli | A fully non-linear model for three-dimensional overturning waves over an arbitrary bottom[END_REF].

In case of body motions, the velocity v b is non zero and the present model requires to evaluate the high-order term v b • (∇∇φ•n) in Eq.( 4), which is not trivial. This calculation may be simplied by transforming the second-derivative when computing the integral Γ b (t) Gφ tn dS y . This integral is decomposed into

I 1 = Γ b (t) Gφ (1)
tn dS y and

I 2 = Γ b (t) Gφ (2) tn dS y with φ (1) tn = dn dt • (v b -∇φ)+a b •n and φ (2) tn = v b • (∇∇φ•n).
Following Bai and Teng [START_REF] Bai | Simulation of second-order wave interaction with xed and oating structures in time domain[END_REF], the use of the Stokes formula and basic vector analysis manipulations lead to the relationship:

I 2 = Γ b (t) G (v b •∇(∇φ)) • ndS y = ∂Γ b (t) G(∇φ × v b ) • dy+ Γ b (t) G [(Ω × ∇φ) + (∇G • ∇φ)v b -(∇G • v b )∇φ] • ndS y (9)
assuming the rigid body velocity is given by v

b (x) = v G + Ω × (x -x G ) with v G the
translational and Ω the rotational velocity vectors of the rigid body. We recall that, by convention, the unit normal vector n is oriented towards the inside of the rigid body. In this respect, the line integral ∂Γ b (t) G(∇φ × v b ) • dy must be evaluated considering that the tangent vector to the waterline dy is such that the vector dy × n(y) points outside of the rigid body.

Time-stepping

As in the original 2D-FNPF code of [START_REF] Grilli | Corner problems and global accuracy in the boundary element solution of nonlinear wave ows[END_REF], we update both the position and the potential on the free surface Γ f (t) by an explicit scheme based on a second-order explicit Taylor series expansion. In this scheme, the values of the potential φ and the position x at time t i+1 may be expressed as follows:

             φ i+1 = φ(x i+1 , t i+1 ) = φ(x i , t i ) + dφ dt (x i , t i )∆t + d 2 φ dt 2 (x i , t i ) ∆t 2 2 x i+1 = x i + dx dt (x i , t i , φ i , φ i n )∆t + d 2 x dt 2 (x i , t i , φ i , φ i n ) ∆t 2 2 (10)
From system [START_REF] Mola | A stable and adaptive semi-lagrangian potential model for unsteady and nonlinear ship-wave interactions[END_REF], several numerical schemes can be devised according to the choice of the advection velocity used in the material derivative d• dt . Let v p be a velocity eld which is chosen to advect the free surface particles. When using the velocity v p for moving the free surface nodes, we obtain the following rst-order derivatives:

     dφ dt = φ t + v p • ∇φ dx dt = v p (11)
Then, if we take again the material derivative of the rst-order coecients along the velocity vector v p , we obtain:

       d 2 x dt 2 = ∂v p ∂t + ∇v p • v p d 2 φ dt 2 = dφ t dt + dv p dt • ∇φ + v p • d∇φ dt (12)
Given the free surface dynamic boundary condition φ t = -gz -1 2 ∇φ • ∇φ, we can derive:

dφ t dt = g • v p -∇φ • (∇φ t + ∇∇φ • v p ) (13)
We also have:

dv p dt • ∇φ + v p • d∇φ dt = ∂v p ∂t • ∇φ + (∇v p • v p ) • ∇φ + v p • ∇φ t + (∇∇φ • v p ) • v p (14)
Summing all of these contributions, we obtain:

d 2 φ dt 2 = g • v p + (v p -∇φ) • ∇φ t + ∂v p ∂t • ∇φ + (∇∇φ • v p ) • (v p -∇φ) + (∇v p • v p ) • ∇φ (15)
One can easily check that setting v p = ∇φ in Eq. [START_REF] Kring | Accelerated nonlinear wave simulations for large structures[END_REF] gives the Lagrangian second-order development:

d 2 φ dt 2 = g • ∇φ + ∇φ t • ∇φ + (∇∇φ • ∇φ) • ∇φ (16)
as can be found e.g. in Appendix A of [START_REF] Dold | An ecient surface-integral algorithm applied to unsteady gravity waves[END_REF].

We recall that the nonlinear conditions on the free surface read:

     ∂η ∂t = φ z -η x φ x -η y φ y ∂φ ∂t = -gη - 1 2 ∇φ • ∇φ (17)
When dealing with xed vertical cylinders or fully submerged bodies (or cases without any body), it can be useful to consider a semi-Lagrangian scheme, allowing only vertical motion of the uid particles. For this, we consider the denition v p = ∂η ∂t e z . (For the equivalent material derivative with curved structures, see Zhang and Kashiwagi [START_REF] Zhang | Application of ALE to nonlinear wave diraction by a non-wall-sided 29 structure[END_REF].)

While expanding this scheme, we rst make the assumption that the function (x, y) ∈ R 2 → η ∈ R is single-valued (i.e., the waves are not overturning), which allows to get the relationships:

∂η ∂t = ∂φ ∂z - ∂η ∂x ∂φ ∂x - ∂η ∂y ∂φ ∂y = ∂φ ∂z + n x n z ∂φ ∂x + n y n z ∂φ ∂y (18)
where n = (n x , n y , n z ) is the outward unit normal vector to the free surface. We can thus work out the following formula:

∂η ∂t = φ n n z (19)
which is simpler to evaluate numerically as the present model relies on the distribution of the variables φ and φ n on the boundary.

Discrete velocity on unstructured grids

For rst-order elements, the velocity in a triangle is computed using a relationship found in Meyer et al. [START_REF] Meyer | Discrete dierential-geometry operators for triangulated 2-manifolds[END_REF]. In the triangle number j, called T j , the gradient of any linear function is calculated as follows:

∇φ Tj = 1 2A j ((φ j,i+1 -φ j,i )n j × (x j,i -x j,i-1 ) +(φ j,i-1 -φ j,i )n j × (x j,i+1 -x j,i )) (20)
where A j denotes the area of T j , n j the unit outward normal vector to T j and the index (j, i) corresponds to the local index of the node number i located on the triangle T j . The indices (j, i -1) and (j, i + 1) refer to local indices of the nodes in the triangle T j such that the arc x j,i-1 x j,i x j,i+1 is positively oriented with respect to the local normal vector n j . According to these denitions, we can further derive the following relationships:

n j = (x j,i+1 -x j,i ) × (x j,i-1 -x j,i ) (x j,i+1 -x j,i ) × (x j,i-1 -x j,i ) A j = 1 2 (x j,i+1 -x j,i ) × (x j,i-1 -x j,i ) (21)
We then take an average of Eq. ( 20) over the 1-ring neighborhood of the vertex i, which is the set of triangles containing x i and denoted by S i (see Fig. [START_REF] Ferrant | Fully nonlinear interactions of long-crested wave packets with a three dimensional body[END_REF]). Taking into account the contribution of φ n to the gradient, we obtain an approximation of the velocity vector:

x j,i x i T j x j,i+1 x j,i-1
∇φ i = 1 j∈Si A j j∈Si A j ∇φ Tj + φ n 1 j∈Si A j j∈Si A j n j (22)
Eq. ( 20) has been, for example, used in the eld of computer graphics [START_REF] Meyer | Discrete dierential-geometry operators for triangulated 2-manifolds[END_REF] in order to evaluate geometric quantities such as the principal curvature on irregular meshes. It is

shown in [START_REF] Meyer | Discrete dierential-geometry operators for triangulated 2-manifolds[END_REF] to achieve a comparable accuracy to nite dierences schemes and has the advantage to be less sensitive to the mesh conguration.

Discrete acceleration on unstructured grids

For evaluating the terms related to second-order derivatives of φ, we try to take advantage of the previous methodology. This is made possible by expressing the secondorder tensor ∇∇φ in a local basis attached to the vertex under consideration. In matrix form, this term reads:

V i kl = (∇∇ B(x i ) φ) kl = ∇(∇φ • s k ) • s l (23)
with B(x i ) = (s 1 , s 2 , s 3 ) an orthonormal basis such that s 3 is the normal vector to the discrete surface at the node x i . We dene the vector:

β ijk = 1 2A j ((∇φ j,i+1 • s k -∇φ j,i • s k )n j × (x j,i -x j,i-1 ) +(∇φ j,i-1 • s k -∇φ j,i • s k )n j × (x j,i+1 -x j,i )) (24)
The components of the matrix V i can be calculated as:

               V i kl = 1 j∈Si A j ( j∈Si A j β ijk ) • s l , k ∈ {1, 2, 3}, l ∈ {1, 2} V i 13 = V i 31 V i 23 = V i 32 V i 33 = -V i 11 -V i 22 ( 25
)
Figure 2: Irregular triangular grid of a torus with a mean radius of R 1 = 2m and a cross-section of radius R 2 = 1m generated by the algorithm NETGEN2D in the SALOME platform [START_REF] Cea | The Open Source Integration Platform for Numerical Simulation[END_REF].

The relationships V i 13 = V i 31 and V i 23 = V i 32 are obtained by symmetry of the matrix V i which comes from the fact that its components are second-order partial derivatives. The

last relationship V i 33 = -V i 11 -V i
22 stems from the Laplace equation.

This scheme consists therefore in interpolating linearly the projections of ∇φ on the local basis B(x i ) over the 1-ring neighbourhood of x i . The outward normal vector at x i is evaluated by means of the following weighted sum:

s 3 (x i ) = j∈Si A j n j j∈Si A j (26)
with the vector n j previously dened on each triangle of S i . This vector is then scaled to dene a unit normal vector, and from this we dene two tangential vectors to the discrete surface at x i . For this purpose, we need to dene a rst tangential vector by projecting any point x j belonging to a triangle of the set S i and dierent from the point x i onto the tangential plane oriented by s 3 and passing through x i by means of the following formula:

s 1 (x i ) = x j -x i -s 3 (x i ) • (x j -x i ) x j -x i -s 3 (x i ) • (x j -x i ) (27)
A last vector of the basis B(x i ) is then given by s 2 (x i ) = s 3 (x i )×s 1 (x i ). The performance of these discrete dierential operators is nally assessed on the case of the irregular mesh of a torus with a mean radius of R 1 = 2m and a cross-section of radius R 2 = 1m. An example of discretization for the torus is represented on Fig. 2. We use the function φ = exp k z z sin(k x x + k y y) as a test case. We consider the following error indicators:

The normalized maximum error of the vector ∇φ (denoted by L ∞ in the gures)

and dened as =

max i=1..N ∇φ-∇φ ref i max i=1..N ∇φ ref i
The normalized average error of the vector ∇φ (denoted by L 1 in the gures) and dened as =

k x = 1m -1 , k y = 0.5m -1 and k z = k x 2 + k y 2 .
It can be concluded that the proposed method converges in L 1 norm on a smooth surface. In this norm, the error is scaling, in the root mean square sense, as ∆x R1 -1.940

for the velocity vector and as ∆x R1 -1.785

for the acceleration vector.

Parallel assembly of the system matrices

In order to speed-up the simulations for large-size problems, the nodes which make up the BEM mesh are subdivided into nearly equal subsets based on a domain decomposition. That is to say, the main simulation (e.g., time-stepping) is handled on a single processor, but when many processors are available, the workload for solving the Laplace equation is divided over the available processors. This is done by taking the linear system of equations shown in Eq. 7, and only computing (and storing) part of this coecient matrix on any given processor. The resulting system is then solved solved with an iterative linear solver, BiCGSTAB. The Message Passing Interface (MPI) is used to exchange data across the distributed computer cluster.

We This approach is in contrast to Bai and Eatock Taylor [START_REF] Bai | Numerical simulation of fully nonlinear regular and focused wave diraction around a vertical cylinder using domain decomposition[END_REF], for example, where they implement domain decomposition by adding extra boundary elements within the domain to create individual subdomains. They therefore only scale up to around a 10x speed-up compared to their solution on a single processor. This means that the discrete equations that are being solved can change depending on the number of processors used. Here, on the other hand, only the coecient matrix is shared, so the same equation is being considered. This is possible because the amount of information that needs to be shared between processors is small compared to the size of the overall system matrix.

Mesh deformation

Time-updating of deformable surfaces

The velocity potential φ and the position x on the free-surface are updated using Eq. [START_REF] Mola | A stable and adaptive semi-lagrangian potential model for unsteady and nonlinear ship-wave interactions[END_REF]. In order to avoid an incompatibility of velocity on the far-eld edges, the freesurface boundary conditions are corrected in a certain area near the exterior boundary of the free surface. For a cylindrical uid domain, as those which will be considered in the validation presented hereafter, we consider an inner cylindrical domain of radius R λ and an outer cylindrical domain of radius R ext (Fig. [START_REF] Bai | Fully nonlinear simulation of wave interaction with xed and oating ared structures[END_REF]). Then, the absorbing area is dened by the set of the assembling step plus the solution of the system matrix with respect to the number of processors np, invoked with MPI, with a mesh of 16,743 collocation nodes. Tn p is the CPU time of the process when using np processors.

D λ = {x = (r, θ, z) ∈ Γ f such that (r, θ) ∈ [R λ , R ext ] × [0, 2π]}. In
the set D λ , the free surface and the potential are corrected as:

z * ← z -∆tω i r -R λ R ext -R λ 2 (z -D(t)η i (r, θ, t)) φ * ← φ -∆tω i r -R λ R ext -R λ 2 (φ -D(t)φ i (r, θ, z * , t)) (28)
where η i and φ i correspond to a theoretical wave prole. For the absorption of the incident waves, other techniques like that developed by Clamond et al. [START_REF] Clamond | An ecient model for three-dimensional surface wave simulations. Part II: Generation and absorption[END_REF] could be more ecient for irregular waves. The implementation of this method would require to solve (in parallel with MPI) an additional linear system corresponding to the vertices of the free surface mesh. This is left for future works.

For the bottom mounted cylinder, we will use a rectangular domain. In this case, the absorbing beach dened above is divided into two parts. The rst part lies in front of the entrance of the tank and is dened as:

D λ,1 = {x = (x 1 , x 2 , x 3 ) ∈ Γ f such that (x 1 , x 2 ) ∈ [0, L y ] × [0, l λ ]
} where l λ is the length of the absorbing beach. In this domain, η i and φ i will be chosen to correspond to a desired incident wave prole. At the end of the domain, a second beach is dened:

D λ,2 = {x = (x 1 , x 2 , x 3 ) ∈ Γ f such that (x 1 , x 2 ) ∈ [0, L y ] × [L x -l λ , L x ]}.
In this second domain, we set η i = 0 and φ i = 0.

Time-updating of rigid surfaces

The nodes associated to the uid particles on the body surface are updated following the rigid body motion. We restrict ourselves to the forced motion of a rigid body. We recall that, with a rigid body, for any material points x A and x B belonging to the body, we have the relationship

{ d(x B -x A )
dt } B = 0 in the reference frame of the body B.

Rotation matrices around the axes of the xed global reference frame are dened in the following fashion:

R θx =   1 0 0 0 cos θ x -sin θ x 0 sin θ x cos θ x   (29)
R θy =   cos θ y 0 sin θ y 0 1 0 -sin θ y 0 cos θ y   (30)
R θz =   cos θ z -sin θ z 0 sin θ z cos θ z 0 0 0 1   (31)
At any time, the link between the position expressed in R and the one expressed in B is made with the equation:

(x(t) -x G (t)) R = R θx R θy R θz (x(t) -x G (t)) B (32)
with (x(t)x G (t)) B a time-invariant vector. We dene the xed basis R = (e x , e y , e z )

and the moving basis B = (e xb , e yb , e zb ). The latter can be obtained by the transformation R θx R θy R θz applied to the basis (e x , e y , e z ).

Merging of the meshes

As the free surface vertices are tracked in a Lagrangian manner during their motion, the meshes of the free-surface and the body need to be reconnected. Since the distance between the displaced waterline and the displaced body appears to be small, similarly to Bai and Eatock Taylor [START_REF] Bai | Higher-order boundary element simulation of fully nonlinear wave radiation by oscillating vertical cylinders[END_REF], we make a projection of the nodes located on the waterline onto the body geometry. Here we consider an application of this technique in the case of a cylindrical body surface. The following procedure is used for each node of the waterline:

1. Apply the MEL time-stepping to the position x * f :

x * f (t + ∆t) = x f (t) + dx f dt ∆t + d 2 x f dt 2 ∆t 2 2 (33)
2. Express the vector x * f in the body reference frame B:

x * f (t + ∆t) B = x G (t + ∆t) B +R -1 θz R -1 θy R -1 θx (x * f (t + ∆t) -x G (t + ∆t)) R (34)
3. Express x * f (t+∆t) B in a coordinate system suited to the shape under consideration (for a cylinder, we choose the polar coordinate system):

x * f (t + ∆t) B = x G (t + ∆t) B + r f (cos θ f e xb + sin θ f e yb ) + z f e zb (35)
4. Replace the radial coordinate r f by the radius of the cylinder R b in the reference frame B:

x f (t + ∆t) B = x G (t + ∆t) B + R b (cos θ f e xb + sin θ f e yb ) + z f e zb (36)
5. Transform the new coordinate x f (t + ∆t) B into the global reference frame R:

x f (t + ∆t) R = x G (t + ∆t) R + R θx R θy R θz (x f (t + ∆t) -x G (t + ∆t)) B (37)
For each geometric node of the waterline, we denote by x f the position vector of the node belonging to the mesh of the free surface and by x b the position vector of the node belonging to the mesh of the body surface. Prior to updating the position of the nodes for the next time-step, x f and x b are equal. An interpolation of the free surface elevation angular distribution found with the procedure described above, is made on a xed uniform angular distribution. In other words, for each double-node of the intersection, we rst express x f and x b in the basis (e xb , e yb , e zb ):

(x f -x G ) B = R b (cos θ f e xb + sin θ f e yb ) + z f e zb (x b -x G ) B = R b (cos θ b e xb + sin θ b e yb ) + z f e zb (38)
Then, for each angle θ bj , we compute the indices:

i 1 = argmin i∈[[1,N d ]] [|θ f i -θ bj |; θ bj ≥ θ f i ] i 2 = argmin i∈[[1,N d ]] [|θ f i -θ bj |; θ bj < θ f i ] (39)
with N d the number of double-nodes at the intersection of the mesh of the body surface and the mesh of the free surface. Afterwards, we interpolate linearly the values of the potential φ j and the values of the local free-surface elevation

(x f -x b ) B • e zb in the interval [θ i1 , θ i2 ].

Filtering of the waterline

When considering steep sea states, one often observes instabilities near the waterline, which eventually propagate throughout the domain. As in many publications, we thus lter, at each time step, the variables of the waterline dened previously. We apply, similarly to [START_REF] Zhou | Numerical investigation of wave radiation by a vertical cylinder using a fully nonlinear HOBEM[END_REF][START_REF] Tanizawa | The state of the art on numerical wave tank[END_REF], a ltering formula with a moving stencil of 7 points for the variables φ and (x fx b ) B • e zb along the waterline only. For the potential φ, this formula reads:

φ i = 1 32 (-φ i-3 + 9φ i-1 + 16φ i + 9φ i+1 -φ i+3 ). (40)

Remeshing step

In time-domain simulations with oating bodies, the geometry is always changing.

As a consequence it is necessary to change the mesh of the geometry. On the one hand, in order to prevent the triangles to be distorted in the vicinity of the piercing surface cylinder, a Laplace smoothing technique is applied, for the mesh of the free surface, on the projection of the nodes onto an horizontal plane. This technique has also been used by [START_REF] Ma | Finite element simulation of fully non-linear interaction between vertical cylinders and steep waves. Part 1: methodology and numerical procedure[END_REF][START_REF] Bai | Higher-order boundary element simulation of fully nonlinear wave radiation by oscillating vertical cylinders[END_REF]. On the other hand, the mesh of the body is generated with regular quadrangles divided into triangles by their diagonal. For this part of the mesh, we modify only the vertical position of the nodes to get a uniform vertical distribution.

Numerical results

4.1. Sway motion of a truncated vertical cylinder

Problem setup

We are concerned in this section with the imposed periodic motion in sway (i.e. along x 1 axis) of a truncated cylinder. The position vector of the center of gravity is given the form:

x G (t) = (x 1 (t), x 2 (t), x 3 (t)) = (A sin ωt, 0, 0) (41) 14

Rext -R λ R λ -R R dφ dt = -gz + v 2 2 -γ(x)(φ -φ ref ) dx dt = ∇φ -γ(x)(η -η ref ) dφ dt = -gz + v 2 2 dx dt = ∇φ

Mesh convergence

In this section, we study the mesh convergence of the model for one of the frequencies tested hereafter. We note that by having a damping region such as Eq. 28, as used before by Cointe [START_REF] Cointe | Numerical simulation of a wave channel[END_REF], volume conservation is not inherently guaranteed, as the freesurface boundary conditions are modied. We consider the non-dimensional wavenumber kR = 1.4 (with a cylinder radius of 1 m), which gives a wavelength λ ≈ 4.48m and a period T ≈ 1.72s. The radius of the external domain is chosen as R ext = 8m. An absorbing beach of radial length 4.5m is used. We check in this section the convergence by looking at the conservation of the discrete water volume in the wave tank. The characteristics of the meshes used for the convergence study are recalled in Table 1. 1: Truncated cylinder in sway motion along the axis Ox: characteristics of the meshes used in the volume conservation study. ∆r e (respectively ∆r i ) is the space-step of the mesh on the circumference of the outer (respectively inner) cylinder. ∆z e (respectively ∆z i ) is the space-step of the mesh on the vertical direction of the outer (respectively inner) cylinder. N F S is the number of nodes on the free surface. N body is the number of nodes on the truncated cylinder surface. l 1D is the length of the edges on the mesh of the free surface.

Mesh a

In Fig. 7, we represent the volume error dened as:

(t + ) = V (t + ) -V (t + = 0) V (t + = 0) (42)
as a function of the non-dimensional time t + = t T and the spatial discretization. For each mesh, we enforce a CourantFriedrichsLewy (CFL) condition by automatically adapting the time-step as follows:

∆t = C 0 min (i,j);i<j;(xi,xj )∈Γ f (t) 2 x i -x j √ gd ( 43
)
where we choose C 0 = 0.4, following Grilli et al. [START_REF] Grilli | An ecient boundary element method for nonlinear water waves[END_REF].

It is seen in Fig. 7 that the error V stabilizes after 6 periods and decreases with the time-step. More specically, after 7 periods, the mean relative error is around 8.10 -5

for the time-step ∆t a = 0.0122T , while it is around 4.10 -5 for the time-step ∆t b = 0.0069T ≈ 1 2 ∆t a . This shows that the error is decreasing linearly with the time-step.

Moreover we can observe that the amplitude of the high frequency oscillations decreases also with the time-step.

Force coecients

The radius of the outer cylinder is set to R ext = 2λ where λ is the wavelength expected from the linear theory at the water depth d and the angular frequency ω (Fig. 8). The entrance of the absorbing beach is located at a distance R λ = λ from the origin (0, 0, 0). This problem has been recently addressed by Zhou et al. [START_REF] Zhou | Numerical investigation of wave radiation by a vertical cylinder using a fully nonlinear HOBEM[END_REF], using a HOBEM devised for structured quadrangular meshes. Our results are compared to the linear BEM of Yeung [START_REF] Yeung | Added mass and damping of a vertical cylinder in nite-depth waters[END_REF], the second-order frequency BEM model of Teng et al. [START_REF] Teng | Simulation of second-order radiation of 3D bodies in time domain by a B-spline method[END_REF], and the nonlinear method of Zhou et al.

For the amplitude of motion under consideration, A R = 0.15, we observe on Fig. [START_REF] Shao | A Harmonic polynomial cell (HPC) method for 3D Laplace equation with application in marine hydrodynamics[END_REF] a fair agreement between the added-mass computed with our model and the results of Zhou et al. [START_REF] Zhou | Numerical investigation of wave radiation by a vertical cylinder using a fully nonlinear HOBEM[END_REF]. The discrepancy with the linear results of Yeung is reduced as the nondimensional wavenumber kR decreases. This behavior is not surprising as, for a given choice of non-dimensional amplitude A R , the associated wave steepness will increase with increasing value of kR, resulting in larger nonlinear eects.

For the linear radiation damping coecient B 11 , represented in Fig. 10, a still closer agreement with the linear theory is found. A similar agreement is observed for the nondimensional rst-order coecient of the horizontal force

F (1)
x ρAR 2 , represented on Fig. 11, except a slight discrepancy which appears at kR = 1.

Despite the apparently linear behavior of the horizontal force, this case clearly shows the interest of using a nonlinear model. As mentioned by Zhou et al. [START_REF] Zhou | Numerical investigation of wave radiation by a vertical cylinder using a fully nonlinear HOBEM[END_REF], according to a theoretical result demonstrated by Wu [START_REF] Wu | A note on non-linear hydrodynamic force on a oating body[END_REF], the vertical force oscillates at twice the frequency of the motion 2ω. As before, the discrepancy between the weakly nonlinear model of Teng et al. [START_REF] Teng | Simulation of second-order radiation of 3D bodies in time domain by a B-spline method[END_REF] and the fully nonlinear model of Zhou et al. increases with increasing values of the parameter kR. The existence of the second-order coecient F

(2) z is well captured with our model as can be seen in Fig. 12. We see that our method allows us to nd a very good agreement with the frequency analysis, for wavenumbers such that second-order Fourier coecient associated to the vertical force Fz with respect to the non-dimensional wavenumber kR for various numerical models.

- We now consider a more complicated test case with a semi-submersible structure, hereafter referred to as the Dutch Tri-oater [START_REF] Bulder | Study to feasibility of and boundary conditions for oating oshore wind turbines[END_REF], which is a proposed structure for oshore wind turbines. The shape of the oater involves three surface piercing cylinders, rigidly connected. Moreover, each cylinder is equipped at its base with a heave plate.

This geometry is a slight simplication of original the Dutch Tri-oater geometry [START_REF] Bulder | Study to feasibility of and boundary conditions for oating oshore wind turbines[END_REF],

as we do not consider here the smaller supporting struts, but they should not constitute a signicant source of hydrodynamic forces with an inviscid model. A sketch of the geometry adopted in this work, also studied by Antonutti et al. in [START_REF] Antonutti | An investigation of the eects of wind-induced inclination on oating wind turbine dynamics: heave plate excursion[END_REF], is represented on Fig. [START_REF] Machane | High-order schemes in boundary element methods for transient non-linear free surface problems[END_REF].

For the numerical wave tank, we again use a cylindrical domain such as the one coecients requested to evaluate this solution are computed at the beginning of the simulation.

For the size of the domain, trial and error showed that a choice of parameters R ext = 3λ + Λ, R λ = λ + Λ with Λ = 70m and γ 0 = 1, leads to a stable value of wave loads after a duration of 40 wave periods, except for nondimensional wavenumbers kr ≥ 2.5, for which longer time simulations are needed. On the contrary, by choosing (R λ , R ext ) = (2λ + Λ, 3λ + Λ) or (R λ , R ext ) = (λ + Λ, 2λ + Λ), that is, with an absorbing beach of radial length λ, we cannot obtain a convergence of the wave loads towards a periodic state.

We study below the diraction of monochromatic waves on the oater. The waves are generated by imposing the velocity and acceleration of fth-order Stokes waves on the far-eld lateral sides. In addition, as already mentioned earlier, the wave prole Above the mean free surface level, each cylinder is assumed to be innite. is damped out near the external boundary. With the methodology exposed above, we were able to simulate successfully the waves whose characteristics are given in Table [START_REF] Liu | Computations of fully nonlinear three-dimensional wave-wave and wave-body interactions. Part 2. Nonlinear waves and forces on a body[END_REF] of the Appendix A. We note that for the steepest waves, the convergence of the wave loads toward a periodic state is slower, which suggests that in the future, it could be worth computing the Response Amplitude Operator (RAO) curve with a xed steepness rather than with a xed wave height (as done here in the traditional manner). The meshes used in this study are always adapted to the wavelength under consideration.

The unstructured mesh of the free surface is done with the algorithm NETGEN2D in the SALOME platform [START_REF] Cea | The Open Source Integration Platform for Numerical Simulation[END_REF]. Each simulation used a mesh with a number of degrees of freedom in the range: N dof ∈ [[18000, 23000]].

Zero and rst order loads

We analyze the time-series of the forces F x and F y , and the overturning moment M y , over a window of 6 periods. Details on the time intervals considered for the Fourier analysis are given in Table 3 of Appendix A. For each frequency, we dene the RAO of the temporal signal F as the ratio between the rst-order Fourier coecient F (1) and the wave height H, and we assume that the density of the water is ρ f = 1025 kg/m 3 .

The RAO of the longitudinal force, F x , agrees very well with the open-source linear wave model NEMOH [START_REF]Open-source code NEMOH[END_REF] (Fig. 18). In particular, we nd that the oscillations of the curve are well reproduced for the small periods. Similar agreements are found for the vertical force (Fig. 19), using reference values published in [START_REF] Antonutti | An investigation of the eects of wind-induced inclination on oating wind turbine dynamics: heave plate excursion[END_REF]. There are some deviations between the two results, even for these moderate wave amplitudes, at the peaks, as one might expect between linear and nonlinear models. Similarly, the nonlinear model described here predicts a lower maximum overturning moment, as compared to the linear model NEMOH (Fig. 20).

In addition, we show in Fig. [START_REF] Bai | Simulation of second-order wave interaction with xed and oating structures in time domain[END_REF] the horizontal drift force, dened as the zero-order coecient F

x of the Fourier series F x , with respect to the wave period T . Our results are compared to a frequency model computing the Quadratic Transfer Functions described in [START_REF] Philippe | Introducing Second Order Low Frequency Loads in the Open-Source Boundary Element Method Code Nemoh[END_REF]. Except for the smallest periods, a very good agreement with the weakly nonlinear model is observed. For the highest frequencies tested in this study, our semi-Lagrangian scheme exhibits important deviations from the theory. This calculation provides further condence in the ability of our numerical model to represent nonlinear features of the wave-structure interaction problem.

As an illustration of the importance of nonlinear eects, we choose to represent on Fig. [START_REF] Grilli | Corner problems and global accuracy in the boundary element solution of nonlinear wave ows[END_REF], the ratio of second-order to rst-order coecients for the quantities of interest:

F x , F z and M y . Whereas this ratio does not exceed 20% for F x , it may reach 60% for F z and M y , for the smallest periods of this study. This shows that nonlinearities become important for these shorter wavelengths with the Dutch Tri-oater.

Conclusions

We presented in this paper an implementation of a fully nonlinear potential wave model to simulate wave-structure interactions using unstructured triangular meshes, important for being able to handle future industrial applications with arbitrary problem geometry. The assembling of the system matrix is made with an ecient use of parallelization on distributed computer systems. Two time-stepping schemes, based on discrete derivatives with rst-order shape functions, are derived. The accuracy of the whole algorithm could be easily enhanced by using high-order elements.

The model is applied to various problems involving surface piercing cylinders. The forces on a truncated cylinder in nite water depth, subjected to a prescribed sway motion, using a fully Lagrangian motion of the free surface mesh, is found to agree reasonably well with reference results. Similarly, forces on a bottom-mounted vertical cylinder resulting from diraction of regular waves, using a semi-Lagrangian time-stepping scheme, capture higher-order nonlinear eects (i.e., up to the fourth harmonic of the horizontal force).

Finally, in order to show the potential of the method in dealing with complex structures, we also successfully compute both linear diraction loads and nonlinear drift forces for a geometry inspired by the Dutch Tri-oater [START_REF] Bulder | Study to feasibility of and boundary conditions for oating oshore wind turbines[END_REF]. Other important features, such as improving computational speed (i.e., with the fast multipole method), such as started by Harris et al. [START_REF] Harris | Fast integral equation methods for fully nonlinear water wave modeling[END_REF], will be considered in upcoming works, as well as using improved accuracy (i.e., cubic B-spline elements), or including more physics (i.e., coupling to Navier-Stokes solvers). Moreover, the case of arbitrary geometries for the rigid body needs to be addressed. |F (2) x | |F |F (2) z | |F |M (2) y | |M

Figure 1 :

 1 Figure 1: Sketch of the set of elements S i , also called the 1-ring neighborhood of the vertex x i .

Figure 3 :

 3 Figure 3: Maximum normalized error (left) and average normalized error (right) for the velocity vector v associated to the potential φ = exp(kzz) sin(kxx + kyy) with kx = 1m -1 , ky = 0.5m -1 and kz = kx 2 + ky 2 on a torus.

Figure 4 :

 4 Figure 4: Maximum normalized error (left) and average normalized error (right) for the acceleration vector a associated to the potential φ = exp(kzz) sin(kxx + kyy) with kx = 1m -1 , ky = 0.5m -1 and kz = kx 2 + ky 2 on a torus.

  test the eciency of our parallelization by computing the solution of a mixed Boundary Value Problem for the velocity potential φ(x, y, z) = x. The geometry is a box of size L x × L y × L z = 1m × 1m × 1m. On the top surface, we impose a Dirichlet boundary condition φ = x. On the remaining surfaces, we impose a Neumann boundary condition φ n = n x , n x being the x-component of the outward normal vector to the box. We test the solution process for a moderate grid of 16,743 collocation nodes with dierent numbers of processors n p ∈ {1, 4, 16, 64}. The results are represented on Fig. 5.

Figure 5 :

 5 Figure 5: Speed-up T np =1 Tn p

Figure 6 :

 6 Figure6: Sketch of the decomposition of the domain with respect to the free surface boundary conditions, noting boundary conditions for fully Lagrangian time-updating. In this case, Laplace smoothing is used to maintain a well formed mesh on the free-surface.

Figure 7 :

 7 Figure 7: Truncated cylinder in sway motion along the axis Ox: computation of the relative volume error V (t + ) for various spatial and time discretizations.

kR ≤ 0. 8 .

 8 For decreasing values of kR, the present results exhibit a closer convergence towards the results of the second-order model in comparison with the model of Zhou et al. For larger kR, our nonlinear model deviates from the weakly nonlinear model and predicts higher values of F (2) z , which is not the case of the model of Zhou et al.

Figure 8 :

 8 Figure 8: Computational domain (left panel) and close-up of grid near cylinder (right panel) for test cases of waves moving cylinder.

Figure 9 :

 9 Figure 9: Truncated cylinder in sway motion along the axis Ox: computation of the added-mass with respect to the non-dimensional wavenumber kR for various numerical models.

Figure 10 : 15 Figure 11 : 8 F ( 2 ) z ρgA 2 R

 101511822 Figure 10: Truncated cylinder in sway motion along the axis Ox: computation of the radiation damping coecient with respect to the non-dimensional wavenumber kR for various numerical models.

Figure 12 :

 12 Figure12: Truncated cylinder in sway motion along the axis Ox: second-order Fourier coecient associated to the vertical force Fz with respect to the non-dimensional wavenumber kR for various numerical models.

Figure 13 : 15 ,

 1315 Figure13: Truncated cylinder in sway motion along the axis Ox: zero-order Fourier coecient associated to the vertical force Fz with respect to the non-dimensional wavenumber kR for various numerical models.

Figure 14 :

 14 Figure 14: Truncated cylinder in sway motion along the axis Ox: rst-order Fourier coecient associated to the overturning moment My with respect to the non-dimensional wavenumber kR for various numerical models.

Figure 15 :

 15 Figure 15: Computational domain (upper panel) and close-up of grid near cylinder (lower panel) for test cases of waves interacting with a xed cylinder, in the case of kR = 0.245 and kA = 0.10.

4. 3 .

 3 Diraction on a modied xed Dutch Tri-oater platform 4.3.1. Setup of the model

  represented in Fig. (6). This time, the reference free surface elevation η ref and the reference potential φ ref correspond to the solution of the fth-order Stokes theory. The Design draft (m) 12.0 Column centre-to-center spacing (m) 68.0 Column diameter (m)

4 )Figure 16 :

 416 Figure16: First through fourth-order Fourier coecients associated to the time-series of Fx for a bottommounted cylinder with a period T = 2.03s, kR = 0.245 and the water depth kd = π.

Figure 17 :

 17 Figure 17: Geometry of the modied Dutch Tri-oater taken into account in the potential ow model.

  linear model Nemoh Antonutti et al. (2014) present FNPF model, H = 2m

Figure 18 :

 18 Figure 18: Nonlinear diraction of a wave of height H = 2 m around the modied Dutch Tri-oater for various periods T with a water depth of d = 50 m: RAO of the horizontal force Fx.

Figure 19 :

 19 Figure 19: Nonlinear diraction of a wave of height H = 2 m around the modied Dutch Tri-oater for various periods T with a water depth of d = 50 m: RAO of the vertical force Fz.

Figure 20 :

 20 Figure 20: Nonlinear diraction of a wave of height H = 2 m around the modied Dutch Tri-oater for various periods T with a water depth of d = 50 m: RAO of the overturning moment My.

Figure 21 :

 21 Figure 21: Nonlinear diraction of a wave of height H = 2 m around the modied Dutch Tri-oater for various periods T with a water depth of d = 50 m: horizontal drift force Fx (0) .

Figure 22 :

 22 Figure 22: Nonlinear diraction of a wave of height H = 2 m around the modied Dutch Tri-oater: ratio between second-order to rst-order coecients for the quantities Fx, Fz and My.

Table

  Mesh b Mesh c Mesh d Mesh e

	2πRext ∆re 2πR ∆re h ∆ze B ∆zi	100 20 4 8	150 30 6 12	200 40 6 14	250 50 8 18	300 60 10 20
	min l 1D [m]	0.313	0.209	0.157	0.126	0.105
	max l 1D [m]	0.505	0.335	0.251	0.201	0.168
	N F S	856	3141	3286	12437	12576
	N body	210	535	760	1116	1940

Table 2 :

 2 Geometric parameters of the Dutch Tri-oater chosen as in the study[START_REF] Antonutti | An investigation of the eects of wind-induced inclination on oating wind turbine dynamics: heave plate excursion[END_REF].

	|f (1) ) |/(ρgAR 2	0.00 6.0 6.2 6.4 6.6 6.8 7.0 7.2	0.05	0.10	0.15	0.20	BEM Exp.	0.25
	|f (2)	0.00 0.0 0.1 0.7 0.6 0.5 0.4 0.3 0.2	0.05	0.10	0.15	0.20		0.25

Table 3 :

 3 Nonlinear diraction for a wave height H = 2 m around the modied Dutch Tri-oater: characteristics of the waves simulated with respect to the non-dimensional wavenumber. The parameter r denotes the radius of the heave plates and is set to r = 9 m. d is the water-depth, such that H d = 0.04. ω is the angular wave frequency. k is the wave number. T is the wave period. t + 1 = t 1 T and t + 2 = t 2 T with [t 1 , t 2 ] the time interval on which the Fourier analysis is performed.

Acknowledgements

The contribution of E. Dombre to this work was partly funded by the French ANRT Agency (CIFRE agreement # 2011-1724) and by the Fondation des Ponts through the Barré de Saint-Venant research fund. The work of J. C. Harris and M. Benoit was partially funded as a part of the French ANR (Agence Nationale de la Recherche), project ANR11-MONU-018-01 MONACOREV.

On Fig. 13, the non-dimensional dierence between the zero-order coecient F (0) z and the hydrostatic force ρgV b , V b being the discretized wet volume of the body, is represented. Our method predicts very well the load for kR = 0.4 and the trend of this quantity is correctly reproduced.

Finally, we represent on Fig. 14 the variation of the rst-order coecient associated to the overturning moment M (1) y ρAR 3 . Our results are very close to those of Zhou et al. [START_REF] Zhou | Numerical investigation of wave radiation by a vertical cylinder using a fully nonlinear HOBEM[END_REF] for wavenumbers kR smaller than kR = 1. For larger values of kR, we observe with our model a clear deviation from the linear theory, suggesting a signicant inuence of nonlinear eects at high frequencies, although this is contrary to the results of Zhou et al..

Diraction of long waves on a bottom-mounted cylinder

In this section, we analyze the diraction of a long wave on a bottom mounted vertical cylinder of radius R = 0.25m (Fig. 15). The wavenumber k is chosen such that kR = 0.245, which gives for linear theory a wave period T = 2.03s. Simulations are repeated for 6 dierent wave amplitudes such that kA = {0.025, 0.05, 0.075, 0.10, 0.125, 0.150}.

The signal is analyzed with a discrete Fast Fourier Transform on the time-interval t + = t T ∈ [START_REF] Bai | Higher-order boundary element simulation of fully nonlinear wave radiation by oscillating vertical cylinders[END_REF][START_REF] Wu | Simulation of nonlinear interactions between waves and oating bodies through a nite-element-based numerical tank[END_REF]. For this case, the wavelength in deep-water may be estimated as λ ≈ 6.41m. wave conditions are imposed using the fth-order Stokes theory. Fourier coecients of the potential up to fourth-order are represented on Fig. 16, where they are compared to the experiments carried out by Huseby and Grue [START_REF] Huseby | An experimental investigation of higher-harmonic wave forces on a vertical cylinder[END_REF]. Results match the experiments well, with a slight dierence seen in the second harmonic, but this is quite consistent with other fully nonlinear computations, such as those by Ferrant [START_REF] Ferrant | Fully nonlinear interactions of long-crested wave packets with a three dimensional body[END_REF] and Shao and Faltinsen [START_REF] Shao | A Harmonic polynomial cell (HPC) method for 3D Laplace equation with application in marine hydrodynamics[END_REF]. Obtaining such an agreement for harmonics of the wave force signal up to the fourthorder clearly demonstrates the nonlinear capabilities of the proposed modeling approach.

The domain is a box with the dimensions

Note that contrary to Ferrant [START_REF] Ferrant | Fully nonlinear interactions of long-crested wave packets with a three dimensional body[END_REF], we do not decompose the velocity potential into an incident component plus a diraction component.

Note that the simulations performed here lead to simulated values for the forces which vary smoothly in time. The free surface pattern, can be noisy near the waterline without the ltering described by Eq. 40. As this test case is for a xed body, we do not need to invoke the Laplace smoothing or regridding described earlier, simply this ltering along the waterline. Recent experimental tests of interactions between regular waves and vertical cylinders [START_REF] Swan | The interaction between steep waves and a surface-piercing column[END_REF], however, suggest that waves near the waterline may physically be causing local breaking, or be damped by viscous eects. Calibrating the amount of dissipation required will be further investigated in future work; here we simply remove high-frequency waves, similar to Longuet-Higgins and Cokelet [START_REF] Longuet-Higgins | The deformation of steep surface waves on water. I. A numerical method of computation[END_REF], which does not signicantly aect waves which are well-resolved.


Patrick Savelli

Vincent Savaux

Pauline Desnos

Ali Zeineddine

Matthieu Kanj

Christophe Delacourt

Flexible Multi-standard Digital Front-End for LPWA Technologies

This paper proposes a flexible multi-standard digital frontend (DFE) hardware architecture designed for three main low-power wide-area (LPWA) technologies: LoRa, Sigfox, and narrowband-Internet of things (NB-IoT). We demonstrate the feasibility of an unified DFE architecture that fits the requirements of these LPWA technologies. The proposed DFE architecture has been implemented on a Spartan-6 field programmable gate array (FPGA) within a base station receiver, and tested using a universal software radio peripheral (USRP) based platform. We show that the expected performance can be achieved with low hardware complexity, in terms of memory and logic requirements.

Introduction

The development of the Internet of things (IoT) market has led to a great diversity of wireless connectivity solutions. This paper addresses the low-power wide-area (LPWA) solutions, targeting IoT networks with coverage in the kilometer order, using low-power end devices having a battery life autonomy of multiple years. Various LPWA technologies and standards have emerged, with different air interfaces : LoRa, Sigfox, NB-IoT, LTE-M, DASH7, Weightless, NB-Fi, RPMA, MY-THINGS, and others [START_REF] Technologies | Internet of Things (IoT) (Poster)," Table 7. FPGA DFE complexity figures Registers LUTs RAM DSP Available[END_REF][START_REF] Ayoub | Internet of Mobile Things: Overview of LoRaWAN,DASH7, and NB-IoT in LPWANs standards andSupported Mobility[END_REF][START_REF] Mekki | A comparative study of LPWAN technologies for large-scale IoT deployment[END_REF]. In this context where multiple options are possible, there is a need for technical components that are flexible enough to support different technologies. In this paper, we focus on the main LPWA technologies that have been widely adopted and deployed by operators, as part of large ecosystems [START_REF]LPWAN Market Report 2018-2023[END_REF]: LoRa, Sigfox and NB-IoT The LoRa technology introduced by Semtech uses a chirpbased patented modulation scheme, while the network can be deployed by any third-party entity based on the LoRaWAN open specifications [START_REF]LoRa Alliance Technical Committee[END_REF] maintained by the LoRa Alliance1 . Sigfox technology proposed by Sigfox S.A. also uses the unlicensed spectrum, and is based on the differential binary phase-shift keying (D-BPSK) and Gaussian frequency shift keying (GFSK) modulation schemes. Moreover, the network deployment of this technology is controlled by Sigfox. NB-IoT, initially published in Release 13 of the third generation partnership project (3GPP) organization 2 , is based on cellular mobile network, using resources within frequency bands allocated to LTE (inband mode), or dedicated to NB-IoT (standalone mode). This diversity increases the LPWA network deployment cost in case each technology is independently deployed. To solve the problems of a heterogeneous network, and to efficiently exploit virtualization techniques, the interoperability between these LPWA technologies should be maximized. On the radio front-end level, this is achieved by providing LPWA gateways that are capable of supporting the largest number of technologies. This type of gateways requires a digital front-end (DFE) implementation that provides the flexibility needed to meet the requirements of these different technologies. This paper proposes a generic approach to build a flexible multi-standard DFE for LPWA technologies. A focus on the hardware implementation of the DFE is presented. A common architecture for the three considered LPWA solutions is proposed, whilst respecting the constraints imposed by their specifications. The rest of the paper is organized as follows. In section 2, the flexible DFE architecture is presented based on the requirements of the different LPWA technologies. Section 3 presents the hardware implementation of this flexible DFE. Section 4 concludes this work.

2 Multi-standard DFE design approach

Overview

A preliminary analysis phase of the three protocols allowed us to design the DFE around three blocks shown in Fig. 1, and described as follows:

• Rate adaptation: the sampled signal at the output of the analog to digital converter (ADC) is sub-sampled from sampling frequency f 0 (ADC frequency) to frequency f 1 , with f 1 ≤ f 0 . A very efficient and low-complexity rate adaptation is implemented based on [START_REF] Zeineddine | Reconfigurable newton structure for sample rate conversion[END_REF][START_REF] Zeineddine | Efficient arbitrary sample rate conversion for multistandard digital front-ends[END_REF].

• Filtering: A low-pass finite impulse response (FIR) filter operating at frequency f 1 , allows to reject neighbouring interferers and blocking signals, as well as performing a pre-equalization in the band-pass area. Since the computational cost of the filtering step is directly related to the length of the filter (i.e. its number of taps), it is designed as short as possible while respecting the requirements. We used our own FIR filter design tool3 to generate equiripple filters [START_REF] Paquelet | On the symmetry of FIR filter with linear phase[END_REF]. • Channels demultiplexing: the aim of the channels demultiplexing block is to split the whole band into different channels, e.g. channels of 125 kHz in LoRa (i.e. f 2 ≥ 125 kHz) or of 180 kHz in NB-IoT (i.e. f 2 ≥ 180 kHz). It must be noted that there is no default channel split defined in Sigfox, so we propose to apply several channel demultiplexing blocks in parallel, using the same input signal with a different frequency shift, so as to use a similar architecture.

These blocks are common to the three technologies, and only require specific parameters to be adapted to each one of these. It then allows for rapid hardware developments of flexible DFEs. Fig. 2 shows an example of frequency response for the rate adaptation, the low-pass filter, and the combination of the latter (referred to as "global"). In this example, the DFE is adapted to the reception of NB-IoT within an LTE band, considering a 5MHz bandwidth. The ADC frequency is set to f 0 = 60.8 MHz and the output frequency of the rate adaptation and filtering stage is f 1 = 7.68 MHz, corresponding to the sampling frequency of the 5 MHz LTE band. The filter has been designed for a 100dB rejection and a ripple below ±1 dB. Furthermore, it can be observed that the filter applies an equalization in the pass-band to counterbalance the attenuation due to the rate adaptation block.

Design constraints

We provide in this section the list of RF performance parameters that impact the DFE design for the targeted air interfaces.

LoRa: LoRa networks are operated in unlicensed bands using a proprietary physical layer (PHY) designed by Semtech. The LoRa PHY main characteristics are listed in Table 1. The European industrial, scientific and medical (ISM) frequency band used for LoRa operation is indicated for the sake of example, but other unlicensed frequency bands are also used according to the geographical area [START_REF]Lorawan®regional parameters[END_REF].

Based on these characteristics, we have derived the design parameters of our DFE for the LoRa mode of operation as indicated in Table 2. The provided parameters have been carefully adjusted to achieve high channel selectivity, assuming the parallel reception of 8 LoRa RF channels with 200kHz channel raster. The output rate of the DFE in this case is 250kSamples/s, which corresponds to a 1MB/s per channel data rate when using 16 bits for each of the I and Q paths.

Sigfox: Sigfox has released their proprietary radio specifications in [START_REF]Sigfox connected objects: Radio specifications[END_REF], including the air interface description and the different regional parameters. fixed channel raster, and must be able to demodulate an uplink modulated signal present at any frequency within the usable bandwidth. In addition, the average RF channel frequency drift observed during a transmission is as high as ±20 Hz/sec due to the reference frequency source drift, and the typical burst duration exceeds 1 seconds. This constraint needs to be addressed when designing the filtering and channel demultiplexing in the DFE In order to fulfill the above constraints, our proposal is to consider 300Hz wide channels every 200Hz, using two parallel filtering and channels demultiplexing entities. Each parallel entity processes a different version of the received signal shifted by a 200Hz offset. In total, there are 750 channels at the output of the DFE to cover the 150kHz band. The DFE parameters set to implement this Sigfox configuration are provided in Table 4. Note that at the DFE output, the channel bandwidth is 300Hz, so it is required to further process the output of the DFE to demodulate the 100 bauds D-BPSK signal.

NB-IoT: NB-IoT is a cellular LPWA technology introduced in the Release 13 of the 3GPP specifications. In this contribution, we target the NB-IoT inband operation mode, where some LTE resource blocks of 180 kHz bandwidth are allocated to NB-IoT. For instance, four specific resource blocks are defined when NB-IoT is transmitted in the 5MHz LTE band. Our proposal for the NB-IoT operation mode is to use the DFE for the extraction of 180kHz wide physical resource blocks (PRB) within the LTE band. Tables 5 and6 summarize the NB-IoT constraints and the corresponding configuration, respectively, considering the LTE 5MHz bandwidth, and inband mode.

FPGA implementation and results

Table 7 summarizes the Rx DFE hardware resources utilization for each of the targeted IoT standards. Four key resources are detailed: registers used to save temporary results and variables, look-up tables (LUTs) that determine the combinatorial behaviour, memories (RAMs) to store parameters such as the filter taps, and digital signal processors (DSPs) that are specialized units for multiplication/addition computation. This table shows that any of the LoRa, Sigfox or NB-IoT DFEs fits in the small cost-effective B205mini board. The different technologies use similar resources, except for the RAM requirement that is much higher for the Sigfox mode. This is mainly due to the number of channels to be processed in parallel, requiring a large amount of filter taps. The RAM requirement is also higher in case of NB-IoT compared to LoRa. The worst case in NB-IoT inband mode is indeed the 20MHz LTE bandwidth case, where 18 PRBs are dedicated to NB-IoT and need to be extracted in parallel. This impacts the number of filter

Figure 1 .

 1 Figure 1. General DFE block diagram, including rate adaptation, filtering, and channels demultiplexing (demux.).

Figure 2 .

 2 Figure 2. Filter shapes at each step and global filtering. Note that FIR filter compensates the rate adaptation attenuation in pass band.

Table 3

 3

	summarizes the radio

It is also important to note that the channel bandwidth occupied by the Sigfox signal (100 Hz) is low compared to the channel frequency inaccuracies of the transmitter and receiver ends. For this reason, the Sigfox gateway receiver cannot assume a

Table 1 .

 1 LoRa PHY characteristics

		Europe
	Frequency band	863-870MHz
	Duplex mode	half duplex
	Channel bandwidth	125kHz
	Channel raster	200kHz
	Number of channels 3 mandatory, typically 8
	Gateway Tx power	20dBm
	Modulation	Chirp Spread Spectrum
	Spreading Factor	7-12
	Coding rate	4/5, 4/6, 4/7, 4/8

Table 2 .

 2 DFE configuration for LoRa, 8 channels

	DFE design parameters
	ADCs output Rate f 0	60.8MHz
	Rate adapt. output rate f 1	1.6MHz
	Channels demultiplexing 250kSamples/s
	output rate f 2	per channel
	number of filter taps	500
	Achieved channel selectivity
	carrier offset	rejection
	±85kHz	75dB

Table 3 .

 3 Sigfox PHY characteristics

	Europe

Table 4 .

 4 DFE configuration for Sigfox

	DFE design parameters
	ADCs output Rate f 0	60MHz
	Rate adapt. output rate f 1	150kHz
	Channels demultiplexing 400Samples/s
	output rate f 2	per channel
	number of filter taps	24000
	Achieved channel selectivity
	carrier offset	rejection
	±200Hz	60dB
	3 Hardware Implementation
	3.1 Test platform
	Our gateway experimentation platform is composed of an off-
	the-shelf NI USRP B205mini board. This software defined
	radio (SDR) platform contains the RF transceiver and an Xilinx
	Spartan-6 FPGA with a 83.3 × 50.8 × 8.4 mm form factor. The
	DFE function has been implemented in the user programmable
	FPGA. The receiver also includes, prior to the B205mini RF
	receiver, an analog front-end board of our own design, that
	includes a low noise amplifier (LNA) achieving 30 dB gain
	with 3 dB noise figure (NF) in the 868 MHz ISM band and the
	LTE band 8. The functional and performance tests are done
	using IoT (LoRa, Sigfox or NB-IoT) commercial devices.

https://https://lora-alliance.org/

https://www.3gpp.org/news-events/1785-nb_iot_ complete

available at https://wirelesslibrary.labs.b-com.com/ FIRfilterdesigner/

taps that need to be stored within the block. In case of LoRa, the considered number of parallel channels is set to 8, hence the lower memory needs. Moreover, it must be emphasized that a small (cost-effective) FPGA has been used for our test, allowing for the embedding of at most two technologies on the same board. However, the three technologies could be easily embedded on a larger board, allowing for a true multi-standard DFE.

The integration of the DFE within the complete network frastructure has been initiated, and preliminary measurements have been conducted. Fig. 3 shows some of the sensitivity measurements results obtained for the LoRa case. The curves represent the block error rate as a function of the input signal level, for different spreading factor configurations (SF).

Conclusion

In this paper, we present a flexible multi-standard digital DFE for the main LPWA technologies: LoRa, Sigfox and NB-IoT. The solution is based on three hardware subsystems that can be configured and tuned according to the targeted technology: rate adaptation, filtering and channel demultiplexing. In the future, we plan to assess and further improve the full RF receiver performance for the three technologies, and particularly immunity to adjacent signals and blockers. We are also working on the integration of hardware detection blocks dedicated to the above mentioned technologies with encouraging preliminary results.


Christophe Delacourt

Patrick Savelli

Vincent Savaux

A Cloud RAN Architecture for LoRa

This paper deals with a cloud radio access network (CRAN) architecture for the LoRa system. In the suggested design, the gateway embeds a limited remote radio head (RRH), including the analog radio-frequency (RF) analog part, the digital-to-analog and analog-to-digital conversion, and a digital front-end (DFE). The other LoRa network functions, including the physical (PHY) layer, the LoRaWAN medium access control (MAC) layer, and the application and customer servers are implemented as cloud resources. The presented approach leads to a flexible RAN that is robust to the variations of capacity needs. Furthermore, it allows us to test very specific LoRa features, such as the detection or demodulation, while bypassing the other ones including the hardware RRH. The methodology and tools we used to deploy a LoRa cloud RAN are detailed, and results concerning the performance indicator (CPU load, memory consumption) are provided as well.

Introduction

LoRa is a low power wide area network (LPWAN) operated in unlicensed frequency bands for Internet of things (IoT) applications. This cost effective, long range and energy efficient technology is being actively deployed globally by network operators, both for public or private networks coverage. The typical architecture of a LoRa network relies on gateways connected through an Internet protocal (IP) based backhaul interface to the centralized LoRa network servers. The gateway is the network equipment that exchanges data with the IoT devices through the air interface and relays messages between these end-devices and the network server. The gateway typically implements the RF and PHY layer functions, including the encoding of the transmit signal, the radio frame generation and up-conversion in the uplink direction, and the received signal down-conversion, demodulation and channel decoding in the downlink direction. Note that the PHY layer is a proprietary solution designed by Semtech. Basics of the modulation (chirp spread spectrum), channel coding, PHY header and payload construction are described in [START_REF] Augustin | A Study of LoRa: Long Range & Low Power Networks for the Internet of Things[END_REF]. The gateway is connected to the network server (NS), that implements the Lo-RaWAN®MAC layer [START_REF]Lorawan™specification[END_REF], and to the application server (AS) that is involved in some of the LoRaWAN security procedures (e.g. management of join requests, encryption of application payload). The AS provides an interface to the end user to collect and send data to the devices, typically through a web application. The LoRaWAN protocol is an open specification maintained by the LoRa Alliance®. The architecture proposed in this contribution consists in offloading some of the processing functions of the gateway in the cloud infrastructure. In this approach, the network is deployed with limited hardware equipment installed on the field, to implement the remote radio head (RRH) functions with a software defined radio (SDR) approach. Basically, this hardware contains the antennas, the RF analog parts, the high speed analog-to-digital and digital-to-analog converters, and a field programmable gate array (FPGA) reconfigurable hardware to perform functions such as rate adaptation, filtering, channels multiplexing (in the transmit side) and de-multiplexing (in the receive side), and the generation of the modulation signal. An overview of the possible RRH solutions is provided in [START_REF] Desnos | Flexible multi-standard digital front-end for lpwa technologies[END_REF]. The other gateway functions, such as the PHY layer baseband processing (demodulation and channel encoding and decoding), as well as the NS and AS functions, are implemented on cloud resources managed through orchestration services, using virtualization technology. This functional split between the centralized baseband functions, also named the baseband unit (BBU), and the RRH installed on the top towers close to the antenna is an architecture evolution referred to as C-RAN. This approach offers easier deployment and scaling capability through dynamic shared resources allocation, where the software network functions are instantiated on the fly in cloud infrastructures according to the network load. It allows for easier network upgrades, enhancements, testing, monitoring, and maintenance. In additionn the centralization of the PHY processing natively enables advanced signal processing techniques such as the joint demodulation of frames received by multiple gateways The rest of the paper is organized as follows: Section 2 is dedicated to the description of our proposed architecture. Sections 3 and 4 present the methodology and the results, respectively, and Section 5 concludes this paper.

2 Cloud RAN Architecture

Overview

The cloud-RAN IoT architecture considered in this paper is depicted in Fig. 1, where the receiver blocks only are il- lustrated (simpler blocks are used in the transmitter). It is composed of two main components: hardware RRH located on the antenna towers, and the other functions of the architecture that are offloaded in the cloud infrastructure. The hardware RRH is, in turn, is composed of analog and digital functions: the analog radio-frequency (RF) blocks of amplification, down conversion, filtering, and analog-to-digital conversion (ADC) in the receiver. A reconfigurable FPGA embeds the DFE functions of rate adaptation [START_REF] Zeineddine | Efficient arbitrary sample rate conversion for multi-standard digital front-ends[END_REF], low-pass filtering [START_REF] Paquelet | On the symmetry of FIR filter with linear phase[END_REF], and the channel demultiplexing. Note that the RRH also contains the transmitter functions.

As illustrated in Fig. 1, the LoRa protocol stack is fully implemented in software in the cloud. This includes the PHY layer and packet forwarder, usually embedded in the antenna tower hardware in state-of-the-art solutions. The PHY layer contains the Rx digital processing chain: synchronization, symbol demodulation and de-mapping, channel decoding. It also implements the channel encoding for the Tx path. The processing is distributed in "dockerized" micro-services 1 , interfaced through remote procedure callbased interfaces, which allows for flexibility and dynamic scaling according to the network processing load. The PHY layer is interfaced with the LoRa network server through the packet forwarder interface. The network server, application server and customer server are based on open source software. The software functions are orchestrated with Kuber-1 https://www.docker.com/ netes2 , which offers deployment automation facilities, and efficiently handles the maintainability and scalability of the platform.

Detailed Functions

The PHY layer software architecture relies on a micro service approach where each different functional component of the Tx and Rx chains is a separated entity that provides a specific service with a dedicated interface. Each of the different micro services run on a separate docker container, to provide an abstraction that packages code and dependencies together. Multiple containers can run on the same machine and share the OS kernel with other containers, each running as isolated processes in user space. Docker containers are self contained software environments for development that include all elements needed to run: code, runtime, system tools and libraries. The use of containers ensures portability between different platforms and clouds, and enable easier scaling of the application.

The different micro services that compose the suggested LoRa PHY software are summarized as follows, and illustrated illustrated in Fig. 2: • Driver: the driver is in charge of configuring the RF and the FPGA components (DFE for the Rx path, and modulator for Tx), and encapsulates IQ data stream from each channel in network packets.

• Replay server: The replay server is intended to buffer the incoming IQ packets stream for retransmission when the synchronization or demodulation service requests it. Its main goal is to overcome the network latency introduced by the microservices. In any case the introduced latency is negligible (of order 10 -2 s) compared to the duration of the LoRa transmission window.

• Detection: this component implements the received LoRa frame detection for all spreading factors. The detection function is a continuous process as required by LoRa standard and therefore it generates a constant CPU load per channel.

• Triage: this service coordinates all the reception chain. Triggered by a detection event, it manages synchronization and demodulation processes.

• Synchronization: the synchronization service starts by requesting retransmission of the packets to the replay server, then it processes the frame time and frequency synchronization.

• Demodulation: the demodulation service is responsible for demodulating and decoding the LoRa PHY frame. This frame includes the PHY header and payload parts. The service starts by requesting retransmission of the packets to the replay server, applies synchronization corrections, and then decodes the header part to determine the parameters required for the payload decoding (channel coding rate, payload length, presence or not of a cyclic redundancy check (CRC) in the payload). It then performs the channel decoding of the PHY payload part, prior to the decoded frame transmission to the MAC layer.

It must be precised that each LoRa gateway channel has exactly one replay server, detection and triage services. However, as depicted in 2, synchronization and demodulation services are stateless and can be scaled/mutualized among multiple gateways, depending on the load of the system. This feature allows for a very flexible CRAN implementation.

In addition to the previous full-CRAN architecture, we also suggest an alternative implementation of LoRa PHY layer hardware/software split, where less blocks are deported in the cloud. This architecture, depicted in Fig. 3, is especially relevant in use cases where the bandwidth between the RRH and the cloud services is limited. Thus, we keep the same hardware DFE, and we propose an ARM-based cost-effective board embedding parts of the micro-services. In this configuration, only synchronization and demodulation are performed in the cloud, and thus the bandwidth between the gateway and the cloud is only required when a message is detected. The peak load in this case is 1MB/s per channel when a frame is detected, which is the average load in the first functional plit proposed on the Fig. 2.

Methodology & Tools

In order to ease development and tests of the microservices, we developed a hardware simulator called LiveRF. Given an input scenario described in json format, LiveRF is able to generate, on the fly, the IQ data stream for each LoRa channel. This microservice can replace the driver mi- croservice in order to perform receiver performance simulations, including signal generation, addition of noise and RF effects and impairments such as frequency error and drift, phase rotation, AGC, quadrature error, etc. LiveRF is used to perform exhaustive message parameters testing, load tests, and complex test such as frames collision, which is difficult to reproduce in real enviromnent conditions. It is also used for continuous integration, as it allows to perform functional non regression testing, and sensitivity performance evaluation of the RX chain. Moreover, the tool Monitor has been developed, a specialized monitoring microservice which consolidates all the information and events provided by the RX chain. Based on Monitor, we are able to develop a number of tools, most notably a dedicated PHY monitoring web application showed in figure 4. This application is a viewer of the I and Q We have also developed a frame capture tool that has proven to be very useful for debugging purposes. This tool captures I and Q signals of frames that can be replayed using LiveRF for further analysis. This process can be triggered according to a set a filtering events. For example, one may capture only the frames that have been unsuccessfully decoded, for a given SF, or on a given frequency channel. Finally, Monitor is also used within our automated GW sensitivity measurement framework, used both in real condition with a controlled LoRa RF generator, or in simulation mode with LiveRF. In addition, during the whole development phase, we have performed regular tests with real LoRa devices and run numerous sensitivity performance campaigns to validate the gateway and calibrate the simulator.

Results

Based on the proposed architecture and methodology, we have developed an efficient, robust and perfectly interoperable platform with existing LoRa devices. The network architecture has been deployed and tested using commercial devices in a real environment. depending on the network load. Furthermore, there is still a lot of room for fine level software optimizations. Anyway, we have performed preliminary load tests on a Kubernetes cluster composed of 5 virtual machines with 4 cores and 8GB Ram, for a total of 20 cores and 40GB RAM. The underlaying hardware is an Intel Xeon Gold 6138 CPU @ 2.00GHz. Results are given in terms of CPU load (100% corresponds to 1 full loaded core) and used memory in MB. They are measured after 10 minutes of a constant load of 6, 14 and 24 LoRa frame/s as depicted in Table 1. It must be noticed that 24 frames/s greatly overestimates a usual LoRa devices deployment. We conclude that the proposed architecture is able to support practical use cases on small size clusters.

Conclusion & Future work

In this paper we have presented a new cloud RAN architecture for the LoRa system. This architecture is flexible and has been derived into two configurations: a full-RAN solution and a cost-effective ARM based board to limit the required bandwidth between the RRH and the cloud.

For the LoRa standard, we have planned to work on joint demodulation of messages between multiple gateways to further improve reception performance. Finally, our goal is to extend this flexible and reconfigurable CRAN solution to support multiple IoT protocols.

Figure 1 .

 1 Figure 1. Cloud RAN Architecture: the gateways (GWs) only embed the analog RF and the DFE, the whole LoRa protocol stack is processed in external servers.

Figure 2 .

 2 Figure 2. A possible implementation of LoRa PHY Layer hardware/software split.

Figure 3 .

 3 Figure 3. Alternative implementation of LoRa PHY Layer hardware/software split: driver, detection, triage, and replay server are deported in a low-cost ARM CPU based card.

Figure 4 .

 4 Figure 4. LoRa PHY monitoring Web Application

Table 1 .

 1 Preliminary load test results

	The inherent flexibility of the
	architecture enables quick prototyping of any kind of new
	processing algorithms and tools. Addition of new features
	and upgrades are also greatly simplified. Intensive auto-
	mated testing significantly reduces the development cycle,
	thanks to continuous deployment and Kubernetes.
	The major drawback of the approach is the network
	overhead mechanically introduced by the microservices.
	Futhermore, comparatively to a monolithic implementa-
	tion, a lot of memory optimization are prevented, leading to
	a slightly slower execution time. However, this limitation is
	largely mitigated by the excellent scaling capabilities of the
	architecture in a cluster environment. The frame detector
	has a fixed CPU cost and memory usage (continuous detec-
	tion), whereas synchronization and demodulation services
	load depends on the number of detected messages. Those
	two services are stateless and can be scaled or mutualized

https://kubernetes.io/