Dorian Pena
email: dorian.pena@etu.unistra.fr

Yannick Hoarau
email: hoarau@unistra.fr

Eric Laurendeau
email: eric.laurendeau@polymtl.ca

A Single Step Ice Accretion Model Using Level-Set Method

Keywords: Aircraft Icing, Level-Set Method, Eulerian Droplet Transport Equations, PDE Icing Model 5

A new approach using the Level-Set framework is developed in the NSMB (Navier-Stokes Multi-Block) compressible solver for modeling the ice/air interface evolution through time during in-flight icing. Droplet distribution and impingement efficiency are computed by an Eulerian approach and the accreted ice is calculated by a PDE model. An icing velocity field is introduced and the Level-Set equations are solved on body-fitted multi-block structured grids. The whole process is parallelized with the MPI library for efficient calculations. Single step icing is simulated on NACA23012 and NACA0012 airfoils and on the ONERA-M6 and the GLC-305 swept wings. In all the studied cases, the results are in good agreement with existing and available data validating the feasibility of the approach.

Introduction

Ice accretion and its effects represent a wide research area within of the aeronautic community. According to the US National Transport Safety Board, icing is one of the major causes of flight accidents [START_REF] Reehorst | National transportation safety board aircraft accident investigation supported[END_REF]). By impacting on aircraft structures, the supercooled water droplets present in a significant quantity in clouds can freeze or flow as a water film depending on the surface temperature. Two main types of ice known as rime ice and glaze ice are formed depending on the environmental conditions. The dry rime ice usually happens at low flow speeds and low flow temperatures when the entire droplet freezes at the impact point. When the flow temperature and the flow speeds are higher, wet glaze ice is usually formed. In that case, a fraction of the droplet freezes at impact point and the other portion flows as a thin water film runback. Typical forms of rime and glaze icing are shown in Fig. 1. The modification of the geometry induced by icing often results in large aerodynamic degradation. Because of a stringent certification process and since all configurations cannot be experimentally assessed in wind tunnels, aircraft in-flight icing needs to be simulated via computational fluid dynamics.

Ice accretion simulation codes used by industries for the certification process are, still today, commonly based on 2D inviscid panel methods to simulate the fluid, a boundary layer method based on empirical correlations [START_REF] Kays | Convective heat and mass transfer[END_REF]) to evaluate the heat transfer coefficient , a Lagrangian formulation to predict the droplets trajectories and the so-called Messinger model to compute the ice thickness [START_REF] Messinger | Equilibrium temperature of an unheated icing surface as a function of air speed[END_REF]). At this point a new geometry is obtained and this one-way coupling process can be repeated until the desired exposure time is reached (multi-step approach) [START_REF] Hasanzadeh | Quasi-steady convergence of Multistep Navier-Stokes Icing Simulations[END_REF]).

In recent years, icing codes based on the solutions of the Reynolds-Averaged Navier-Stokes (RANS) equations have been developed such as LEWICE3D (NASA) [START_REF] Bidwell | Users manual for the NASA Lewis three-dimensional ice accretion code[END_REF]),FENSAP-ICE (McGill University) [START_REF] Beaugendre | Development of a second generation in-flight icing simulation code[END_REF]), CANICE2D-NS ( École Polytechnique Montréal) [START_REF] Hasanzadeh | Quasi-steady convergence of Multistep Navier-Stokes Icing Simulations[END_REF]), ONICE3D (ONERA) [START_REF] Montreuil | ECLIPPS: 1. three-dimensional CFD prediction of the ice accretion[END_REF]). Because of a seamless integration to Navier-Stokes solvers, an Eulerian formulation for droplet trajectories initially introduced by [START_REF] Scott | Navier-stokes solution to the flowfield over ice accretion shapes[END_REF]) have been implemented in several icing codes [START_REF] Bourgault | A finite element method study of Eulerian droplets impingement models[END_REF]; [START_REF] Jung | Numerical modeling for eulerian droplet impingement in supercooled large droplet conditions[END_REF]; [START_REF] Kim | Ice accretion modeling using an eulerian approach for droplet impingement[END_REF]; [START_REF] Montreuil | ECLIPPS: 1. three-dimensional CFD prediction of the ice accretion[END_REF]; [START_REF] Zhu | 3D ice accretion simulation for complex configuration basing on improved messinger model[END_REF]; [START_REF] Cao | Numerical simulation of ice accretions on an aircraft wing[END_REF]; [START_REF] Jung | Development of eulerian droplets impingement model using hllc riemann solver and pod-based reduced order method[END_REF]). Benefits of an Eulerian approach are numerous, including the direct determination of the impingement efficiency without requiring a statistical averaging process, the computation of the droplet state vector in the entire domain without defining appropriate initial positions of the water droplets and the possibility to compute the droplet field on the same mesh used for the flow.

The traditional 1D Messinger Model have been updated to 2D using PDE formulation by [START_REF] Myers | Extension to the Messinger model for aircraft icing[END_REF] and [START_REF] Beaugendre | FENSAP-ICE's threedimensional in-flight ice accretion module: ICE3D[END_REF]. An iterative algebric 2D formulation of the Messinger model has been developed as well by [START_REF] Zhu | 3D ice accretion simulation for complex configuration basing on improved messinger model[END_REF] and implemented in a previous version of NSMB. One of the advantage of the iterative formulations from Beaugendre et al. and Zhu et al. is that it does not require the explicit detection of the stagnation points.

In these 2D formulations the water film properties are mainly governed by the friction and the pressure gradient. In FENSAP-ICE and ONICE3D the heat transfer coefficient is determined with the Navier-Stokes flow solver replacing the boundary layer method historically used in icing simulation. Mesh deformation methods have been employed to track the ice/air interface in an automatic multi-step icing calculations by [START_REF] Beaugendre | FENSAP-ICE's threedimensional in-flight ice accretion module: ICE3D[END_REF] and [START_REF] Montreuil | ECLIPPS: 1. three-dimensional CFD prediction of the ice accretion[END_REF]. Grid deformation techniques are usually used because of their easiness of implementation in an automatic process of icing simulation. Therefore, their applicability strongly depends on the type of simulated ice and the exposure time. The ice accretion at the leading edge of an airfoil induces cells stretching and thus results in a lower cells density on the profile and a lower accuracy of the simulated turbulent flow. Another common difficulty of such methods when simulating glaze ice shapes is the generation of cells intersections and negative volumes requiring a specific treatment. In LEWICE, a procedure detects if two segments of the updated body surface intersect. If an intersection is found the corresponding points are removed in the revised geometry. Consequently, in most icing simulation cases, mesh deformation algorithms fail quickly to maintain a good quality of the grid especially in the boundary layer region. For three-dimensional configurations the difficulties are even more serious requiring a time-consuming and painful manual remeshing. There exists also some body of literature on the use of global remeshing techniques for multi-steps icing calculations [START_REF] Kreeger | SmaggIce 2.0: additional capabilities for interactive grid generation of iced airfoils[END_REF]; [START_REF] Koomullil | Iced airfoil simulation using generalized grids[END_REF]).

The present study is an effort to develop a new method of icing simulation in order to avoid these limitations and improve the models. Among the panel of available methods, the Level-Set framework introduced by [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF] reveals to be the most promising alternative to the mesh deformation techniques. In a Level-Set framework, a passive scalar is defined equal to zero at the interface, positive outside and negative inside. The Level-Set method have been used by [START_REF] Beaugendre | Computation of ice shedding trajectories using cartesian grids, penalization, and level sets[END_REF] to compute the ice shedding trajectories. [START_REF] Beaugendre | Computation of ice shedding trajectories using cartesian grids, penalization, and level sets[END_REF] consider a constantvolume block of ice and apply a penalisation method around the level-set to simulate the wall. In our view, the level zero could represent the time evolution of the ice/air interface, the inside zone could be seen as the accreted ice allowing, for example, the resolution of heat transfer equation in this domain coupled with an external flow computation at the outside zone. When simulating multi-steps icing, a benefit of using the Level-Set method to track the evolution of the ice/air interface is that the output contour of an icing step can be directly used as input for mesh generators programs without needing additional treatment. The Level-Set method can also follow with ease shapes that change topology allowing for example the simulation of ice cracking and ice shedding which constitute recurrent issues for safety in aircraft design. Finally, one can avoid manual remeshing to track the air/ice interface by using adaptive mesh refinement and penalization method around the Level-Set interface. In such framework the use of wall functions may be useful to model the turbulent boundary layer. Future work will investigate these concepts.

In a previous paper [START_REF] Pena | Development of a threedimensional icing simulation code in the NSMB flow solver[END_REF]), a three-dimensional icing code has been developed and validated in the NSMB (Navier-Stokes-Multi-Block) flow solver. The icing code implemented consists of an Eulerian droplet transport equations and an iterative Messinger model. A mesh deformation algorithm previously implemented in NSMB was used to track the ice/air interface through time for multi-steps ice accretion simulations. In this paper, a novel approach for tracking the ice/air interface using the Level-Set method is introduced and the feasibility of the approach is demonstrated. The Level-Set equations are solved on structured body-fitted grids and only single step icing calculation are performed. The latest updates to the icing code are discussed in sections 2.2 & 2.3, including a droplet module for collection efficiency calculation solved with a total variation diminishing (TVD) scheme via deferred correction and the implementation of a PDE icing model for the calculation of the mass of ice on the body surface. The Level-Set framework developed for this study is detailed in section 2.4. The entire process is parallelized with Message Passing Interface (MPI) for complex multi-block configurations. In section 3 we perform test cases on two-dimensional NACA23012 and NACA00012 airfoils and simulate glaze icing on the ONERA-M6 and the GLC-305 swept wings.

Numerical method and governing equations

The icing code is developed in the Navier-Stokes solver NSMB (Navier-Stokes-Multi-Block) [START_REF] Vos | Recent Advances in Aerodynamics inside the NSMB (navier-stokes multi-block) Consortium[END_REF]). NSMB is a numerical software solving the Navier-Stokes equations via finite-volume. NSMB is a structured multi-block code parallelized in MPI that includes numerous numerical schemes and turbulence closures (URANS, LES, Hybrid approaches as Detached Eddy Simulation and DDES), several high-order spatial discretization schemes and remeshing techniques. Chimera grids and Immersed boundary method are available as well. The code was initially developed by EPFL (Lausanne, Switzerland), KTH (Stockholm, Sweden), AEROSPA-TIALE, SAAB and CERFACS (Toulouse, France) and is now used and upgraded by CFS Engineering (Lausanne, Switzerland), ICUBE (Strasbourg, France), EPFL, ETH (Zurich, Switzerland), TU (Dresden, Germany), IMFT (Toulouse, France) and Ecole Polytechnique de Montreal (Montreal, Canada).

Methodology

Single step icing is simulated by a one-way procedure. The automatic process requires firstly the computation of the turbulent and compressible flow. The droplet module calculates the droplet field and the impingement efficiency and the liquid water film and ice properties are evaluated by the thermodynamic module. At this point, the ice shape is determined either by the Lagrangian mesh deformation algorithm or the proposed Level-Set method.

Modeling of the Droplet Field

One key element in icing simulation is the evaluation of the amount of water droplets on the exposed body. To model the droplet trajectories and their impingements on a solid structure an Eulerian approach is developed. The Eulerian model solves a droplet velocity field and a volume fraction distribution in the complete computational domain and on the same mesh than the one used by the flow solver. Droplets are considered as spheres, with constant diameter and not subjected to deformation or breaking. The droplet phase is sufficiently diluted to neglect interactions between droplets. Temperature of the droplet phase is set to be constant and equal to free-stream temperature, heat transfer with the surrounding air is neglected. The only forces acting on droplets are the drag, the gravity and the buoyancy. The governing equations for the conservation of mass and momentum of the droplets are written as :

         ∂α ∂t + ∇.(αu) = 0 ∂αu ∂t + ∇. (αu ⊗ u) = α C D Re d 24K (u a -u) + α 1 - ρ a ρ 1 F r 2 g (1)
The variable α represents the non-dimensionalized volume fraction of water and u the non-dimensionalized velocity field of droplets. u a is the nondimensionalized velocity of air, ρ the density of water, ρ a the density of air, g the gravity vector, L is the reference length, K = ρdU ∞ /18Lµ is an inertia parameter, µ a is the dynamic viscosity of air and d is the median diameter of the droplets. The first term on the right-hand-side of the momentum equation accounts for the drag acting on the droplet considered as low Reynolds number spheres. The droplets Reynolds number (Re d ) is defined based on the slip velocity between the air and droplet and the droplet diameter.

Re d = ρd|u a -u| µ a (2)
The drag coefficient behaviour of the droplets C D is given by [START_REF] Schiller | A drag coefficient correlation[END_REF] :

     C D = 0.44 Re d , if Re d > 1000 C D = 24 Re d , if Re d < 0.1 C D = 24×(1+0.15×Re 0.657 d ) Re d (3)
To model the impact of incoming droplets on the body surface, the boundary conditions evolve during the iterative resolution so that Neumann boundary conditions are applied when the incoming droplet flux, defined as the normal product of the droplet velocity vector and the surface vector, is positive and 155 Dirichlet boundary conditions otherwise.

Freestream values of droplet velocity and volume fraction are imposed as far field boundary conditions. The freestream volume fraction is set equal to α ∞ = 1 and the droplet velocity vector equals to air freestream velocities. The ability of the configuration to capture incoming droplets on solid walls is defined by the collection efficiency β. Using the non-dimensionalized formulation for droplet velocities and volume fraction as defined above, the impingement efficiency factor is given by :

β = αu.n (4)
The equations are treated as a generic transport problem on structured grids and recast in a finite volume form. A third order TVD incompressibletype Upwind differencing scheme with a large choice of limiters is used to evaluate the convective fluxes at the interfaces. The source term is split and the negative part treated implicitly allowing a better stability and faster convergence rate of the iterative scheme. A First order backward Euler scheme is used to discretize the time derivative term allowing implicit resolution of the equations. The non-linearity of the convective term is treated with a Picard type approach and the system is solved with a parallel BiGSTAB algorithm.

Thermodynamic Modeling

To calculate the properties of the ice and the liquid water film on the exposed body one need to solve the mass and energy conservations. The input to the equations are the wall shear stress τ wall and the heat transfer coefficient h c evaluated by the flow solver NSMB, the collection efficiency β and the droplet velocity field u computed by the droplet module. The PDE icing model implemented in the NSMB flow solver is based on the Messinger model [START_REF] Messinger | Equilibrium temperature of an unheated icing surface as a function of air speed[END_REF]). The Mass and energy conservation equations are solved for each cell located on the wall. Radiation is neglected in the energy conservation equation as it accounts for a very low contribution in most icing cases.

According to [START_REF] Beaugendre | FENSAP-ICE's threedimensional in-flight ice accretion module: ICE3D[END_REF] the system of equation writes :

         ∂ρ w h ∂t + ∇.(ρ w ūf h) = ṁimp -ṁevap -ṁice ∂ρ w C w hT ∂t + ∇.(ρ w C w ūf hT ) = qimp -qevp -qice -qcnv (5)
h is the water film height, ρ w the water density, C w the specific heat of water, T the surface temperature in Celsius, ṁevp the mass rate of evaporation, ṁimp the impinging mass rate of water droplets and ṁice the resultant mass rate of ice accretion. qimp is the heat rate of impacting droplets, qevp the evaporating heat rate and qice the heat rate of freezing water. The main contribution in the energy balance comes from the heat convection term qcnv . The impinging mass rate contribution is written as ṁimp = LW C. V ∞ . β and the evaporating mass rates ṁevp are determined by a parametric model developed by [START_REF] Kreith | Principles of heat transfer[END_REF] and used in LEWICE.

The heat rate contributions are :

           qimp = ṁimp C w T ∞ + |u| 2 2 qevp = ṁevp Le+Ls 2 qice = ṁice (C i T -L f ) qcnv = h c (T -T rec ). ( 6
)
T rec is the adiabatic wall temperature, C i is the specific heat of ice, L e is the latent heat of evaporation, L s is the latent heat of sublimation and L f is the latent heat of fusion.

ūf is defined as the mean velocity of the water film and is obtained by averaging a linear profile of the water film velocity across the water depth. The mean velocity of the water film ūf writes [START_REF] Beaugendre | Development of a second generation in-flight icing simulation code[END_REF]) :

ūf (x) = h 2µ w τ wall (x), (7)
with τ wall the wall shear stress, x the wall surfaces coordinates and µ w is the viscosity of water.

To close the system and insure a physical solution, the following compatibility relations must be verified [START_REF] Beaugendre | Development of a second generation in-flight icing simulation code[END_REF]) :

         h ≥ 0, ṁice ≥ 0, hT ≥ 0, ṁice T ≤ 0. (8)
The first two relations state that the water depth and the mass of ice can be only positive. The relation hT ≥ 0 states that the temperature of the water film can only be equal or greater than the freezing temperature of water and the last relation states that ice can only be formed when the surface temperature is below the water freezing point.

The heat transfer coefficient is a key parameter in icing simulations.

Most icing codes such LEWICE, CANICE, [START_REF] Cao | Numerical simulation of ice accretions on an aircraft wing[END_REF], [START_REF] Sang | Numerical simulation of icing effect and ice accretion on three-dimensional configurations[END_REF]. . . use boundary layer methods for its determination. Solving the Navier-Stokes equations permits to compute the boundary layer properties. Thus, in NSMB, this coefficient is calculated from the heat fluxes q wall by first computing the recovery wall temperature distribution T rec . A second flow calculation is performed by imposing a constant and arbitrary wall temperature T wall . The computed heat flux q wall obtained from the RANS solver allows the calculation of the heat transfer coefficient by the following formula (see [START_REF] Montreuil | ECLIPPS: 1. three-dimensional CFD prediction of the ice accretion[END_REF] for more details about the calculation) :

h c = q wall -q adiabatic T wall -T rec = q wall T wall -T rec . ( 9
)
The convective terms of the Eq. 5 are discretized via finite volume method with a Roe scheme :

φ Roe (i, j) = 1 2 (F(U i ) + F(U j )) .n - 1 2 |J(U i+1/2 ).n|(U j -U i ) ( 10
)
with φ Roe the evaluated Roe fluxes at interfaces, F the fluxes in each wall surface cell, J the Jacobian of the system, U the icing state vector and n the normals vector. At node i the following resulting semi-discrete system is solved :

           Ω i ρ w ∂h i ∂t -S h,i + j =i φ Roe h (i, j) = 0 Ω i ρ w C w ∂h i T i ∂t -S hT,i + j =i φ Roe hT (i, j) = 0 (11)
j stands as the neighbour index of node i, S as the source term of the system of equations and Ω the volume of the cell.

The partial differential icing equations have apparently three unknowns (h, T, ṁice ) but we can reduce these unknowns to two by defining three types of surfaces and solve the system in each node by a trial and error process similarly to the traditional Messinger model. For a dry surface, corresponding both to rime icing and no-icing case, we can set h = 0 as there is no water runback. For a wet surface, corresponding to glaze icing, there are both ice accretion and liquid water runback therefore we can set T = 0. Finally, for a liquid surface, there is no ice accretion and we set ṁice = 0.

In each case, several terms of the system of equation 11 disappear. The type of the surface being an unknown, a trial and error process must be used to solve the equations. Therefore, the type of surface is guessed and the updated state vector U n+1 is calculated. If the compatibility relations are verified the assumption is correct. If the compatibility relations are not satisfied we switch to another type of surface. We repeat the process until the compatibility relations are satisfied and then switch to the next cell. In NSMB, the system of two equations is solved explicitly in time via a Runge-Kutta scheme and is parallelized with MPI.

A Level-Set Approach for Modeling the Ice/Air interface evolution

The usual process to determine the new iced geometry for the next icing step is to consider the ice growths towards the normal direction to the surface and compute the updated wall coordinates. At this point, the grid can either be remeshed entirely or just deformed and smoothed. We develop in this section a novel approach to track the ice air interface evolution by using the Level-Set method on body fitted grids.

The interface is described by the Level-Set function Φ defined in the computational domain as Φ < 0 is the solid, Φ > 0 is the air and Φ = 0 is the interface [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]).

In the present study, the Level-Set function is initialized in the computational domain as follow :

• Φ = 1, the outside zone (air), in the computational domain • Φ = -1, the inside zone (ice), in the boundary conditions halo cells As boundary conditions, Dirichlet boundary conditions are applied at the wall with values Φ BC,wall = -1 and Neumann boundary conditions are applied on the computational domain far field.

The evolution of the interface is known by solving the so-called Level-Set equation [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]) :

∂Φ ∂t + V∇.Φ = 0 (12)
with V the Level-Set velocity field.

Considering that the ice grows in the direction normal to the wall, an icing velocity can be defined at the wall v wall with the following relation :

v wall = ṁice ρ ice .n (13)
To compute the Level-Set velocity field V, the wall velocity v wall is projected perpendicular to the wall in the computational domain by the following relation :

V(i, j, k) = v wall (i , j , k ) (14)
i , j and k being the indices of the closest wall cell to the computational domain cell of indices i,j,k.

As our goal for this study is to demonstrate the feasibility of the approach, a first order scheme in time and space is chosen to discretize the Level-Set equation. Such low order scheme is expected to be dissipative but it will be shown that this approximation is acceptable in the context of a single step icing model as long as the mesh is sufficiently refined in the Level-Set path and the total Level-Set deformation is moderate.

The Level-Set equation is recast in the following conservative form :

∂Φ ∂t + ∇.(VΦ) = Φ.∇V (15)
The convective terms are discretized similarly to section 2.2 by an incompressibletype first order UDS scheme and the time derivative is evaluated by a first order backward Euler scheme. The following matrix is build for each variable φ : (VΦ):

a P + Ω dt φ n+1 P + i=N,S,E,W,B,F a i φ n+1 i = S n + Ω dt φ n (16)
with (a P , a E , a W , a N , a S , a F , a B ) the fluxes at current point P , east, west, north, south, forward and backward locations. The variable φ is defined as φ : (V x Φ, V y Φ, V z Φ), S is the right hand side and Ω the volume of the cell.

The right hand side part is treated explicitly breaking the unconditional stability of the implicit scheme. It will be shown that this method is still able to handle large timesteps at the cost of higher numerical dissipation. The system of three equations is solved with preconditioned BiGSTAB. The whole process is parallelized with MPI. At the end of the single step calculation, the Level-Set interface, i.e. the contour zero Φ = 0 is extracted by the post processing software Tecplot R . If a multilayer icing calculation is desired, the output contour can be directly used as geometry input for a mesh generator program. The main objective of the framework currently developed in NSMB is to avoid remeshing by using penalization method and wall functions in such case the wall is modelled by a forcing term inside the flow and droplet equations. In such situation the following Level-Set property is useful to track the interface through time :

n = ∇Φ |∇Φ| . ( 17
)

Results and Discussions

In this section, single step icing simulations are performed on two-dimensional NACA23012, NACA0012 aifoils and on the three-dimensional ONERA M6 and GLC-305 swept wings configurations. The flow field is solved implicitly using LU-SGS method. Spatial discretization is achieved using second order central scheme with artificial dissipation or third order Roe Scheme. Spalart-Allmaras [START_REF] Spalart | A One-Equation Turbulence Model for Aerodynamic Flows[END_REF]) and its ONERA roughness extension [START_REF] Aupoix | Extensions of the spalart-allmaras turbulence model to account for wall roughness[END_REF]) is used as turbulence model. Table 1 summarizes the icing cases presented here.

Two-Dimensional cases

Case 1 : NACA23012

Simulation of rime icing is performed on a NACA23012 airfoil. Three body fitted o-grid meshes have been generated. The coarse mesh contains 128 × 128 cells, the medium mesh contains 512 × 256 cells and the fine mesh contains 512 cells. The airfoil chord is c = 1.84m and the computational domain is extended to 50 × c. The flow is simulated at an angle of attack of 2 • , a Mach number of M = 0.32, a Reynolds number of Re = 16.23 × 10 6 , a static temperature of T = 252.45K and with a constant equivalent sand grain number of k s = 0.78mm. The parameters of the icing simulation are a constant and mono-dispersed droplet diameter of d = 15µm, a liquid water content of LW C = 0.30g/m 3 and a total exposure time of t = 1200s. Simulation is compared with results from [START_REF] Broeren | Effect of High-Fidelity Ice-Accretion Simulations on Full-Scale Airfoil performance[END_REF] obtained experimentally in icing tunnel.

Simulations are performed on a quad core 3.5Ghz CPU. Convergence of the flow to Resl2 = 10 -6 is achieved in approximately 2 × 450s on the medium grid, the droplet equations to ResL2 = 10 -5 requires less than 150s of CPU time, the thermodynamic module requires only a single step to converge as the configuration is fully rime and the Level-Set computation takes approximately 5s using a timestep of dt = 2s. The collection efficiency and the instantaneous mass of ice per second per square meter are represented on extruded mesh for more visibility in Fig. 2. One can note that in rime configurations, the collection efficiency and the mass of ice distributions are strongly correlated because the impacting droplets freeze immediately. The collection efficiency peak obtained in this configuration is β = 0.38. The icing velocity field V in 2D coordinates is shown in Fig. 3. One must note that while the icing velocity field is computed in the entire domain in the present study, the velocity field is only necessary in a narrow band from the initial Level-Set location to the expected icing limit location. The Level-Set function after a calculation time of t = 1200s on the medium grid is shown in Fig. 4.

The contour Φ = 0 corresponding to the ice/air interface is plotted in Fig. 5 and compared with the Lagrangian mesh deformation and experimental data for the three meshes. One can show that less accreted mass is obtained on the coarse mesh and noticeable discrepancies are obtained between the deformed grid points and the level-set contour. The medium and fine meshes shows a good agreement with the experimental data and between the Level-Set contour and the deformed grid points. By zooming it can be observed that the Level-Set and lagrangian mesh curves only differs at a maximum of 10 -4 × c at the leading edge. In this rime configuration, one can see that the single step approximation is quite acceptable to reproduce a realistic ice shape in terms of amount and localisation of ice.

The Level-Set contour at t = 1200s is compared with five Level-Set timesteps of dt = (0.1s; 2s; 120s; 240s; 1200s) in Fig. 6. The Level-Set contour computed with timesteps of dt = 0.1s and dt = 2s are in a excellent agreement. Visible discrepancies are obtained with dt = 120s and numerical dissipation increases for larger timesteps.

Case 2 : NACA0012

Glaze icing is simulated on a single block NACA0012 airfoil. The body fitted O-mesh contains 512 × 128 cells. The chord of the airfoil is c = 0.53m and the computational domain extends to 50 × c. A subsonic turbulent flow is simulated at an angle of attack of 4 • , a Mach number of M = 0.32, a Reynolds number of Re = 4.44 × 10 6 , a static temperature of T = 262.04K and with a constant equivalent sand grain roughness number of k s = 0.55mm. The droplet field is simulated with a mono-dispersed droplet diameter of d = 20µm and a liquid water content of LW C = 0.55g/m 3 . The timestep of the Level-Set equation is set to 1s and a total exposure time of t = 420s is calculated.

The heat transfer coefficient against the distance from the stagnation point computed by the flow solver is compared with LEWICE. In despite of the different type of method used by LEWICE and NSMB-ICE a good agreement is obtained. The collection efficiency against the distance from the stagnation point is compared with LEWICE in Fig. 8. The collection efficiency peak is β = 0.72 . In a glaze configuration there is no direct correlation between the collection efficiency and the instantaneous mass of ice per second distributions. As less ice forms near the stagnation point, the remaining water film flows on the surface as runback and freezes downstream mainly on the top of the leading edge. The convergence of the thermodynamic module is plotted in Fig. 9. A L2-Residual drop of ResL2 = 10 -11 is achieved in less than 10s of CPU time. The spurious oscillations of temperature residual during the first steps are due to the change of the type of surface (from dry type to wet type). The evolution of the Level-Set equation is shown in Fig. 10 at exposure times t : (140s, 280s, 420s). One can note that the dissipation errors increase with simulation time and that the results show a distinct mass loss between the Lagrangian mesh deformation and the Level-Set function near the top horn. The visualization of the Level-Set function exhibits dissipation of the solution in this area. This phenomena is expected to be reduced by using mesh refinement in this area and higher order schemes to discretize the Level-Set equations. The level-zero of the Level-Set equation Φ = 0, determined by the post processor visualisation software, is plotted in Fig. 11 and compared with the Lagrangian mesh deformation and the experimental data (Lewice-Cdrom (1996)). It can be noted that a singlestep icing process is only able to reproduce roughly the main top horn. A common observation is that glaze icing simulation requires several icing steps to reproduce a more realistic ice shape [START_REF] Hasanzadeh | Quasi-steady convergence of Multistep Navier-Stokes Icing Simulations[END_REF]).

3.2. Three-Dimensional case 3.2.1. Case 3 : ONERA-M6

The test case is a three-dimensional case performed on the ONERA-M6 swept wing. The mesh contains 884736 cells divided into 4 blocks. The mean aerodynamic chord is c = 0.53m, the semi-span is b = 1m and the computational domain extends to 20 × c. The flow is simulated at an angle of attack of 6 • , a Mach number of M = 0.15, a Reynolds number of Re = 2.17 × 10 6 , a static temperature of T = 265K and with a constant equivalent sand grain roughness number of k s = 0.55mm. The parameters of the icing simulation are a constant and mono-dispersed droplet diameter of d = 20µm, a liquid water content of LW C = 1.0g/m 3 and a total exposure time of t = 500s. The timestep of the Level-Set equation is set to 0.5s. Simulation is performed on four CPU at 3.5Ghz.

As summarized on Table 2, convergence to ResL2 = 10 -5 is achieved in a total CPU time of 16.6 hours for the flow and 3.6 hours for the droplet computations. A total CPU Time of 55min is necessary to achieved a convergence of ResL2 = 10 -12 of the icing equations and a total CPU Time of 6.11 min is necessary for the Level-Set interface to reach the desired exposure time of t = 500s. In the context of a complete icing simulation, the computational time of the Level-Set equation is negligible compared to the calculation times of flow, droplet and thermodynamic as it accounts for only 0.47% of the total time. Several test cases performed previously show the computational times of the Level-Set equations and the mesh deformation (without mesh smoothing) are of the same order of magnitude.

The calculated heat transfer coefficient, the collection efficiency, the water liquid water film height and the instantaneous mass rate of ice per square meter per second distributions are represented in Fig 12 . The heat transfer coefficient distribution is higher near the wing tip inducing more accreted ice in this area and a thinner liquid water film height h. On the other hand, the liquid water film depth h increases at the bottom side of the leading edge near the wing root where the heat transfer coefficient obtained gets lower inducing less accreted ice in this area. These observations are consistent with the fact that the leading edge radius decreases with span location.

Impingement curves are plotted on sections A, B & C respectively at 20%, 60% and 90% of span length and compared with [START_REF] Sang | Numerical simulation of icing effect and ice accretion on three-dimensional configurations[END_REF] in Fig. 13. The present results show a good agreement with those from Sang et al. which were obtained via a Lagrangian droplet tracking approach. The three-dimensional ice interface corresponding to the zero level of Φ is represented at wing tip and wing root in Figs. 14 & 15 at exposure time t = 500s. Ice shape corresponding to Φ = 0 is compared on sections A, B & C with mesh deformation approach in Fig. 16. Ice shape is compared near the wing tip on an alternative cut corresponding to the plane y = 0 in Fig. 17. A global good agreement is obtained between both approaches, the discrepancies are mainly observed on top horn formed near the wing tip. This phenomena is due again to dissipation errors and the coarse mesh impacts the accuracy of the determination of the Level-Zero of the Level-Set function. According to our knowledge there is unfortunately no experimental or numerical data regarding ice shapes to compare on this configuration. Overall, the present results validates the 3D implementation the Level-Set equations.

Case 4 : GLC-305 Swept Wing

The final test case is a simulation of icing on the GLC-305 swept wing. This test case is reported in [START_REF] Papadakis | Experimental investigation of ice accretion effects on a swept wing[END_REF]. Details of the swept wing model and specifications of the sections are provided in Fig. 18. The mean aerodynamic chord of the configuration is MAC = 0.4755m. The mesh contains 786432 cells divided into 4 blocks. The subsonic flow is simulated at an angle of attack of A0A = 4 • , a freestream temperature of T = 267.02K, a Mach Number of Ma = 0.2047, a constant sand-grain roughness of k s = 0.53mm and a Reynolds number (based on the MAC) of Re = 2.432 × 10 6 . The droplet spray is simulated with a mono-dispersed median volume diameter of d = 20µm, a liquid water content of LWC ∞ = 0.65g/m 3 and a total exposure time of t = 600s. A Level-Set timestep of dt = 0.5s is chosen and the calculation is performed on a four-core CPU. The Level-Set interface is compared on three cuts with LEWICE and experimental data on Fig. 19 for cut A, Fig. 20 for cut B and Fig. 21 for cut C. The multilayer ice shapes obtained by LEWICE have been obtained by calculating icing on 2D sections in the streamwise direction and thus neglecting the runback in the transverse direction of the flow. One can note although only single step is simulated, the Level-Set contour calculated by NSMB-ICE fits fairly the experimental data and validate the feasibility of the approach for three dimensional and complex configurations.

Conclusions

A single step icing model using Level-Set methodology has been presented in detail and introduced as a first step in the development of a new generation ice accretion code. This new approach has been implemented in the multiblock structured finite volume compressible flow solver NSMB.

The icing is computed by a one-way procedure on body-fitted multi-block structured grids and parallelized with the MPI library for efficient calculations. Droplets tracking are computed implicitly with a TVD scheme and BiGSTAB algorithm. The PDE icing model is solved explicitly with Roe scheme using a trial and error process. The icing velocity field is projected from the wall to the computational domain and the Level-Set equations are solved implicitly with a first order scheme in time and space.

Rime icing has been simulated on a NACA23012 airfoil and glaze icing has been simulated on a NACA00012 airfoil and on the ONERA-M6 and the GLC-305 swept wings. All results show that the approach is feasible, fast and provides results in good agreement with the traditional Lagrangian mesh deformation technique as long as the mesh is sufficiently refined and the Level-Set deformation is moderate. Higher order schemes will be implemented to discretize the level-set equations in order to reduce the discrepancies observed between both approaches. Local mesh refinement will also be investigated to reduce the numerical dissipation errors. In future work, the promising method introduced in this paper will be extended to an automatic multi-steps approach using penalization method and wall functions. This framework will tackle a recurrent limitation of current ice accretion codes namely the simulation of complex and three-dimensional icing at high exposure time.

16 Case 3 -Comparison of ice shapes obtained with the Level-Set method and the Lagrangian mesh deformation on sections

  & C at total exposure time t = 500s. . . . . . . . . . . . . 17 Case 3 -Comparison of ice shapes obtained with the Level-Set method and the Lagrangian mesh deformation at the tip on plane Y = 0 at total exposure time t = 500s. . . . . . . . . . . 18 Case 4 -Specifications of the GLC-305 swept wing Papadakis

Figure 1 :

 1 Figure 1: Typical forms of icing simulated in NSMB. On top : rime icing on a NFL0414 airfoil at Re = 6.95 × 10 6 , T = 257.59K. On bottom : glaze icing on a SA13112 airfoil at Re = 3.91 × 10 6 , T = 263.19K.

Figure 2 :

 2 Figure 2: Case 1 -Collection Efficiency β and mass of ice per second per square meter ṁice on the medium grid.

Figure 3 :

 3 Figure 3: Case 1 -Visualisation of the icing velocities on the medium grid.

Figure 4 :

 4 Figure 4: Case 1 -Visualisation of the LevelSet function at total exposure time of t = 1200s on the medium grid.

Figure 5 :

 5 Figure 5: Case 1 -Comparison of the Level-Set contour Φ = 0, the deformed mesh and the experimental data at total exposure time t = 1200s on the coarse, medium and fine grids with a fixed Level-Set timestep of dt = 2s.

Figure 6 :

 6 Figure 6: Case 1 -Effect of the Level-Set timestep on the the Level-Set contour Φ = 0 at total exposure time of t = 1200s on the medium grid.

Figure 7 :

 7 Figure 7: Case 2 -Comparison the heat transfer coefficient with LEWICE.

Figure 9 :

 9 Figure 9: Case 2 -Convergence of the icing equations with Roe Scheme.

Figure 10 :

 10 Figure 10: Case 2 -Visualisation of Level-Set function at time t = 140s, t = 280s and t = 420s. The black line is the zero level of the Level-Set function ( i.e. the ice/air interface).

Figure 11 :

 11 Figure 11: Case 2 -Comparison of the Level-Set contour Φ = 0, the Lagrangian mesh and experimental data (Lewice-Cdrom (1996)) at total exposure time t = 420.

Figure 12 :

 12 Figure 12: Case 3 -Visualisation of the heat transfer coefficient HT C, the collection efficiency β, the height of the water film h and the mass of ice per second per square meter ṁice .

Figure 14 :

 14 Figure 14: Case 3 -Visualisation from the wing tip of the ice shape obtained with the Level-Set method at exposure time t = 500s.

Figure 15 :

 15 Figure 15: Case 3 -Visualisation from the wing root of the ice shape obtained with the Level-Set method at total exposure time t = 500s.

Figure 16 :

 16 Figure 16: Case 3 -Comparison of ice shapes obtained with the Level-Set method and the Lagrangian mesh deformation on sections A,B & C at total exposure time t = 500s.

Figure 17 :

 17 Figure 17: Case 3 -Comparison of ice shapes obtained with the Level-Set method and the Lagrangian mesh deformation at the tip on plane Y = 0 at total exposure time t = 500s.

Case 3 -Visualisation from the wing root of the ice shape obtained with the Level-Set method at total exposure time t = 500s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

et al. (2005). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 Case 4 -Comparison of the Φ = 0 contour with LEWICE and experimental data on section A of the GLC-305 swept wing. . 20 Case 4 -Comparison of the Φ = 0 contour with LEWICE and experimental data on section B of the GLC-305 swept wing. . 600 21 Case 4 -Comparison of the Φ = 0 contour with LEWICE and experimental data on section C of the GLC-305 swept wing. .

Acknowledgments

The authors would like to thank HPC Strasbourg and CINES for computational resources, Kazem Hasanzadeh (Polytechnique Montreal) and Dominique Charbonnier (CFS-Engineering) for grid generation. The partial financial support from NSERC is acknowledged.

List of Figures

Typical forms of icing simulated in NSMB. On top : rime icing on a NFL0414 airfoil at Re = 6.95 × 10 6 , T = 257.59K. On