Zhaobin Li

Ganbo Deng

Patrick Queutey

Benjamin Bouscasse
email: benjamin.bousscasse@ec-nantes.fr

Guillaume Ducrozet

Lionel Gentaz

David Le Touzé

Pierre Ferrant

Comparison of wave modeling methods in CFD solvers for ocean engineering applications

Keywords: Wave modeling in CFD, Relaxation Zone, Internal Wave Generator, Potential/viscous flow coupling, SWENSE

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

Introduction

Reproducing the open sea environment in a finite domain is a common challenge faced by both physical and Numerical Wave Tanks (NWT) for ocean engineering applications. The key technique is to generate incident waves and to prevent the reflection at the domain boundaries at the same time.

A physical experimental wave tank usually generates waves with wavemakers and uses an absorption beach to prevent wave reflection on the other end.

In NWTs, reproducing a wavemaker using moving mesh and a wave-absorption beach is computationally expensive and thus is often replaced by numerical wave generation and absorption models. The simplest method is to use the incident wave velocity and the free surface elevation as boundary conditions. This kind of boundary allows incident waves to travel in and out but causes wave reflections when the real waves are different from the incident waves, e.g., diffracted waves by structures. Without proper treatments, large domains should be used to secure a time window long enough for data analysis before the reflected waves travel back (same as in physical wave tanks, but with a high computational cost).

Efforts in two major directions to overcome this difficulty can be found in the literature.

The first direction focuses on developing generationabsorption boundary conditions, aiming to use the boundary to generate incident waves and absorb the diffracted waves at the same time, for example, the weakly reflecting boundary condition proposed by [START_REF] Petit | Numerical simulation and validation of plunging breakers using a 2D Navier-Stokes model[END_REF] and the active absorption boundary condition proposed by [START_REF] Higuera | Simulating coastal engineering processes with OpenFOAM R[END_REF]. However, these methods use linear superposition to combine the incident and the diffracted waves and are thus limited to small diffraction cases.

Methods in the second direction treat the wave generation and the absorption problems separately.

The Internal Wave Generator (IWG) method generates incident waves in a specific zone inside the computational domain and absorbs all the outgoing waves with sponge-layers on boundaries. The IWG method is originally proposed by [START_REF] Larsen | Open boundaries in short wave simulations-a new approach[END_REF] in a Boussinesq equation model, and has been extended to Navier-Stokes models by adding source terms either in the mass conservation equation [START_REF] Lin | Internal wave-maker for Navier-Stokes equations models[END_REF] or in the momentum equation [START_REF] Choi | Numerical simulations using momentum source wave-maker applied to RANS equation model[END_REF] to generate unidirectional regular or irregular waves. [START_REF] Chen | Generation of 3D water waves using mass source wavemaker applied to Navier-Stokes model[END_REF] and [START_REF] Ha | Generation of 3d regular and irregular waves using navier-stokes equations model with an internal wave maker[END_REF] have enabled multi-directional wave generation by using two perpendicular IWG zones. The IWG is transparent to the flow, allowing any waves passing through the wavemakers without interference, so this is particularly good for cases with significant wave reflection, such as simulating wave shoaling and breaking on the beach [START_REF] Lara | Rans modelling applied to random wave interaction with submerged permeable structures[END_REF] and waves interacting with structures, such as, wave energy converters [START_REF] Schmitt | On the use of openfoam to model oscillating wave surge converters[END_REF][START_REF] López | Optimisation of turbineinduced damping for an owc wave energy converter using a rans-vof numerical model[END_REF], floating offshore wind turbines [START_REF] Calderer | Fluid-structure interaction simulation of floating structures interacting with complex, large-scale ocean waves and atmospheric turbulence with application to floating offshore wind turbines[END_REF], etc.

Another widely used method is the Relaxation Zone (RZ) method. Differently from the IWG method where waves are generated anywhere in a zone inside the computational domain, the RZ method blends the numerical solution obtained by the solver to a reference solution in a zone near the boundaries of the domain. The velocity and the free surface position computed are "relaxed" to an incident wave solution, this generates waves out of the RZs and also avoid outgoing waves being reflected. This method was initially used by [START_REF] Mayer | A fractional step method for unsteady free-surface flow with applications to non-linear wave dynamics[END_REF] to only absorb outgoing waves by setting the incident solution to zero. After being extended for wave generation and implemented in OpenFOAM by [START_REF] Jacobsen | A Wave Generation Toolbox for the Open-Source CFD Library: OpenFoam R[END_REF], the RZ method is now widely used by the ocean and coastal engineering community. It has been applied to simulate regular and irregular wave propagation [START_REF] Choi | Generation of regular and irregular waves in Navier-Stokes CFD solvers by matching with the nonlinear potential wave solution at the boundaries[END_REF], wave breaking on a steep beach [START_REF] Chella | Breaking characteristics and geometric properties of spilling breakers over slopes[END_REF], and waves interacting structures, such as, substructures of offshore wind turbines (Paulsen et al., 2014a,b), wave energy converters [START_REF] Kamath | Numerical investigations of the hydrodynamics of an oscillating water column device[END_REF][START_REF] Palm | Coupled mooring analysis for floating wave energy converters using cfd: Formulation and validation[END_REF], coastal bridges [START_REF] Seiffert | Experiments and computations of solitary-wave forces on a coastal-bridge deck. part i: Flat plate[END_REF][START_REF] Hayatdavoodi | Experiments and computations of solitary-wave forces on a coastal-bridge deck. part ii: Deck with girders[END_REF], and ships [START_REF] Seo | Effect of wave periods on added resistance and motions of a ship in head sea simulations[END_REF].

Moreover, coupling viscous flow solver with potential wave theory is also used to model waves in NWTs, as in the Spectral Wave Explicit Navier-Stokes Equations (SWENSE) [START_REF] Ferrant | A potential/RANSE approach for regular water wave diffraction about 2-D structures[END_REF]Vukčević et al., 2016a;[START_REF] Li | Progress in coupling potential wave models and two-phase solvers with the SWENSE methodology[END_REF][START_REF] Li | Spectral Wave Explicit Navier-Stokes Equations for wave-structure interactions using twophase Computational Fluid Dynamics solvers[END_REF] method. The SWENSE method solves the incident wave propagation problem with fully non-linear potential wave models, and uses the CFD solver only to solve the complementary field. Since the incident waves are no more treated by the CFD solver, one needs only to prevent the complementary field (the diffraction, radiation, etc) from reflecting in the far-field. This can be easily achieved since the complementary waves decrease naturally as they travel away. An absorption approach may also be used to damp out the complementary waves in the far-field. This method has been successfully applied to simulate wave-structure interactions, such as ships in waves [START_REF] Luquet | Viscous flow simulation past a ship in waves using the swense approach[END_REF][START_REF] Reliquet | Simulation of wave-body interaction using a single-phase level set function in the swense method[END_REF] and vertical cylinder in waves (Vukčević et al., 2016b;[START_REF] Li | Calculation of high-order wave loads on a vertical circular cylinder using the swense method[END_REF]. [START_REF] Miquel | Analysis of different methods for wave generation and absorption in a cfdbased numerical wave tank[END_REF] compared the active wave absorption method (from the first category) and the RZ method (from the second category) and concluded that the latter is more effective.

The present work aims to compare the wave modeling methods from the second category, namely, the IWG, the RZ, and the SWENSE methods, for ocean engineering applications. The three methods are implemented in a similar numerical framework: using the Finite Volume Method and simulating the twophase air water flow by the Volume of Fluid (VOF) technique [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF]. The IWG method is implemented in ISIS-CFD1 , a CFD code developed at Ecole Centrale de Nantes and distributed commercially by NUMECA International as a part of the FINE/Marine computing suite2 ; the RZ and the SWENSE methods are implemented in foamStar, developed by Ecole Centrale Nantes and Bureau Veritas, based on Open-FOAM3 , a widely used open source CFD software package. This comparison is achieved by studying a real ocean engineering application case, concerning a fixed Catenary Anchor Leg Mooring (CALM) buoy exposed to regular incident waves in open sea. The experiment carried out in the ocean engineering tank at Ecole Centrale de Nantes is also provided for the comparison. The wave generation and absorption performance and the computational efficiency of the three methods are mainly concerned. A problem related to the turbulence modeling in two-phase flow is also remarked.

The paper is structured as follows. After this introduction, the theory part describes the governing equations of NWTs and the principles of the three wave generation and absorption methods. Sect. 3 provides information of the model test in the experimental wave tank. Sect. 4 includes the numerical study of the same case, the data analyses, and the comparison. The conclusion is given in the end.

Theory

This section provides the theory about establishing a NWT using two-phase CFD solvers. In Sect. 2.1, the basic governing equations of two-phase flow is given. In Sect. 2.2, the theory of the IWG method, the RZ method, and the SWENSE method is described.

Governing equations

In classical ocean engineering application, the fluid is assumed as incompressible and the surface tension is ignored. In this condition, the two-phase air-water flow in a NWT is modeled with the Reynolds Averaged Navier-Stokes Equations (RANSE), as follows:

∇ • u = 0 (1) ∂ (ρu) ∂t + ∇ • (ρu ⊗ u) = -∇p + ρg + ∇ • µ eff ∇u + ∇u T (2)
where u is the velocity vector field, ρ is the density of the fluid, p is the pressure field, g is the gravitational acceleration vector, and µ eff is the effective dynamic viscosity defined as µ eff = µ + ρν t , with µ the fluid dynamic viscosity and ν t the turbulent eddy viscosity. k -ω SST model is used to determine the turbulence eddy viscosity [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF]. The fluid properties ρ and µ remain to be determined with interface capturing techniques, such as the VOF method [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF] in the present implementations.

The VOF method defines a scalar field α representing the volume fraction of water in a cell. This field is equal to 1 when the cell is fully occupied by water, and is equal to 0 when the cell is full of air. Cells on the interface have 0 < α < 1. The fluid properties are defined with an interpolation between water and air using α as weight.

ρ = αρ water + (1 -α)ρ air (3) µ = αµ water + (1 -α)µ air (4)
The VOF field is advected by the velocity field u, as follows,

∂α ∂t + ∇ • (uα) = 0 (5)
These equations are the common and fundamental part of two-phase VOF solvers and are implemented in ISIS-CFD.

In foamStar, the momentum equation (Eqn. 2) is expressed in a slightly different form as follows,

∂ (ρu) ∂t +∇•(ρu ⊗ u) = -∇p d -g•X∇ρ+∇• µ eff ∇u + ∇u T (6)
where p d is the dynamic pressure defined as p d = p + ρg • X, and X = (x, y, z) is the vector field of the cell center coordinates. Eqns. ( 2) and ( 6) are fully equivalent.

The VOF transport equation in foamStar also contains a small modification, as follows

∂α ∂t + ∇ • (uα) + ∇ • (u r α(1 -α)) = 0 (7)
Compare to Eqn. ( 5), the extra term represents an artificial compression to keep the interface sharp [START_REF] Rusche | Computational fluid dynamics of dispersed two-phase flows at high phase fractions[END_REF].

Wave generation and absorption methods 2.2.1. IWG

The IWG method produces waves in a specific region inside the computational domain by adding a momentum source function in the governing equations. The implementation in ISIS-CFD is inspired by the method of [START_REF] Choi | Numerical simulations using momentum source wave-maker applied to RANS equation model[END_REF]. A momentum source term is added in Eqn. (2). The modified momentum equation reads,

∂ (ρu) ∂t + ∇ • (ρu ⊗ u) = -∇p + ρg + ∇ • µ eff ∇u + ∇u T + rρa I (8)
where rρa I is the extra momentum source term calculated using incident wave kinematics. r is a scalar field that is equal to one inside the wave generation zone and to zero elsewhere, as shown in Fig. 1(a).

At the outlet boundary, a damping zone is applied to damp the waves to zero with numerical wave damping techniques.

RZ

The RZ approach [START_REF] Jacobsen | A Wave Generation Toolbox for the Open-Source CFD Library: OpenFoam R[END_REF] is used together with wave generating and absorbing boundary conditions. NWTs using such boundary conditions alone suffer from the reflection problem when the simulated waves and the imposed incident waves are different at the boundary. To avoid such differences, RZs are added (see Fig. 1(b)) to gradually blend the computed value to the target value and force them to match at the boundary. The relaxed velocity and VOF fields in these regions are defined as a linear combination of the CFD solution and the target incident wave value, as follows:

α relax = ωα I + (1 -ω)α (9) u relax = ωu I + (1 -ω)u (10)
where the variables without subscripts denote the solution of CFD solver and the subscript I denotes the incident wave values. In foamStar, the incident waves are defined by an external wave solver based on fully non-linear potential theory [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF][START_REF] Ducrozet | CN-Stream: Open-source library for nonlinear regular waves using stream function theory[END_REF][START_REF] Ducrozet | A modified highorder spectral method for wavemaker modeling in a numerical wave tank[END_REF][START_REF] Ducrozet | HOS-ocean: Opensource solver for nonlinear waves in open ocean based on High-Order Spectral method[END_REF]. ω is a space-dependent weight function chosen here as:

ω = 1 -3d 2 + 2d 3 ( 11
)
where d is the normalized distance in the RZ defined as the distance to the boundary divided by the length of the RZ [START_REF] Monroy | Développements et validation de l'outil CFD OpenFOAM pour le calcul de tenue à la mer[END_REF]. This blending procedure is applied at the end of each time step after a converged solution has been obtained. The relaxed values are then used in the further simulation. The pressure field is not relaxed, since it is obtained directly from the Poisson equation in the incompressible flow case.

This technique is used to generate incident waves at the inlet of the CFD domain, and can also prevent wave reflections at the boundary by letting the incident part of the waves go out of the domain and absorbing the other parts with the RZ. A thorough study on the optimal parameter setup can be found in the work of [START_REF] Choi | Generation of regular and irregular waves in Navier-Stokes CFD solvers by matching with the nonlinear potential wave solution at the boundaries[END_REF].

SWENSE

The SWENSE method [START_REF] Ferrant | A potential/RANSE approach for regular water wave diffraction about 2-D structures[END_REF] is a specific method designed to solve the wave-structure interaction problem using a potential/viscous flow coupling approach. It treats the wave-structure interaction problem by decomposing the total fields into: i) the incident part, and ii) the complementary part. The incident part describes the incident wave propagation in the computational domain without the structure; the complementary part represents the diffracted, radiated waves, and viscous effects. This decomposition is illustrated by Fig. 2. The decomposition is applied in the entire computational domain, no specific wave generation or absorption zone is defined as in the IWG and the RZ.

The two decomposed parts are treated separately.

The incident wave part is addressed by non-linear potential flow theory, under the inviscid and irrotational flow assumption and solved by spectral wave models. Regular incident waves are solved by the stream function wave theory [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF][START_REF] Ducrozet | CN-Stream: Open-source library for nonlinear regular waves using stream function theory[END_REF], while irregular incident waves are solved with High-Order Spectral (HOS) method [START_REF] Ducrozet | A modified highorder spectral method for wavemaker modeling in a numerical wave tank[END_REF][START_REF] Ducrozet | HOS-ocean: Opensource solver for nonlinear waves in open ocean based on High-Order Spectral method[END_REF].

The complementary part is treated by a viscous CFD solver. The governing equations of the complementary field are derived using the following idea: a variable χ (velocity, pressure) in Navier-Stokes equations (Eqns. 1 and 2) is decomposed into an incident variable χ I and a complementary variable χ C , as follows,

χ = χ I + χ C ( 12
)
The incident variable χ I represents the inviscid flow and is described by Euler equations. The governing equations of χ C are mathematically derived by substituting Eqn. ( 12) into the Navier-Stokes equations (Eqns. 1 and 2) and simplifying them by Euler equations.

The original SWENSE method proposed by [START_REF] Ferrant | A potential/RANSE approach for regular water wave diffraction about 2-D structures[END_REF] was made for single-phase solvers. It has been extended for two-phase CFD solvers in two recent works using different decomposition strategies. In the first work (Vukčević et al., 2016a,b), the SWENSE decomposition is applied only to the velocity field, leaving the entire pressure field to be solved by the viscous solver. In contrast, in the second work [START_REF] Li | Progress in coupling potential wave models and two-phase solvers with the SWENSE methodology[END_REF][START_REF] Li | Spectral Wave Explicit Navier-Stokes Equations for wave-structure interactions using twophase Computational Fluid Dynamics solvers[END_REF], the SWENSE decomposition is applied to both the velocity and pressure fields. This more complete incident and complementary decomposition leads to less solution effort of the viscous solver. For this reason, the method in the second work is used in the present comparison.

The governing equations of the complementary fields read:

∇.u C = 0 (13) ∂u C ∂t + u C .∇u C + u C .∇u I + u I .∇u C = - ∇p C ρ - p I ρ I ∇ρ ρ + ∇ • µ eff ∇u C + ∇u T C ρ ( 14
)
where ρ I = ρ water is the water density, and p I is the incident wave pressure provided by nonlinear potential theory. All the incident fields are extended to the viscous solver's computational domain. After solving the complementary field, the total field is reconstructed to obtain the final solution.

The merit of such decomposition is on the mesh requirement regarding the incident wave propagation. This one can be loosened since the incident wave information is explicitly known. A high mesh density is only necessary near the structure to solve the complementary field with a high level of accuracy.

The boundary condition of the complementary field at farfield boundary is defined to be zero. To ensure the complementary field vanish at far-field and to prevent its reflection, the outgoing complementary waves are either absorbed with coarse mesh or with RZs in the far-field zone.

Experimental study

The experimental study deals with a fixed CALM buoy [START_REF] Rousset | Model tests for Principia R&D[END_REF] in regular incident waves. The complexity of the geometry and of the flow represents a realistic ocean engineering application case, thus this study is suitable for comparing the different NWT models.

Experimental setup

The experiment has been conducted in the ocean engineering basin of ECN, which is rectangular and is 50m long, 30m wide, and 5m deep. One side of the tank is equipped with a segmented wavemaker, composed of 48 hinged flaps distributed over the width of the basin. Each flap is controlled separately. The wave generator system is equipped with an active wave absorption control feature. The opposite end of the basin is equipped with a passive wave breaking beach to prevent wave reflections. The side walls are reflective.

The geometry of the CALM buoy contains a truncated cylinder and a thin skirt near the bottom to provide additional heave damping force through vortex shedding, as shown in Fig. 3 and Tab. 1. In the experiment, the buoy is captured by a tripod and fixed at its design draft position. The center of the buoy locates on the basin's center-line and is 17.19m away from the wavemaker.

Uni-directional regular incident waves are used in this case. The wave parameters are provided in Tab. 2. Wave forces and moments on the buoy are measured by a 5D-dynamometer (Mz is not measured). The free surface elevations are measured at four positions with resistance wave gauges. The data is acquired at 100Hz.

Three probes are located near the structure (see Fig. 4). Probe No.4 is located at 8m aside from the buoy and 17.

Sliding window Fourier transform

The amplitudes of the zeroth (the mean value), the first and the second harmonic components of the measurement are obtained by the sliding window Fourier transform.

The sliding windows have a size equal to one wave period (T = 1.8s), the time increment between two windows is ∆t = 0.01T . After such a transform, the amplitudes are recorded in function of time. Their time-averaged values and the standard deviations can be calculated in a chosen analysis zone.

Non-dimensionalization

The data is non-dimensionalized to take into account the probable discrepancies between the target and the measured incident waves. Such discrepancies are related to the use of a linear transfer function for the control of the wavemaker.

Figure 5 shows the free surface elevation measured by wave probe 4 (η 4 ), together with its first harmonic amplitude (η (1) 4 ) extracted with the sliding window Fourier transform.

Ideally, its amplitude should remain constant and be equal to that of the target incident wave, if i) the wave generation and absorption are perfect and ii) this position is not affected by the diffracted waves. However, the measured wave amplitude contains several variations (see Fig. 5(b)), which indicates that the above assumptions are not fully satisfied. Such variations are interpreted as:

(A) At t ≈ 45s, the arrival of the wave front: where small amplitudes waves are followed by steeper waves before the amplitude becomes nearly constant.

(B) At t ≈ 55s, the arrival of diffracted waves: at first the diffracted waves from the buoy reach the probe; several periods later the waves get reflected at the side wall travels back and reach the probe again;

(C) At t ≈ 85s, the arrival of reflected waves from the absorbing beach.

As a result, only the interval after the peak A and before the reflection B corresponds to the pure incident waves generated by the wavemaker at this probe. The measured wave height (crest to trough) in this interval is:

H m ± σ H m = 0.152m ± 0.0015m
where H m is the mean value and σ H m is the standard deviation. It is worth noting that the measured mean wave height H m is slightly smaller than 0.16m, the target value.

In wave tank experiments, this difficulty in generating the exact target incident waves is overcome by converting the raw experimental data to non-dimensional values with a correction taking the measured incident wave information into account.

The non-dimensional wave forces are defined as

F nondim = F A m kρgV ( 15
)
where F is the measured force, A m is the measured amplitude, k is the wave number, and V is the immersed volume. The measured wave amplitude A m is defined as

A m = 1 2 H m (16)
The non-dimensional wave elevation is defined as

η nondim = η A m (17)
with η the measured free surface elevation.

Wave forces and free surface elevations

The first harmonic amplitudes of horizontal force and vertical forces (dimensionless values) are shown in Fig. 6.

Similarly as for the incident wave measurement, the three variations are also observed in the force amplitudes:

(A) At t ≈ 45s, the arrival of the wave front;

(B) At t ≈ 65s, the arrival of reflected diffracted waves from side wall;

(C) At t ≈ 85s, the arrival of reflected waves from the beach. It is worth noting that the variation (B) happens about 10 seconds later for the force than that for the wave elevation at probe 4, since only the reflection caused by the side walls has an influence on the force and these waves need more time to propagate back to the center of the basin.

To obtain final results describing the true response to the incident waves, only data in the effective time window (after the peak A and before the variation B, where the amplitudes are almost constant) are used.

The experimental results are summarized in Tab 3. In the table, the superscript i represents the ith harmonic component. The results of first harmonics are given in the first place, followed by the second harmonics. The zeroth harmonic amplitude, representing the drift force, is also provided in the end. For each component, the mean value (x) corresponds to the time-averaged value during the zone of analysis and the standard deviations (σ x ) include already the error on the incident waves in the non-dimensionalization. The relative standard deviation (σ x /x) represents the difficulty to obtain a precise result.

In this experiment, all the first harmonic amplitudes have relative standard deviations about 1%, showing a good data quality. However, the second harmonics are more difficult to obtain and contain larger relative standard deviations, since the mean values are much smaller. Only the horizontal force (F (2) x ) has a fairly small standard deviation of about 5%. The zeroth harmonic amplitude of the horizontal force (F (0) x ) has also a relative standard deviation of about 5%.

Numerical Study

This section reproduces the wave tank experiment using the three NWTs to compares the three wave modeling methods. The following features are considered:

• the accuracy;

• the computational efficiency;

• the wave generation and absorption ability.

The section begins by introducing the CFD solvers implementing the three wave modeling methods, and continues with a description of the numerical setup. The numerical results and the discussions are presented in the end.

CFD solvers

The wave generation and absorption methods presented in Sect. 2.2 are implemented in two-phase CFD solvers. The IWG method is implemented with ISIS-CFD, while the RZ method and the SWENSE method are implemented in foamStar. Both ISIS-CFD and foamStar are based on the second-order accurate Finite Volume method. Unstructured polyhedral mesh with colocated variable arrangement is used. The governing equations are solved with the SIMPLE algorithm [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF]. The free surface is treated with the VOF method [START_REF] Hirt | Volume of fluid (VOF) method for the dynamics of free boundaries[END_REF].

foamStar is based on interDyMFoam, the native OpenFOAM solver for incompressible two-phase flow. The solver's main structure is maintained; only slight modifications are made to implement the wave generation and absorption methods. Details can be found in [START_REF] Monroy | Développements et validation de l'outil CFD OpenFOAM pour le calcul de tenue à la mer[END_REF] for the RZ implementation and in [START_REF] Li | Progress in coupling potential wave models and two-phase solvers with the SWENSE methodology[END_REF][START_REF] Li | Spectral Wave Explicit Navier-Stokes Equations for wave-structure interactions using twophase Computational Fluid Dynamics solvers[END_REF] for the SWENSE method.

Numerical setup

Two types of mesh are used in this test case. One is rectangular and the other is cylindrical, as shown in Figs. 7 and8.

The rectangular configuration is normally used by most numerical wave tanks because its uniform and regular mesh is usually more accurate for the incident wave propagation.

The cylindrical configuration has fine cells only near the structure and coarse mesh in the far-field. This mesh arrangement contains a lower number of cells while keeping the density near the structure the same as in the rectangular configuration. It is consequently often more efficient. However, due to the coarse mesh in the far field, this configuration is not suitable to simulate wave propagation using classical Navier-Stokes models. But this configuration is compatible with the SWENSE method and shows its best performance, since the SWENSE method only need to solve the complementary fields, which is located close to the structure (in the zone with high mesh density).

In both configurations, the origin of the coordinate system is located at the center of the buoy and on the free surface at rest. The x axis points to the wave propagation direction, and z axis points upward. The depth of the computational domain is the same as in the experimental wave tank z min = -5m. Only half of the domain is simulated by using a longitudinal symmetry plane.

The following part only details the spatial discretizations. The temporal discretization is set to have 360 time steps per wave period, which is considered to be fine enough according to a separate study [START_REF] Choi | Generation of regular and irregular waves in Navier-Stokes CFD solvers by matching with the nonlinear potential wave solution at the boundaries[END_REF].

Rectangular configuration

The rectangular configuration in Fig. 7 has Cartesian far-field mesh. Unstructured refined mesh is used locally near the structure to fit the geometry (see Fig. 9). Three meshes with different far-field mesh densities: 20x, 40x, and 80x are used to compare the three models on their spatial discretization requirement. The number 20, 40, and 80 signify the number of cells per wavelength (λ) along the x direction. The far-field meshes are the only difference. The mesh density near the buoy is kept invariant and is assumed to be adequate to capture the flow details.

The mesh information is detailed in Tab. 4, where λ = 5.05m and H = 0.16m are the target incident wave length and height. Please note that the vertical discretization is kept unchanged in the three configurations and is assumed to be enough. The transverse (y) direction is not varying simultaneously with cells in the x direction. Such a configuration is commonly used to reduce the number of cells in unidirectional regular wave simulations. It has a preference wave propagation direction along which the mesh is more refined than the other direction.

Cylindrical configuration

As shown in Fig. 8, the mesh of the cylindrical configuration has fine cells near the center of the domain and coarse cells in the far-field. A local refinement is applied near the structure. The mesh density in this refinement zone is comparable to that of the rectangular configuration.

Using the cylindrical mesh with the SWENSE method assumes that the complementary fields (diffracted waves) need to be solved accurately only in a region close to the structure. This study analyzes the necessary size of this region.

Three meshes with different refinement zone sizes are used, as illustrated by Fig. 10 and referred to as C1, C2 and C3. Their details can be found in Tab. 5. In the table, ∆r near and ∆r far represent the radial size of the smallest and the largest cells respectively, r refine denotes the width of the local refinement zone near the buoy.

Wave generation and absorption

The wave generation and absorption methods used in the rectangular configuration are detailed as follows.

• The IWG method in ISIS-CFD: A wave generation zone of one wavelength long is set next to the inlet boundary (x ∈ (x min , x min + λ)), as shown in Fig. 7(a). The wave damping zone is 1.5 wavelength long, and is next to the outlet boundary (x ∈ (x max -1.5λ, x max )). The side wall is modeled as slip wall, since using damping zones near the side wall is not feasible. This configuration is very similar to the experimental condition. Reflection from the side wall is expected.

• The RZ method in foamStar uses incident wave boundary condition at the inlet, the outlet and the side wall. Three RZs with 1.5 wavelength long are used next to these boundaries (x ∈ (x min , x min + 1.5λ) and x ∈ (x max -1.5λ, x max ) and y ∈ (y max -1.5λ, y max )), as shown in Fig. 7(b). This setup avoids the reflections at any domain boundary to simulate an open sea condition.

• The SWENSE method in foamStar does not need wave generation zone, but uses RZs to prevent the reflections of complementary waves. These zones are set in the same way as in the RZ method to simulate an open sea condition.

For the cylindrical configuration, the SWENSE method in foamStar uses a RZ of 1.5 wavelength long around the far-field boundary to simulate an open sea condition, as shown in Fig. 8.

The stream function wave theory [START_REF] Rienecker | A Fourier approximation method for steady water waves[END_REF][START_REF] Ducrozet | CN-Stream: Open-source library for nonlinear regular waves using stream function theory[END_REF] is used to calculate the target incident waves as shown in Tab. 2.

Results

Same information is collected in the numerical study as in the experiment. The raw data is non-dimensionalized and processed with sliding windows Fourier transform in the same way as in the experiment. For simplicity, only the analysis procedure for the horizontal wave force is detailed here. The complete results are provided in Tab. 6, and further compared in Figs. 15 and16.

It is worth noting that the non-dimensionalization in the numerical study uses the target incident value and not the measured value, so to see the errors in the wave propagation.

IWG in ISIS-CFD

The simulation using the IWG method starts with calm water condition at t = 0. The waves are generated in the wave making zone and propagates in the NWT. Figure 11 plots the horizontal wave force during the first 10 wave periods and its first harmonic amplitude. The three different curves represent the results using different discretizations. From Fig. 11(b), it is observed that the force amplitudes are not constant, but with two obvious variations. Similar to the experiment, the horizontal wave force exhibits a peak after the arrival of the wave front (point A). It is also worth noting that the wave force amplitudes begins to decrease at t ≈ 12s (point B). Although this decrease is similar to the experiment, it is not due to the side wall reflection, since the reflected waves need more time to travel back. This decrease is caused by the damping zone at the outlet, the distance of which is only 1.5 wavelength away from the buoy. The traveling time estimated with the wave group velocity is equal to 10.8s. For this reason, only a short time window after the peak A and before the reflection B can be used for data analysis.

Besides, the wave forces show a converging trend with the mesh refinement. The results using the 20x mesh is about 10% smaller than the two close results with the 40x and 80x meshes.

RZ in foamStar

At the beginning of the simulation, the incident waves are initialized everywhere in the entire computational domain. Such initial condition is necessary to ensure the compatibility of the flow in the pure CFD zone and in the RZs (blending to incident wave).

The results of such initial condition is shown by the peak A in Fig. 12(b)), which corresponds to the impact caused by the sudden appearance of the buoy in the incident waves. The following trough corresponds to the establishment of the diffracted waves. Note that the amplitude curve starts one period later compared to the time history, since the sliding window Fourier transform uses one wave period as the window length.

After point A (t = 1.8s), the force amplitudes remain almost constant for a short interval. It may be interpreted as the response to the incident waves initialized at the beginning of the simulation. Then force amplitudes decrease between point A and point B (t ≈ 7s). This decrease suggests that the waves generated at the inlet boundary become smaller while propagating in the CFD domain due to numerical errors in the wave modeling method. This loss of wave amplitude is more severe with the coarse mesh (20x) than with the fine mesh (80x), but the duration of such a decrease is the same for the three meshes, which is approximately equal to the time to propagate information from the inlet boundary to the structure at the wave group velocity.

After point B, the results become almost constant and correspond to the converged results with each discretization. Compared to the IWG result, the duration of the stable result is longer, so it offers a longer zone to extract the final results.

The results in different meshes have larger discrepancy compared with the IWG's results in ISIS-CFD. initial condition in the SWENSE simulation is the same as the RZ method, i.e., the incident waves exist everywhere. Thus, the wave force at the beginning of the simulation is similar to that of the RZ method near point A.

The point B in Figs. 13(b) and 14(b) are plotted at the same position as in Fig. 12(b), i.e., at t = 7s, since the information propagation time should be identical. However, the decrease of the force amplitudes between A and B is much smaller than that of the RZ method, showing that the loss of wave amplitude caused by the SWENSE method is also smaller than the RZ method.

After point B the simulation results become almost constant, offering a long zone of analysis to extract the final result.

In Fig. 13, the cylindrical mesh simulations show very close results. The force amplitudes increase slightly when enlarging the size of refinement zone near the structure. It suggests that solving the complementary field accurately in different spatial extents has a slight influence. The quantitative study of this problem will be given in Sect. 4.3.4 and the compromise between the accuracy and the computational cost using the SWENSE method will be further addressed in Sect. 4.4.1. Moreover,in Fig. 14, the difference between mesh configurations are very small, revealing that the far-field spatial discretization has only minor influence on the accuracy of wave

Summary of the results

The simulation results are detailed in Tabs. 6 and 7, including the first two harmonic amplitudes of the wave force (horizontal and vertical) and of the wave amplitudes at the three probes close to the buoy. The zeroth harmonic amplitude of the horizontal wave force is also provided. In the table, the mesh names are followed by the number of cells to indicate the computational cost.

These data are obtained with the sliding window Fourier transform, in the same way as in the experiment. Only the numerical results inside the exploitable zone are used for the data processing, to exclude the transitional regime at the beginning and the influence of the wave reflections (in the IWG case). For example, in the interval wave generation case in Fig. 11(b), the result after peak A and before B is used; in the RZ method and the SWENSE method, the results after point B are used (see .

For each result, the mean value is provided and followed by its standard deviation to give an indication of the periodical uncertainty.

Before going into details, the tables show a general good agreement between the experimental data and the numerical results on the mean values, especially when the finest mesh (80x) is used. This agreement indicates that all the three wave modeling methods are able to simulate correctly the waves and the interaction with the structure. With respect to the standard deviation, the results of the IWG method in ISIS-CFD have the largest deviations, which are comparable with that in experimental data. The results of RZ method and the SWENSE method in foamStar have smaller deviations, suggesting a more stable simulation results, which is in agreement with the observation from Figs. 11-14.

Among these quantities, the first harmonics are one or two orders of magnitudes larger than the second harmonics and thus describe the majority of the physical phenomenon. Besides, they contain relatively small standard deviations, showing better data quality. The second harmonics, especially that of the vertical force has much smaller mean values and much larger standard deviations. For this reason, the comparison shall concentrate on the first harmonics forces, which are the main interest from the engineering point of view.

For this reason, the comparison is at first conducted for the first harmonics of the horizontal and vertical forces, in Fig. 15. The first harmonic amplitude of the free surface elevation at probe 1 is also shown. In each figure, the vertical axis denotes the non-dimensional amplitude. The horizontal axis describes the mesh information in ascending order of number of cells. The mean experimental data is plotted with a black solid line, which is encompassed by a red region covering one standard deviation above and below. The mean values of the simulation results are described by the bars. The error bars represent one standard deviation above and below the mean values.

In Fig. 15(a), the first harmonic amplitudes of the horizontal wave forces (F (1) x ) are plotted. All the three methods achieve good results compared to the experimental data when the finest discretization (80x) is used, but their accuracy varies significantly on coarser meshes. With 80x mesh, the difference to the experimental data is smaller than 3% for the foamStar-RZ, while the results of ISIS-CFD-IWG and the foamStar-SWENSE are very close to the experimental data. When meshes become coarser, the behavior of the three wave modeling methods become different. ISIS-CFD-IWG and foamStar-SWENSE can still produce results with errors smaller than one standard deviation, with 40x mesh, while the foamStar-RZ are about 5% smaller. When the mesh becomes even coarser, only the results of foamStar-SWENSE show no obvious degradation of accuracy. For the foamStar-RZ and ISIS-CFD-IWG, the wave force decreases when enlarging the far-field cell size. This decrease is most due to the numerical diffusion errors induced by the mesh in the incident wave propagation. The IWG method and RZ methods rely on the CFD solver to calculate the incident wave propagation and thus require fine cells everywhere in the domain. In contrast, the incident waves of the SWENSE method are solved with external potential flow solvers and are not affected by the CFD mesh. As a result, the SWENSE method is the least sensitive to the quality of far-field mesh.

In Figs. 15(b) and 15(c), the first harmonic amplitudes of the vertical force (F (1) z ) and the free surface elevation at probe 1 (η 1 ) reveal similar results with the horizontal force, i.e., the three methods show close and the most accurate results with 80x mesh, but the SWENSE method is the least sensitive to the far-field mesh quality. However, for the vertical force (F (1) z ), an obvious larger discrepancy between the numerical results and the experimental data is observed. Even with the finest mesh, this difference still remains and exists for all the three meth-Non-dimensional Amplitude 0.070 ± 0.002 0.165 ± 0.002 0.009 ± 0.005 0.056 ± 0.002 0.034 ± 0.003 0.038 ± 0.003 80x (1.4M) 0.070 ± 0.001 0.166 ± 0.001 0.011 ± 0.005 0.059 ± 0.002 0.034 ± 0.002 0.043 ± 0.002 foamStar-SWENSE C1 (0.72M) 0.067 ± 0.001 0.175 ± 0.001 0.018 ± 0.003 0.067 ± 0.002 0.043 ± 0.002 0.051 ± 0.001 C2 (0.82M) 0.065 ± 0.001 0.178 ± 0.002 0.016 ± 0.003 0.068 ± 0.002 0.043 ± 0.002 0.052 ± 0.001 C3 (1.0M) 0.065 ± 0.001 0.175 ± 0.003 0.015 ± 0.003 0.068 ± 0.002 0.038 ± 0.002 0.050 ± 0.002 20x (1.3M) 0.069 ± 0.001 0.177 ± 0.001 0.014 ± 0.002 0.066 ± 0.002 0.044 ± 0.002 0.041 ± 0.001 40x (1.5M) 0.070 ± 0.001 0.184 ± 0.001 0.011 ± 0.002 0.062 ± 0.002 0.040 ± 0.002 0.039 ± 0.001 80x (2.4M) 0.070 ± 0.001 0.185 ± 0.001 0.011 ± 0.002 0.062 ± 0.002 0.039 ± 0.002 0.040 ± 0.001 Table 7: Comparison between numerical results and experimental data: mean force and second harmonics ods. Provided that the heave damping skirt mainly influences the vertical wave force, it is possible that the local refinement near the skirt (identical in all the meshes) is not enough to capture accurately the complex vortex shedding phenomenon.

F (1) x F (1) z η (1) 1 η (1) 2 η (1)
The zeroth harmonic and the second harmonic of the horizontal force is about 5% and 10% to the first harmonic force. These values can give a more detailed comparison of the accuracy. Moreover, the zeroth harmonic is of engineering interest since it corresponds to the drift force of the waves. Figure 16 provides such a comparison. The SWENSE method still shows a good accuracy on both quantities, regardless the mesh used, while the loss of amplitudes is also obvious in the results of the foamStar-RZ and ISIS-CFD-IWG with mesh 20x, due to the numerical damping caused by large mesh size in the far field. The foamStar-RZ achieves good results with fine meshes (40x and 80x), while the ISIS-CFD-IWG underestimates the drift force and has larger periodical uncertainty in its result.

Discussion

Efficiency of the models

The results reveal that each wave modeling method achieves its best accuracy with the finest discretization (80x mesh). However, the use of coarse mesh in the far-field has different influence on their accuracy. The SWENSE method is able to use much coarser mesh by giving little concession on the accuracy and thus is the most efficient wave modeling method among the three, and is very advantage in cases where the efficiency is more emphasized than the absolute accuracy,

To confirm that, we compare the efficiency of the three models according to their mesh requirement to reach a same level of accuracy.

Taking the first harmonic amplitude of the horizontal wave force as an example: if set an accuracy tolerance of 3% (1.35 < F (1) x < 1.43), then ISIS-CFD-IWG requires the 40x mesh with 1.5 million cells to meet the criteria; foamStar-RZ requires the 80x mesh with 2.5 million cells; while foamStar-SWENSE is able to achieve the same level of accuracy with the coarsest C1 mesh with only 0.72 million of cells.

The efficiency enhancement in the SWENSE method is credited to the use of the potential/viscous coupling approach, so that accurate incident waves are independent of the CFD mesh, only the complementary field needs to be solved.

Moreover, this study confirms that the complementary fields need only to be solved accurately in a close vicinity of the structure, so that a small local refinement zone near the structure with coarse mesh in the far-field is able to produce good results: the cylindrical mesh with a local refinement of 0.04 incident wave length around the structure (C1) is enough to achieve a good estimation of horizontal wave force with -2.3% error compared to the experimental data. Although enlarging this local refinement zone to 0.27 wave length wide (C3) can reduce the error to -1.7%, using larger refinement zone is not always necessary if the efficiency is more emphasized than the accu- To ensure the accuracy of the simulation, especially to validate the result of SWENSE on the coarse mesh, the flow details of the simulation are compared. The IWG method uses the medium mesh 40x; the RZ method uses the finest mesh 80x; the SWENSE method uses the coarsest mesh 20x. With such discretizations, all methods should have a same level of accuracy. Figure 17 plots the Q-criterion and the pressure fields obtained by the three methods, showing a good agreement.

Wave reflection prevention ability

The wave reflection is prevented using two strategies in the present study.

• Full absorption condition: damping zones are used to absorb the entire wave trains at the outlet boundary, as in the IWG method and in the experiment. • Diffraction absorption condition: this condition allows the incident waves to leave the computational domain and absorb the diffracted waves only, as in the RZ method and in the SWENSE method.

It is observed that the diffraction absorption condition prevents better wave reflection, comparing with the full absorption condition.

For example, when the full absorption condition is used, both the experiment and numerical methods show difficulty in damping out the full wave trains at the end of the domain. Reflected waves travel back to the structure and affect the results (see point C in Fig. 6 and point B in Fig. 11(b)). As a result, the exploitable window for the data analysis is reduced. A secondary effect of the full absorption condition is that it can not be used to prevent side wall reflections. In the present study, the numerical simulation with the IWG method is not long enough to make this phenomenon appear, but the experimental data confirms the necessity of avoiding side wall reflections.

The RZ method and the SWENSE method both use the diffraction absorption condition at all domain boundaries. With both methods, no apparent wave reflections are observed, the simulations are more stable and provide longer exploitable time for the data analysis.

Inspired by the RZ method, forcing zones are added to the IWG method to achieve the diffraction absorption condition. The forcing zone is similar to the RZ, but it uses a momentum source in the RANS equations to force the CFD solution equal to the incident wave solution, and thus it provides the absorption of the diffracted waves. In Fig. 18 the results of IWG using forcing zones at side walls and the outlet are plotted by red dash lines and compared to the original method shown by black lines. The forcing zones avoid the wave reflections and show a much more stable result. In the present numerical study, the standard k -ω SST turbulence model is selected to calculate the eddy viscosity in the RANS equations. However, an excessive wave damping and a wave force diminishing are observed with the RZ method implemented in OpenFOAM. Figure 19 shows the 0th and 1st harmonic amplitudes of the horizontal wave force. A loss in amplitude is observed from t = 15s . This problem is purely numeric and limits the effective time window to analyze the physics.

This problem has been reported in the literature [START_REF] Devolder | Application of a buoyancymodified k-ω SST turbulence model to simulate wave run-up around a monopile subjected to regular waves using OpenFOAM R[END_REF]) on using RANS model for two-phase flow simulations and it is attributed to the limitation of the present turbu-lence models, which are mostly derived and tuned for singlephase flows. These models are usually not adequate to consider the discontinuity of physics near the free surface and create excessive and nonphysical turbulent eddy-viscosity, which leads to the wave damping. This problem is often severe in a wave-structure interaction problem without forward speed, since the turbulent viscosity always remains in the computational domain. Figure 20 shows the turbulent viscosity (ν t ) at t = 1T and t = 15T in the simulation with the RZ method. After one wave periods, the turbulent viscosity appears not only close to structure but also in the far-field. The far-field turbulence is created near the free surface, which is shown as a white line, due to the interface discontinuity. The maximum ν t has an order of 10 -3 m 2 /s at t = 1.8s. At the end of the simulation, the turbulent viscosity has a much larger value (order of 10 -1 m 2 /s) due to the accumulation and appears in a larger zone.

The result of the SWENSE method is not affected by this problem. Despite using the same turbulence model, the wave force amplitudes are maintained constant even after a long simulation time. This has to be credited to the potential/viscous coupling used in SWENSE method, so that the viscous effects is only applied to the complementary field and does not affect the incident waves.

Conclusion

The present work compared three wave modeling methods in two-phase CFD solvers: the IWG method in ISIS-CFD, the RZ method and the SWENSE method in foamStar. The accuracy and the performance of these methods are compared by simulating a fixed CALM buoy in regular waves and compared with experiments.

On the accuracy side, all the three methods are able to achieve accurate comparison with the experiment when the mesh is fine enough to propagate the incident waves accurately. Compared to the SWENSE method in foamStar, the IWG method in ISIS-CFD and the RZ method in foamStar show more dependency on fine meshes to obtain accurate results. Only the SWENSE method in foamStar still obtain good results with coarse far-field mesh.

On the performance side, the SWENSE method in foamStar is the most efficient. It is able to achieve a significant gain in efficiency with a tiny concession on the accuracy, by using coarse mesh in the far-field. For a same level of accuracy, the SWENSE method requires only about 25% cells compared to the RZ method, with both methods implemented in Open-FOAM.

The diffraction absorption condition prevents better wave reflection than the full absorption condition. The RZ and the SWENSE method in foamStar, using the diffraction absorption condition, are able to prevent reflections at all the domain boundaries, and thus provide more stable result and longer time window for the data analysis. The full absorption method (damping), used by the IWG method in ISIS-CFD and the experiment, creates larger reflection at the outlet and is inapplicable for side wall absorption. Using forcing zones is able to achieve the diffraction absorption condition in the IWG method.

A last remark shows the k -ω SST turbulence model is not adapted for two-phase flow simulation in standard RANSE solvers, since it causes excessive numerical damping of incident waves. The use of the potential/viscous coupling method as the SWENSE method is able to alleviate this problem and maintain good simulation results for longer simulation time.

  Relaxation zone (RZ) method

Figure 1 :

 1 Figure 1: Wave generation and absorption techniques

Figure 2 :

 2 Figure 2: The SWENSE method decomposes the total field into an incident and a complementary part

  Figure 3: CALM buoy geometry

  19m away from the wavemaker to indicate the far-field wave information.

Figure 4 :

 4 Figure 4: Wave probes near the buoy (top view)

Figure 5 :

 5 Figure 5: Free surface elevation at probe No.4 in the experiment

  Figure 7: Rectangular mesh configurations

Figure 8 :

 8 Figure 8: Cylindrical mesh configuration

Figure 10 :

 10 Figure 10: Local refinement regions near the buoy in the cylindrical configuration

Figure 11 :

 11 Figure 11: Horizontal wave force: internal wave generation method in ISIS-CFD

Figure 12 :

 12 Figure 12: Horizontal wave force: RZ method in foamStar

Figure 14 :

 14 Figure 13: Horizontal wave force: SWENSE method in foamStar with the cylindrical mesh

  Figure 15: Comparison of simulation results: first harmonics of the wave forces and the free surface elevation

Figure 16 :

 16 Figure 16: Comparison of simulation results: the mean horizontal wave force and the second harmonic of the horizontal force

Figure 18 :

 18 Figure 18: Influence of the absorption boundary condition for the IWG method: BC1 wave damping at outlet, and BC2 forcing zones at outlet and at side wall

Figure 19 :

 19 Figure 19: Horizontal wave force result in a long time simulation

  (a) One wave period after the simulation (b) 15 wave periods after the simulation

Figure 20 :

 20 Figure 20: Turbulent viscosity in the computational domain (view from the symmetry plane). Note that the color map uses a logarithmic scale. The free surface is shown as a white line.

Table 1 :

 1 Geometry parameters of the CALM buoy model

	Parameter	Value
	Period (T )	1.80 s
	Height (H)	0.16 m
	Depth (h)	5.00 m
	Steepness (ka)	0.10
	Relative water depth (kh)	6.20
	Table 2: Target incident wave information
	3.2. Data measurement and analysis
	3.2.1. Data measurement

Table 3 :

 3 Experimental results

	Figure 6: The first harmonic amplitude of the wave force in the exper-
	iment

Table 4 :

 4 Rectangular mesh information

	Mesh λ/∆x λ/∆y H/∆z Number of cells
	20x	20	10	16	1.3 M
	40x	40	10	16	1.5 M
	80x	80	20	16	2.5 M

Table 6 :

 6 Comparison between numerical results and experimental data: first harmonics

	3

https://lheea.ec-nantes.fr/software-and-patents/isis-cfd-193387.kjsp

https://www.numeca.com/product/finemarine

http://www.openfoam.com

Acknowledgments

This work has been performed in the framework of the Chaire Hydrodynamique et Structure Marines CENTRALE NANTES -BUREAU VERITAS. Part of this work was granted access to the HPC resources of IDRIS and CINES under the allocations 2018-A0032A01308 made by GENCI (Grand Équipement Na-tional de Calcul Intensif). The first author acknowledges China Scholarship Council (CSC) for the financial support for his PhD study.


Atul D Sontakke

Lisa I D J Martin

Victor Castaing

Bruno Viana

Philippe F Smet

Cathodoluminescence and microstructural analysis of amorphous Yttrium-Aluminum-Borate luminescent powders

We present the cathodoluminescence (CL) and microstructural analysis of amorphous yttrium aluminum/aluminium borate (a-YAB) luminescent powders synthesized by polymeric precursor (PP) and sol-gel (SG) based wet-chemical methods. The a-YAB powders exhibit bright photoluminescence tunable from intense blue to white chromaticity region depending on their calcination temperature. Here, the influence of calcination temperature on CL and microstructural properties of the a-YAB powders has been investigated. The a-YAB powders showed weak but homogeneous CL. Within the optimum calcination temperatures for photoluminescence (650-740°C for PP and 450-600°C for SG series), the CL observed to be more intense for powders synthesized at higher calcination temperature. The microstructure analysis revealed that the SG series powders calcined below 570°C exhibit intraparticle chemical inhomogeneity, whereas no such inhomogeneity could be seen in PP series powders, owing to their higher calcination temperatures. Nevertheless, both series exhibit similar luminescence broadening with calcination temperatures despite these microscopic differences, revealing that the luminescence broadening does not depend on the chemical properties at the microscopic level, but is controlled by more intrinsic effects on the luminescent centers or its immediate surrounding as a function of calcination temperature.

Introduction

The wet-chemically synthesized amorphous yttrium aluminum borate (a-YAB) powders, with a chemical composition close to that of crystalline YAl3(BO3)4, exhibit bright photoluminescence (PL) [START_REF] Sontakke | Afterglow luminescence in wet-chemically synthesized inorganic materials : ultra-long room temperature phosphorescence instead of persistent luminescence[END_REF]. The origin of this luminescence has recently been ascribed to the presence of organic molecules trapped inside the inorganic glassy matrix during the synthesis process and their subsequent evolution with thermal treatment [START_REF] Sontakke | Afterglow luminescence in wet-chemically synthesized inorganic materials : ultra-long room temperature phosphorescence instead of persistent luminescence[END_REF][START_REF] Burner | Evidence of organic luminescent centers in sol-gel-synthesized yttrium aluminum borate matrix leading to bright visible emission[END_REF]. The powders exhibit intense blue to broadband white emission depending on their thermal history as well as their synthesis methods (Figure 1). We have recently reported the polymeric precursor (PP) and sol-gel (SG) based wet-chemical methods to synthesize micrometer sized amorphous a-YAB powders [START_REF] Sontakke | Afterglow luminescence in wet-chemically synthesized inorganic materials : ultra-long room temperature phosphorescence instead of persistent luminescence[END_REF][START_REF] Burner | Evidence of organic luminescent centers in sol-gel-synthesized yttrium aluminum borate matrix leading to bright visible emission[END_REF][START_REF] Guimaraes | Towards a new generation of white phosphors for solid state lighting using glassy yttrium aluminoborates[END_REF]. In PP method syntheses, relatively high calcination temperatures (650-740°C) are necessary to achieve optimum luminescence over the SG method syntheses (450-600°C). Nevertheless, in both series a similar luminescence broadening is witnessed within the investigated temperature ranges.

To further develop this type of luminescence in the framework of luminescent probes [START_REF] Soon | Engineering persistent luminescence nanoparticles for biological applications: from biosensing/bioimaging to theranostics[END_REF][START_REF] Feng | A review of mechanoluminescence in inorganic solids: compounds, mechanisms, models and applications[END_REF] or light emitting devices [START_REF] Kim | Review-phosphor plates for high-power LED applications: challenges and opportunities toward perfect lighting[END_REF], it is imperative to have a good understanding on the origin of the luminescence centers and its broadening phenomenon. While the broad emission spectrum is advantageous for the use of these phosphors in solid state lighting devices with high color rendering, the rather strong emission in the violet to blue part of the visible spectrum limits the overall luminous efficacy [START_REF] Smet | Selecting conversion phosphors for white light-emitting diodes[END_REF]. It is therefore vital to be able to control and optimize the luminescence properties. In contrast to luminescence based on lanthanides (e.g. Eu 2+ , Ce 3+ or the 4f-4f emitting ions) [START_REF] Dutczak | Eu 2+ luminescence in strontium aluminates[END_REF][START_REF] George | Local environments of dilute activator ions in the solid-state lighting phosphor Y3-xCexAl5O12[END_REF][START_REF] Van Hest | Probing the Influence of disorder on lanthanide luminescence using Eu-doped LaPO4 nanoparticles[END_REF] or transition metals (e.g. Cr 3+ or Mn 4+ ) [START_REF] Borisov | Preparation and characterization of chromium(III)-activated yttrium aluminum borate: a new thermographic phosphor for optical sensing and imaging at ambient temperatures[END_REF][START_REF] Sijbom | K2SiF6:Mn 4+ as a red phosphor for displays and warm-white LEDs: A review of properties and perspectives[END_REF], where the emission spectrum can roughly be predicted based on the type of transition and the local environment of the dopant, this structural-luminescence relationship has not yet fully been established in the a-YAB powders, which lack a traditional dopant ion. Moreover, due to the inorganic refractory-type nature of the a-YAB host material, it is difficult to extract the trapped organic species in intact form for chemical analysis [START_REF] Patarin | Mild methods for removing organic templates from inorganic host materials[END_REF][START_REF] Hong | Methods of analysing chemicals associated with microplastics: A review[END_REF]. The complementary solid-state analytical techniques (including NMR, EPR, ENDOR, etc.) revealed the presence of small polyaromatic hydrocarbons (PAH), which exhibit good correlation with the broadening of the luminescence spectral profiles [START_REF] Burner | Evidence of organic luminescent centers in sol-gel-synthesized yttrium aluminum borate matrix leading to bright visible emission[END_REF]. A parallel approach using time-gated spectroscopy and energy level modeling hinted to the presence of multiple emitting centers, which undergo a systematic evolution in relative strengths as a function of calcination temperature, thereby causing the spectral broadening [START_REF] Sontakke | Time-gated triplet-state optical spectroscopy to decipher organic luminophores embedded in rigid matrices[END_REF]. Each emitting center exhibits fluorescence as well as phosphorescence components, and therefore the complete spectrum becomes an admixture of several emission bands. Based on the time-gated spectroscopy and theoretical modelling using DFT calculations, we could construct the energy level structure of first few prominent emitting centers and identify the higher energy violet-blue region emitting center to be a phenalene based compound [START_REF] Sontakke | Time-gated triplet-state optical spectroscopy to decipher organic luminophores embedded in rigid matrices[END_REF]. The results also indicated that the primary energy level structure of the other luminescent centers might be similar but with decreasing energy gaps that lead to the red shifted emissions. This suggested that the broadening might be associated with the bathochromic effect in luminophore compounds. The bathochromic spectral shift can arise due to a π-conjugation in organic molecules, a change in the functional group, or a variation in chemical environment [START_REF] Meier | Extension of conjugation leading to bathochromic or hypsochromic effects in OPV series[END_REF]. During the calcination step in wet-chemical synthesis methods, the primary network structure undergoes densification. Simultaneously, the solid-gas reactions take place inside the powder grains through diffusion, leading to complex molecular reactions. These interactions during calcination process can modify the luminescent compounds as well as chemical properties of the host, responsible for the luminescence properties of the a-YAB powders.

In the present study, the microstructural properties of a-YAB powders as a function of calcination temperatures have been investigated using scanning electron microscopy and cathodoluminescence (CL) in order to seek a correlation between the microstructure and the luminescence properties. CL spectroscopy is a vital technique to unravel the hidden details of functional luminescent materials at microscopic levels, which collectively represents their bulk performances [START_REF] Kuchařová | Microscopy and cathodoluminescence spectroscopy characterization of quartz exhibiting different alkali-silica reaction potential[END_REF][START_REF] Dierre | Role of Si in the luminescence of AlN:Eu,Si phosphors[END_REF][START_REF] Martin | Microscopic study of dopant distribution in europium doped SrGa2S4: Impact on thermal quenching and phosphor performance[END_REF][START_REF] Abe | Hybrid remote quantum dot/powder phosphor designs for display backlights[END_REF][START_REF] Joos | Origin of saturated green emission from europium in zinc thiogallate[END_REF][START_REF] Broxtermann | Cathodoluminescence and photoluminescence of YPO4: Pr 3+ , Y2SiO5: Pr 3+ , YBO3: Pr 3+ and YPO4: Bi 3+[END_REF]. A scanning electron microscope (SEM) was used to study the microstructure evolution, while simultaneously recorded energy dispersive X-ray (EDX) elemental mapping and CL microscopy provided further insights on the structure-property relationship.

Experimental details

Two series of a-YAB powders obtained by PP (calcination temperatures : 650 -740°C) and SG (calcination temperatures : 450 -600°C), based on wet-chemical syntheses were used in the present study. 1,2 The chemical composition of both series of a-YAB powders is stoichiometrically close to crystalline YAl3(BO3)4 (c-YAB). The powders are labeled as PP-x and SG-x, where x represents the final heat-treatment (calcination) temperature in oxidizing atmosphere in degrees Celsius. The powders exhibit a clear white to pale yellow colored tint with the increase in calcination temperature in both series. The Y3Al5O12:Ce 3+ (YAG:Ce), ZnGa2O4:Cr 3+ (ZGO:Cr) polycrystalline phosphors and Ce 3+ : borate glass were also used for comparison. The detailed synthesis procedure for the studied samples is described in the supporting information (SI) and in literature [START_REF] Sontakke | Afterglow luminescence in wet-chemically synthesized inorganic materials : ultra-long room temperature phosphorescence instead of persistent luminescence[END_REF][START_REF] Burner | Evidence of organic luminescent centers in sol-gel-synthesized yttrium aluminum borate matrix leading to bright visible emission[END_REF][START_REF] Sontakke | A comparison on Ce 3+ luminescence in borate glass and YAG ceramic: Understanding the role of host's characteristic[END_REF]. The phase purity of polycrystalline YAG:Ce and ZGO:Cr, as well as the amorphous nature of a-YAB powder was confirmed using X-ray diffraction analysis (supporting information, Figure S1)

The PL spectra of the powders were obtained under 254 nm excitation and the emission signals were recorded using a calibrated CCD detector (Acton Pixis 100, Princeton Instruments). A similar setup was used for X-ray radioluminescence (RL), where 1.54 Å Cu Kα X-ray was used as excitation. For CL spectroscopy, a Hitachi S3400-N scanning electron microscope was used, combined with a fiber-coupled spectrometer (Acton SP2300 monochromator and ProEM 1600 EMCCD camera, Princeton Instruments). The electron microstructure was studied using backscattered electrons. Energy dispersive X-ray imaging of the Y Lα (1.922 keV), Al Kα (1.486 keV) and O Kα (0.525 keV) lines was used for the elemental analysis at the microscopic level. The SEM was operated at 20kV potential and 50 μA current density under 25 Pa pressure to avoid charging of the samples. The samples were prepared by spreading fine powder particles on a carbon tape attached to the sample mount.

Results and Discussion

Figure 2 presents the SEM images of the microstructure and the related CL emission of SG-600 a-YAB powder. All particles exhibit CL emission and the particles' morphology could clearly be reproduced in the CL mapping. The PP series powders also exhibit homogeneous CL emission [START_REF] Guimaraes | Controlled preparation of aluminum borate powders for the development of defectrelated phosphors for warm white LED lighting[END_REF]. It has been observed that the CL signals are weak in a-YAB powders compared to the traditional inorganic phosphor materials. On comparing the CL response with respect to the traditional crystalline phosphors, such as ZnGa2O4 : Cr 3+ (ZGO:Cr), the CL emission intensity in a-YAB powders is several orders of magnitude lower (Figure 3a). The intensity is even weaker for the a-YAB powders synthesized at lower calcination temperatures (Supporting Information, Figure S2). The low intensity for CL signals in a-YAB powders may be due to the amorphous nature of the a-YAB host. Amorphous or glassy hosts are usually less efficient in trapping ionizing radiations and delivering the excitation energy to the luminescent active centers [START_REF] Sontakke | A comparison on Ce 3+ luminescence in borate glass and YAG ceramic: Understanding the role of host's characteristic[END_REF]. The presence of several defects, including the dangling bonds, various point defects, etc., are inherent to the amorphous structure that can easily facilitate a swift non-radiative recombination of the excitonic e-h (electron-hole) pairs created by high energy irradiation. Moreover, the organic luminophores exhibit poor response to high energy ionizing excitations, and even undergo degradation under continuous exposure (Supporting Information, Figure S3), though the respective PL signals are stable with respect to the excitation duration and temperature (Supporting Information, Figures S4,S5). Similarly to the CL, the X-ray radioluminescence (RL) intensity is weak in a-YAB powders compared to the traditional lanthanide doped luminescent materials. Figure 3b presents the RL response of PP-740 a-YAB powder compared with the polycrystalline YAG:Ce phosphor and Ce-doped borate glass. It is interesting to see that the YAG:Ce exhibits almost seven orders of magnitude more intense RL over a-YAB, whereas two orders of enhancement is observed in Ce 3+ doped glass. The intensity difference also reflects the polycrystalline and glassy nature of the respective host materials, i.e. YAG and borate glass, respectively. From Figure 3, it can be noticed that the CL spectral profile of PP-740 a-YAB powder exhibits good resemblance with its PL spectral profile (PP-740 powder in Figure 1), whereas the RL profile is significantly different. The dominant blue region emission in the PL has slightly reduced in CL, but it is significantly reduced in the RL spectrum. This is consistent with the fact that the a-YAB powder exhibits several emitting centers with varied energy levels and the relatively deeper penetration depth of the used X-ray radiation (several tens of μm) might induce radiation trapping of higher energy (violet-blue region) emissions in the RL emission signals [START_REF] Sontakke | Time-gated triplet-state optical spectroscopy to decipher organic luminophores embedded in rigid matrices[END_REF][START_REF] Guimaraes | Controlled preparation of aluminum borate powders for the development of defectrelated phosphors for warm white LED lighting[END_REF]. The penetration depth in CL is expected to be of the order of 3-4 μm at the used electron energy in a-YAB host, so the reabsorption is relatively weaker (Supporting Information, Figure S6).

The multiple emissions in a-YAB powder arise from different emitting centers, whose relative abundance is controlled by the calcination temperature [START_REF] Sontakke | Time-gated triplet-state optical spectroscopy to decipher organic luminophores embedded in rigid matrices[END_REF]. This causes the broadening of predominantly intense violet-blue region emission as the calcination temperature increases. Simultaneously, the longer wavelength emissions become more distinctive at the expense of violet-blue emissions. Although this variation in the relative prominence of different emitting centers could be observed experimentally, the origin of this systematic bathochromic shift in emission is yet to be understood. As discussed earlier, it can either be due to the transformation of emitting centers, such as the π-conjugation of PAH compounds or due to the changes in functional groups, or can also be due to the changes in the chemical environment of the surrounding host matrix. An enlarged π-conjugation reduces the energy of electronic transitions in PAH compounds. Similarly, the changes in functional group or the chemical polarization of surrounding host matrix can modify the energy levels of the emitting centers. The microscopic analysis, such as the powder microstructure and cathodoluminescence could be interesting in this regard to understand the evolution of the a-YAB matrix with temperature and find any relation with luminescence broadening. Figure 4 presents the SEM microstructure images of the PP and SG series powders calcined at various temperatures. For the entire PP series, the images do not exhibit noticeable contrast difference amongst particles suggesting a near homogeneous elemental distribution (Supporting Information, Figure S7) [START_REF] Berz | An Introduction to Beam Physics[END_REF]. However, the SG series powders show interesting microstructural evolution as the synthesis temperature increases. For powders prepared at low temperature, the microstructure reveals particles with dark and bright regions. As the dark/bright color difference in SEM imaging based on backscattered electrons mainly originates from the difference in average atomic weights of constituent elements [START_REF] Lloyd | Atomic number and crystallographic contrast images with the SEM: A review of backscattered electron techniques[END_REF], this suggests an elemental inhomogeneity. The darker particles, which have a relatively lower average weight, slowly vanish as the calcination temperature increases. For a calcination temperature of 570°C, nearly homogeneously bright particles are obtained with the SG synthesis. Figure 5 presents the CL emission spectra of some PP and SG powders. Similar to PL, the CL emission also exhibits spectral broadening and a redshift in the peak position with increasing calcination temperature. Figure 6(a) presents the EDX elemental analysis of particles in the SG-480 powder. It suggests that the dark region particles are exclusively yttrium deficient, whereas aluminum seems well distributed in all the particles. These observations suggest that the dark regions are enriched in boron, associated with aluminum, however boron, due to its low atomic weight, could not directly be detected by EDX. It is worth noting here that despite the presence of such chemical differences, the CL emission profile does not exhibit significant deviation among different particles and peaks at nearly same positions in all the particles of respective samples (Supporting Information, Figures S8,S9). Nevertheless, a more reliable comparative analysis of CL emission profile of different particles is limited due to the different levels of signal to noise ratio, varied penetration depths, exposure based degradation and overall weak CL signals. Interestingly, at 570°C and higher calcination temperatures, particles do not exhibit such distinct elemental inhomogeneity (Figure 4 and Supporting Information, Figure S7). The gradual disappearance of dark regions in particles is accompanied with a systematic appearance of bright spots, whose abundance and dimensions grows for higher calcination temperature. Figure 6(b) shows the elemental mapping of such a particle exhibiting bright region spots representing yttrium rich regions. This points at an increase of diffusion dynamics of atom species within the particles when the calcination temperature increase. Accordingly, the yttrium rich regions appear on the surface of the dark particles, and start to grow bigger with the increase in calcination temperature (Figure 7a,b). Eventually, the bright spots start to coalesce and lead to completely bright particles with near homogeneous elemental distribution (Figure 7c-e). To get further insight on the diffusion reaction, the particles were crushed to see the cross-sectional EDX elemental distribution. Figure 8 shows the microstructure as observed in SEM and EDX mapping of darker particles after gently crushing them. It is clear that mainly the surface layer of the dark region particles is deficient in yttrium. The low density aluminumborate matrix forms the surface region, whereas the dense yttrium-aluminate matrix is segregated in the core region. With increase in calcination temperature, the elements diffuse more effectively leading to the homogeneous yttrium-aluminum-borate matrix. It is evident that the presence of darker particles gradually decreases from SG-450 to SG-570, where the SG-570 does not contain any difference between or within particles. The chemical composition (Al:Y) for SG-570, SG-600 and all the PP-x series powders is similar to c-YAB suggesting a nearly homogeneous microscopic elemental distribution in these powders. Indeed, above 570°C nearly homogeneous particles are obtained with a Y:Al ratio about 1.1 close to that of the crystalline YAB composition (Supporting Information, Figure S10). Consequently, the SG-x series powders show significant structural variations, whereas the PP-x series powders are essentially homogeneous. For both series of powders, the PL and CL emission spectra exhibit spectral broadening as a function of calcination temperature. As the calcination temperature primarily controls the elemental diffusion, an interesting microstructural evolution in SG series a-YAB powders could be observed. But any direct relation between powder microscopic properties and the PL band broadening could not be established. It is therefore evident that the luminescence broadening as a function of calcination process is more intrinsic to the luminescent centers rather than the microscopic level chemical properties of a-YAB powders. It might be controlled by immediate changes in and around the luminescent centers brought about due to thermal treatment during calcination process. The chemical surrounding at immediate neighborhood of emitting PAH compounds can be decisive in this respect, which may be limited to a few atomic layers and therefore could not be resolved in present investigation. Moreover, the other effects such as pi-conjugation enlargement and the effect of functional groups could also be responsible for the bathochromic shift and needs further investigation.

Conclusions

In summary, the a-YAB powders exhibit weak but homogeneous cathodoluminescence emission. The microstructural analysis revealed that the sol-gel synthesized a-YAB powders exhibit elemental inhomogeneity at low calcination temperature (<570°C) and systematically attain a homogeneous a-YAB matrix composition with increase in calcination temperature. This inhomogeneity could not be observed in polymeric precursor synthesized powders. Despite these structural differences, the PL and CL spectra for both syntheses showed a similar spectral broadening with increasing calcination temperature within the respective optimum temperature range. These results suggest that the microscopic level elemental homogeneity and its diffusion dynamics do not exhibit significant influence on the spectral broadening of the luminescence in a-YAB powders during the calcination process, but that this might be controlled by more intrinsic effects in and around the luminescent centers.

Figure 1 :

 1 Figure 1: PL spectra of a-YAB powders obtained by sol-gel (SG) and polymeric precursor (PP) methods. The PL were recorded under 254 nm excitation to achieve complete spectral profiles. The powders are denoted by SG-x and PP-x, where x represents the calcination temperature in °C for both studied synthesis routes.

Figure 2 :

 2 Figure 2: SEM microstructure (a) and corresponding CL intensity map (b) of SG-600 powder.

Figure 3 :

 3 Figure 3: CL spectra (a) and X-ray RL spectra (b) comparison between PP-740 a-YAB powder and representative traditional luminescent materials.

Figure 4 :

 4 Figure 4: Influence of calcination temperature on the microstructure of a-YAB powder in the PP (top) and SG series (bottom), when imaged in the electron backscattered imaging SEM mode.

Figure 5 :

 5 Figure 5: CL emission spectra of PP-series and SG-series a-YAB powders.

Figure 6 :

 6 Figure 6: (a) SEM microstructure and EDX elemental mapping of the elements Al, Y and O in SG-480 powder. (b) The enlarged view of a particle possessing both dark and bright regions. A more saturated color represents a higher density of the elements.

Figure 7 :

 7 Figure 7: SEM microstructure of a-YAB powder particles exhibiting different stages (a-e). (a) is from SG-450 powder, (e) is from SG-570 powder and (b-d) are from intermediate temperature calcined powders. Scale bars represent 25 µm.

Figure 8 :

 8 Figure 8: Microstructure (SEM image, top left) and EDX elemental maps of the elements Y, Al and O in a crushed dark particle. A more saturated color represents a higher density of the elements.

Acknowledgements:

This work was supported by the ANR LuminoPhorLED (ANR-14-CE05-0033) and Lumicor SBO130030 project grants. Authors would like to thank Prof. Isabelle Gautier-Luneau and Dr. Alain Ibanez, CNRS-Inst. NEEL, Grenoble for providing the samples and fruitful suggestions.