domain: phys spi

text:
Pierre Lavoie
email: pierre.lavoie@polymtl.ca

Emmanuel Radenac

Ghislain Blanchard

Eric Laurendeau

Philippe Villedieu

Phd Candidate

Immersed boundary methodology for multistep ice accretion using a level set

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

I. Introduction

Numerical tools for the prediction of in-flight ice accretion are typically based on a quasi-steady assumption where modules are called sequentially and solved to steady state within a time-iterative scheme (e.g., LEWICE [START_REF] Wright | User Manual for the NASA Glenn Ice Accretion Code LEWICE Version 2.2.2[END_REF], FENSAP-ICE [START_REF] Beaugendre | Development of a Second Generation In-Flight Icing Simulation Code[END_REF]). The process is illustrated in Fig. 1 where the modules are: (1) a mesh generation tool, (2) a solver for the aerodynamics, (3) a solver to obtain the droplet trajectories and impingement rates, (4) a solver to obtain the wall convective heat transfer (in the boundary layer), [START_REF] Wright | Validation Methods and Results for a Two-Dimensional Ice Accretion Code[END_REF] a solver to perform a heat and mass balance applied to the deposited water to obtain the ice accretion rate and finally (6) a tool to update the geometry based on the ice thickness evolution. Modules [START_REF] Wright | User Manual for the NASA Glenn Ice Accretion Code LEWICE Version 2.2.2[END_REF] to [START_REF] Petrosino | Ice Accretion Model on Multi-Element Airfoil[END_REF] are embedded in a time loop for which the total ice accretion time is divided in time steps (multi-step) generating successive layers of ice (multi-layer). When using Body-Fitted (BF) meshes, a mesh update is required with each new ice layer generated by the multi-step process, and it can be repeated several times in order to obtain the final ice shape prediction. This leads to additional costs related to the mesh update and additional difficulty in updating the ice shape which can exhibit unphysical surface overlaps in concave regions when using a Lagrangian approach (displacement of surface mesh nodes, a method which is usually employed in ice accretion codes [START_REF] Saeed | Modified Canice for Improved Prediction of Airfoil Ice Accretion[END_REF][START_REF] Bourgault-Côté | Multi-Layer Icing Methodologies for Conservative Ice Growth[END_REF]). Tools for the numerical prediction of ice accretion are constantly evolving. Originally, icing tools were limited to 2D simulations and mostly used potential flow solvers (e.g., panel method) coupled with an integral boundary layer code (e.g., LEWICE [START_REF] Wright | User Manual for the NASA Glenn Ice Accretion Code LEWICE Version 2.2.2[END_REF][START_REF] Wright | Validation Methods and Results for a Two-Dimensional Ice Accretion Code[END_REF] and CANICE [START_REF] Saeed | Modified Canice for Improved Prediction of Airfoil Ice Accretion[END_REF]). While these tools offer fast computation times, they are not well suited for configurations involving separated flows and recirculations (e.g., multi-element airfoil [START_REF] Petrosino | Ice Accretion Model on Multi-Element Airfoil[END_REF]). They are also not easy to generalize for 3D applications and more complex configurations (e.g., full aircraft). Nowadays, most icing software use either a Euler [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF] or Reynolds Averaged Navier-Stokes (RANS) flow solver [START_REF] Beaugendre | Development of a Second Generation In-Flight Icing Simulation Code[END_REF][START_REF] Bourgault-Côté | Multi-Layer Icing Methodologies for Conservative Ice Growth[END_REF][START_REF] Gori | PoliMIce: A simulation framework for threedimensional ice accretion[END_REF][START_REF] Radenac | Validation of a 3D ice accretion tool on swept wings of the SUNSET2 program[END_REF][START_REF] Pena | A single step ice accretion model using Level-Set method[END_REF]. When using a Euler flow solver, it must be coupled with a boundary layer code in order to model the viscous effects and retrieve the convective heat transfer at the wall. This task is not trivial for 3D simulations but is the subject of recent research activities at the ONERA [START_REF] Radenac | Use of a Two-Dimensional Finite Volume Integral Boundary-Layer Method for Ice-Accretion Calculations[END_REF], where a PDE-based formulation is developed for the purpose of icing applications. The modeling of surface roughness is also important as it directly influences the convective heat transfer at the wall. The roughness is often modeled as a constant over the entire wall, but new models with spatially varying roughness heights have been developed recently [START_REF] Fortin | Heat and mass transfer during ice accretion on aircraft wings with an improved roughness model[END_REF][START_REF] Mcclain | Ice accretion roughness measurements and modeling[END_REF][START_REF] Han | Surface Roughness and Heat Transfer Improved Predictions for Aircraft Ice-Accretion Modeling[END_REF]. These models can potentially improve the prediction of glaze ice shapes. Concerning the evaluation of the droplet trajectories and impingement rates, Lagrangian methods were originally used and are still common in current icing tools [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF][START_REF] Gori | PoliMIce: A simulation framework for threedimensional ice accretion[END_REF][START_REF] Bidwell | Icing Analysis of a Swept NACA0012 Wing Using LEWICE3D Version 3.48[END_REF][START_REF] Bellosta | A robust 3D particle tracking solver for in-flight ice accretion using arbitrary precision arithmetic[END_REF]. On the other hand, Eulerian formulations for the droplets have gained interest in the recent years [START_REF] Bourgault-Côté | Multi-Layer Icing Methodologies for Conservative Ice Growth[END_REF][START_REF] Pena | A single step ice accretion model using Level-Set method[END_REF][START_REF] Bourgault | A finite element method study of Eulerian droplets impingement models[END_REF][START_REF] Radenac | IGLOO3D Computations of the Ice Accretion on Swept-Wings of the SUNSET2 Database[END_REF] as they are easy to generalize for complex 3D configurations. Both Lagrangian and Eulerian formulations are also the subject of ongoing research for the modeling of Supercooled Large Droplets (SLD) [START_REF] Honsek | Eulerian Modeling of In-Flight Icing Due to Supercooled Large Droplets[END_REF][START_REF] Trontin | Revisited Model for Supercooled Large Droplet Impact onto a Solid Surface[END_REF] which requires the consideration of additional physics such as droplet breaking, splashing and bouncing. The thermodynamic module is responsible for determining the ice accretion rate (or ice thickness) which is then used to update the geometry (the ice shape). The mass and energy balances usually follow the work of Messinger [START_REF] Messinger | Equilibrium Temperature of an Unheated Icing Surface as a Function of Air Speed[END_REF] but the methods have seen some improvements to be more suitable for 3D simulations or to model more physics.

For instance, iterative algebraic methods (e.g., [START_REF] Zhu | 3D Ice Accretion Simulation for Complex Configuration Basing on Improved Messinger Model[END_REF]) and PDE-based methods (e.g., Shallow Water Icing Model (SWIM) [START_REF] Bourgault | Development of a Shallow-Water Icing Model in FENSAP-ICE[END_REF]) can be used for 3D simulations. Some methods also include more physics such as the conductive heat transfer through the ice thickness (e.g., [START_REF] Myers | Extension to the Messinger Model for Aircraft Icing[END_REF][START_REF] Gori | Local Solution to the Unsteady Stefan Problem for In-Flight Ice Accretion Modeling[END_REF]). These methods assume the ice is continuous, but alternative approaches are also developed where the ice is made of discrete elements and the ice evolution is based on probabilistic arguments.

For instance, the morphogenetic model of [START_REF] Szilder | Comparing Experimental Ice Accretions on a Swept Wing with 3D Morphogenetic Simulations[END_REF] follows this approach and provides an interesting step forward in the prediction of ice scallops and other complex ice shapes specific to 3D simulations. Although significant improvements have been made over the years, some ice shapes are still very difficult to predict even to this day (e.g., glaze ice on swept wings [START_REF] Fujiwara | Comparison of Computational and Experimental Ice Accretions of Large Swept Wings[END_REF]). Thus, the numerical prediction of ice accretion still remains a very active field of research. In addition to improving individual modules constituting the ice accretion process, research must also address the automation and robustness of the framework as a whole. The robustness of ice accretion tools is often limited by the difficulty of generating meshes on complex ice shapes and also by the geometry update which can exhibit overlaps in concave regions if not treated properly. Instead of performing a surface evolution and mesh deformation (e.g., [START_REF] Tong | Three-Dimensional Surface Evolution and Mesh Deformation for Aircraft Icing Applications[END_REF]), another type of strategy is investigated : the use of an Immersed Boundary Method (IBM) combined with a level-set method. This methodology has the potential to alleviate the issues related to the mesh and geometry update. The objective of this paper is thus twofold: confirm the potential of this methodology and assess its accuracy against the usual body-fitted approach.

Immersed Boundary Methods have been applied with success in many fields but have rare applications to the prediction of ice accretion. The CIRA used a discrete IBM to solve compressible inviscid flows on 3D Cartesian grids in [START_REF] Capizzano | A Compressible Flow Simulation System Based on Cartesian Grids with Anisotropic Refinements[END_REF] which was later extended for compressible viscous flows [START_REF] Capizzano | Turbulent Wall Model for Immersed Boundary Methods[END_REF][START_REF] Capizzano | Coupling a Wall Diffusion Model with an Immersed Boundary Technique[END_REF]. A discrete IBM was also applied to an Eulerian droplet solver [START_REF] Capizzano | A Eulerian Method for Water Droplet Impingement by Means of an Immersed Boundary Technique[END_REF] with the intent of performing ice accretion simulations. However, no ice prediction results using these IBMs have been published yet.

Another research team from the University of Strasbourg applied IBMs to a 3D ice accretion code (NSMB-ICE) with a level-set approach [START_REF] Al-Kebsi | Multi-Step Ice Accretion Simulation Using the Level-Set Method[END_REF]. The compressible viscous flow solver employs a penalty method (an IBM) and the Eulerian droplet solver uses a discrete approach similar to the one of [START_REF] Capizzano | A Eulerian Method for Water Droplet Impingement by Means of an Immersed Boundary Technique[END_REF]. The level-set approach initially proposed in [START_REF] Pena | A single step ice accretion model using Level-Set method[END_REF] is used to update the iced geometry, but is applied to multi-step ice accretion (up to 5 steps on a single test case). According to the authors, the implementation is currently limited to laminar flow and rime ice. Furthermore, no comparison is made against a more classical body-fitted approach nor experimental results. The idea of using a level-set for the prediction of ice accretion is also reused in [START_REF] Bourgault-Côté | Multilayer Airfoil Ice Accretion Simulations Using a Level-Set Method with B-Spline Representation[END_REF] where multi-step simulations are performed in 2D combined with the use of NURBS. An explicit tracking of the air/ice interface is also performed in order to enforce mass conservation.

Our initial attempt at applying IBMs to an ice accretion suite is presented in [START_REF] Lavoie | A Penalization Method for 2D Ice Accretion Simulations[END_REF], where a penalization method is applied to the aerodynamic (Euler equations) and droplet solvers (Eulerian formulation). The 1 st step of the multi-step process is performed using a Body-Fitted mesh while for the subsequent steps, the ice shape is immersed on the initial mesh. A geometric approach was used to evaluate the signed distance field required by the penalization method. Furthermore, a Lagrangian node displacement approach was used to update the geometry. This paper extends the contribution of [START_REF] Lavoie | A Penalization Method for 2D Ice Accretion Simulations[END_REF] with two key features. First, an improved penalization method suitable for ice horn accretion is applied to the Euler equations [START_REF] Lavoie | An Improved Characteristic Based Volume Penalization Method for the Euler Equations Towards Icing Applications[END_REF]. Second, a level-set approach [START_REF] Pena | A single step ice accretion model using Level-Set method[END_REF] is implemented in the IBM multi-step ice accretion process to solve the issues related to unphysical geometry update, replacing the Lagrangian geometry update.

The paper first presents the ice accretion suite used as the development platform, IGLOO2D [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF]. Then, implementation details for the penalization and level-set methods are covered: preprocessing, penalization of the Euler equations, penalization of the droplet solver, extraction of the surface data and implementation of the level-set method. A section discusses the benefits of using the level-set method in the IBM ice accretion framework on a manufactured case. Then, rime and glaze ice cases from the AIAA Ice Prediction Workshop (IPW) [START_REF] Broeren | 1st AIAA Ice Prediction Workshop[END_REF] are used for verification where the Body-Fitted and penalized solutions are compared using the multi-step process. Additional verifications are performed for the ice accretion around a NACA0012 airfoil. Finally, ice accretion simulations are performed on a three-element airfoil before conclusions are drawn.

II. Methodology

The 2D ice accretion suite IGLOO2D [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF] is used as the development environment. In IGLOO2D, different types of solvers are available for each module but only the ones used in this paper are discussed. The unstructured mesh generation is handled by GMSH [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The aerodynamic field is evaluated using the Euler equations and the convective heat transfer is evaluated using a simplified integral boundary layer method (SIM2D) [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF]. For the droplet trajectories and impingement evaluation, the Eulerian solver is selected. The ice accretion solver is based on a Messinger-type mass and energy balance to obtain the ice thickness. Finally, the ice geometry is generated by a Lagrangian displacement of the surface nodes. In this paper, a level-set method is also integrated in the ice accretion framework which corresponds to an Eulerian displacement of the geometry.

The modules can be classified either as volume or surface solvers. The aerodynamics (EULER2D) and the Eulerian droplet trajectories (TRAJE2D) are solved on 2D volume meshes. On the other hand, the simplified integral boundary layer method (SIM2D) and the ice accretion (MESSINGER2D) are solved on 1D surface grids.

For the application of the IBM, the suggested approach is to start the multi-step ice accretion process from a standard BF mesh, thus keeping the original BF solution for the 1 st ice layer (as well as for the clean areas of the surface for the following steps). Usually, the BF mesh is updated to match the new ice geometry for each subsequent step. With our IBM, the volume mesh update is avoided and a penalization method is applied to the volume solvers (airflow and droplets trajectory) to impose the correct boundary conditions on the Immersed Boundary (IB) which arbitrarily cuts through the mesh. The surface mesh representing the ice shape (the IB) is however re-meshed in order to retain an adequate discretization of the ice shape (e.g. more refined near high-curvature features). The use of the penalization method requires some modifications to the ice accretion suite : the addition of a preprocessing step, the modification of the volume solvers and the extraction of surface data, as highlighted in red in Fig. 2. These modifications are discussed in the following sections along with the integration of the level-set method in the multi-step process.

A. Immersed Boundary Pre-Processing

In this paper, both an explicit and implicit definition of the IB are required. The multi-step ice accretion process starts from a BF volume mesh, thus initially providing a surface mesh which represents the solid-air interface. It can be interpreted as an Immersed Boundary which correspond to the BF surface for the first step. This explicit definition must be conserved throughout the multi-step ice accretion process in order to use the surface solvers (i.e., ice accretion, boundary layer). On the other hand, the penalization methods (a type of IBM) implemented in IGLOO2D use a signed distance field (implicit definition) to obtain information about the interface at any point in the volume mesh.

The IB preprocessor evaluates the signed distance field (φ) by first detecting the inside (solid) and outside (fluid) cells. Knowing the list of edges defining a closed IB, a ray casting algorithm can be used for this matter [START_REF] Schneider | Chapter 13 -computational geometry topics[END_REF]. Once this information is known, it can be used to determine the sign of the signed distance field where φ > 0 in the fluid and φ < 0 in the solid. The distance is evaluated by taking advantage of the available explicit definition of the interface.

For each cell center, a geometric approach determines the minimum projected distance to the list of edges (or faces in 3D) defining the IB. It comes back to evaluating the minimum distance between a point and a segment in 2D [START_REF] Schneider | Chapter 13 -computational geometry topics[END_REF] or between a point and a triangle in 3D. Then, the normals to the IB (n φ ) and its curvature (κ) can be evaluated from the signed distance field (φ) as:

n φ = - ∇φ ||∇φ|| (1)
κ = ∇ • n φ (2)
Notice that n φ is defined to point towards the solid, contrary to the usual definition, which is useful in the implementation of the penalization methods. An example of signed distance field around a clean NACA23012 airfoil is illustrated in Fig. 3 along with the normals to the wall (n φ ). Here, the signed distance field is strictly positive because it is evaluated on a BF mesh where the contour φ = 0 is the surface of the airfoil.

Fig. 3 Signed distance contours (φ) and surface normals (n φ ) for a clean NACA23012 airfoil

Although the volume mesh update is avoided when using the IBM, the IB is re-meshed at the pre-processing phase for the 2 nd ice accretion step and further. This is possible because the IB discretization (surface mesh) is independent of the volume mesh. The surface re-meshing is done using GMSH where a B-spline is fitted through the discrete list of nodes defining the ice shape. The nodes are then redistributed according to a user-specified characteristic mesh size.

This provides a surface mesh discretization which is very close to what is obtained with the BF ice accretion process.

B. Volume Penalization Method

With a volume penalization method, the boundary conditions are applied by the addition of source terms in the continuous form of the governing equations in order to enforce the desired condition at the IB. The source terms are turned on if the computational volume (a cell) is located inside the solid zone and turned off if in the fluid zone. Hence the governing equations are solved as usual in the fluid but penalized in the solid. The source terms are turned on and off using a mask function (χ) which takes the form of a sharp Heaviside function (Eq. ( 3)).

χ =        0 φ ≥ 0 (fluid) 1 φ < 0 (solid) (3)
For the aerodynamics, the penalization of the Euler equations is performed using the CBVP-Hs method of [START_REF] Lavoie | An Improved Characteristic Based Volume Penalization Method for the Euler Equations Towards Icing Applications[END_REF]. This approach enforces the no-penetration velocity (slip wall, Eq. ( 4)) and uses the normal momentum relation to account

for the wall curvature in the pressure extrapolation (Eq. ( 5)). The conservation of total enthalpy (Eq. ( 6)) and entropy (Eq. ( 7)) are also enforced across the IB to close the system.

(v • n φ )n φ = 0 (4)
n φ • ∇P = κρ||v|| 2 (5)
n φ • ∇H = 0 (6)
n φ • ∇s = 0 (7)
This method was found to provide better flow attachment and a better conservation of entropy and total enthalpy for ice shapes exhibiting high curvature such as ice horns [START_REF] Lavoie | An Improved Characteristic Based Volume Penalization Method for the Euler Equations Towards Icing Applications[END_REF]. The set of penalized Euler equations is given by Eq. ( 8) where the penalization terms enforcing Eqs. ( 4)-( 7) are gathered on the right-hand side (RHS).

∂ρ ∂t + (1 -χ) ∇ • (ρv) = - χ η c n φ • ∇ρ -κ ρ 2 γP ||v|| 2 ∂(ρv) ∂t + (1 -χ) ∇ • (ρv ⊗ v + P I) = - χ η c n φ • ∇(ρv) + κρv 1 - ρ γP ||v|| 2 - χ η ρ(v • n φ )n φ ∂(ρE) ∂t + (1 -χ) ∇ • ((ρE + P )v) = - χ η c ρn φ • ∇H - χ η ρ(v • n φ ) 2 (8)
In Eq. ( 8), v is the air velocity, I is the identity tensor, E is the total energy, η and 1/η c are penalization parameters that can be respectively interpreted as a characteristic time and a characteristic velocity. For more details about the derivation of Eq. ( 8), the interested reader is invited to read reference [START_REF] Lavoie | An Improved Characteristic Based Volume Penalization Method for the Euler Equations Towards Icing Applications[END_REF].

For the Eulerian droplet equations, the penalization method of [START_REF] Lavoie | A Penalization Method for Eulerian Droplet Impingement Simulations towards Icing Applications[END_REF] is used. When droplets impinge the body

(v d • n φ > 0)
, no penalization is applied and the physical equations are solved in the solid zone (φ < 0). The droplets are thus allowed to cross the IB and enter the body. However, when the droplets enter the computational domain from the solid zone (v d • n φ ≤ 0) a boundary condition is applied on the primitive variables (Eq. ( 9)), enforcing a null flux and avoiding re-injection of the droplets.

α = 0 v d = 0        if v d • n φ ≤ 0 (9)
To translate this behavior to the droplet equations using penalization terms, the usual mask function (χ, Eq. ( 3)) is multiplied by a droplet mask function (χ d , Eq. ( 10)), ensuring the penalization term is only active in the solid zone if the droplets are reinjected in the fluid zone.

χ d =        0 αv d • n φ ≥ 0 (impingement) 1 αv d • n φ < 0 (re-injection) (10)
The set of penalized droplet equations is given by Eq. ( 11), where the influence of gravity is neglected and the penalization terms are highlighted in red.

∂α ∂t + ∇ • (αv d ) = - χχ d η α ∂(αv d ) ∂t + ∇ • (αv d ⊗ v d ) = C D Re d 24Stk α(v a -v d )-2 χχ d η αv d (11)
In Eq. ( 11), α is the non-dimensional volume fraction of water, v d is the non-dimensional droplets velocity, v a is the non-dimensional air velocity and C D is the droplets drag coefficient. The droplets Reynolds number (Re d ) and the Stokes number (Stk ) are:

Re d = ρ a ||v a -v d ||D d µ (12)
Stk = ρ d D 2 d U ∞ 18Lµ (13)
where D d is the droplet diameter, µ the dynamic viscosity of air and L a characteristic dimension (e.g., the chord length for an airfoil). The drag model of Schiller and Naumann [START_REF] Schiller | A drag coefficient correlation[END_REF] is used for the droplets which are assumed to remain spherical:

C D =        24 Re d (1 + 0.15Re 0.687 d ) Re d ≤ 1000 0.4 Re d > 1000 (14)
The penalization parameters must be small (η 1, η c 1) to accurately enforce the boundary conditions, which leads to a stiff system of equations. The penalization terms are thus treated implicitly when solving the system of equations for both the aerodynamics and droplet trajectories.

C. Surface Data Extraction

Relevant surface information from the volume solvers (aerodynamics and droplet trajectory) must be communicated to the surface solvers (boundary layer, ice accretion) at each step of the multi-step loop (e.g., pressure, velocity, droplet impingement rate). However, the penalization method does not explicitly provide the data on the Immersed Boundary (IB). Instead, the variables are known in the surrounding cells and an additional extraction step is thus required to recover the surface data.

In this paper, the data is interpolated at the node defining the IB (the surface points) using a weighted least square approach. The nearest cell to the interpolation point is first detected. Then, all the cells sharing a node with the identified cell are flagged as neighbors and included in the interpolation stencil. The penalization methods used in this paper are designed to fill the solid cell with valid data. The solid cells are included in the interpolation stencil, hence the need for methods ensuring a controlled continuity of the solution across the solid-fluid interface, as described in [START_REF] Lavoie | An Improved Characteristic Based Volume Penalization Method for the Euler Equations Towards Icing Applications[END_REF][START_REF] Lavoie | A Penalization Method for Eulerian Droplet Impingement Simulations towards Icing Applications[END_REF].

The interpolation uses an inverse distance weight with a smoothing parameter to avoid dividing by zero when the interpolation point and stencil points are too close. The weight between a cell center J (part of the stencil) and the interpolation point P is evaluated as:

w J = 1 ||r P J || 2 + 2 (15)
where ||r P J || is the distance between P and J. The smoothing parameter is selected as = 0.5∆x J with ∆x J the characteristic size of cell J,

D. Geometry Update via the Level-Set method

A Lagrangian approach can be used to update the geometry according to the normals to the wall (n, pointing towards the fluid) and the ice thickness (h ice ) provided by the ice accretion solver. A simple node update can be performed as:

x new = x old + h ice n (16)
where x new and x old are respectively the new and old node locations. This type of approach does not naturally handle the overlaps that can occur near concave region and requires methods for collision detection and front merging to obtain a usable surface mesh. A simple fix can be implemented in 2D as described in [START_REF] Ruff | Users Manual for the NASA Lewis Ice Accretion Prediction Code (LEWICE)[END_REF]. However, it does not directly translate to a 3D implementation which involves more complex geometric operations on a 2D surface mesh.

Alternatively, the level-set method can be used to update the geometry. This was done for instance in [START_REF] Pena | A single step ice accretion model using Level-Set method[END_REF] where the level-set equation (Eq. ( 17), [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF]) is used with an icing velocity field (V ice ) and solved on the volume mesh.

∂φ ∂t + V ice • ∇φ = 0 (17)
This approach has the benefit of being valid for both 2D and 3D simulations. It also naturally handles the issues related to the geometry update such as geometry overlaps. Here, the level-set method reuses the signed distance field (φ) computed at the IB pre-processing step. The interface (IB or BF) is represented by the contour φ = 0 and is advanced in time (Eq. ( 17)) to generate the ice shape, following the icing velocity vector field V ice . In this paper, the level-set is discretized using a 2 nd order scheme in time (Heun's method) and space (upwind with MUSCL extrapolation). The following sections describe a method to retrieve the icing velocity field and discuss the need for a re-initialization step in the advection of the level-set.

Velocity Propagation

The icing velocity magnitude (V ice,surf ) can be computed from the ice accretion time (∆t ice ) and the ice thickness (h ice ) provided on the surface mesh by the thermodynamics solver.

V ice,surf = h ice ∆t ice (18)
However, V ice,surf must be propagated in the volume mesh in order to perform the level-set advection (Eq. ( 17)). To obtain a behavior similar to the Lagrangian node displacement approach (Eq. ( 16)), the icing velocity is propagated from the surface mesh in the normal direction, producing constant velocity bands. A PDE-based approach (Eq. ( 19)) is used to propagate the information from the surface (V ice,surf ) to the field (V ice ) following the normal direction to the surface.

∂V ice ∂t = sign(φ)n φ • ∇V ice (19)
When the surface mesh corresponds to the body-fitted mesh boundary, V ice,surf is imposed as a Dirichlet boundary condition using ghost cells and is propagated in the fluid zone to obtain (V ice ). When the ice shape is immersed in the mesh (IB), the surface no longer corresponds to the mesh boundaries. For instance, this occurs from the 2 nd ice layer onward in the multi-step icing process. In this situation, a band of cells in the vicinity of the interface is initialized by a nearest neighbor search, taking advantage of the explicit definition of the interface. These cells are then frozen (no update) so they can act as ghost cells when solving Eq. ( 19) on both sides of the IB. The update is prevented by setting the Right Hand Side (RHS) of Eq. ( 19) to zero for the frozen cells. The propagation Eq. ( 19) accounts for the sign of φ in order to propagate V ice,surf from the band of initialized cells towards the fluid (φ > 0) and solid zones (φ < 0).

Once the icing velocity magnitude is known in the volume mesh, the vector field is set as:

V ice = -V ice n φ,0 (20)
where n φ,0 represents the normal to the initial contour φ = 0 (before the advection process begins). In other words, the icing velocity field remains fixed during the advection of the level-set. An example of propagated icing velocity field is illustrated in Fig. 4, showing the constant velocity bands in the normal direction to the interface.

x

Level-Set Advection and Re-Initialization

While the contour φ = 0 is advected using the level-set equation ( 17), φ does not conserve the properties of a signed distance field [START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF]. A re-initialization of the level-set is thus performed (i.e., the signed distance field is re-evaluated).

This could be done by reusing the geometric approach from the pre-processing step. Because the new location of the IB is only known via its implicit definition at this stage, this would imply the application of a contour extraction technique to obtain an explicit definition of the interface (new surface mesh). The signed distance field is instead updated using the re-initialization equation [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF], as follows:

∂φ ∂t = S(φ 0 ) (n φ • ∇φ + 1) (21)
S(φ 0 ) = φ 0 φ 2 0 + 2 (22)
This equation incorporates a smoothed sign function S(φ 0 ) which is based on the signed distance before re-initialization (φ 0 ). According to [START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF], it ensures that φ remains unchanged at the interface during the re-initialization process. In practice, numerical experiments showed the introduction of wiggles in the contour φ = 0 when using this approach, an undesirable behavior as a surface mesh is to be constructed from this extracted interface. To ensure the interface remains exactly at the same location, the idea used for the velocity propagation is repurposed here: freezing the update of a band of cells in the vicinity of the interface. Again, it is done by setting the RHS to zero in Eq. ( 21) for the frozen cells. In this way, the band of cells is not updated and the signed distance remains φ 0 . This approach follows the assumptions that φ remains close to a signed distance field in the vicinity of the interface. In this paper, two iterations of Eq. ( 21) are performed at every time step of the level-set advection process (Eq. [START_REF] Bourgault | A finite element method study of Eulerian droplets impingement models[END_REF]). An example of level-set advection is shown on Fig. 5, where the φ contours are displayed inside the ice shape only. Without re-initialization (Fig. 5a), the signed distance field is distorted inside the solid while activating the re-initialization (Fig. 5b) provides a more regular and sensible solution.

E. Surface Mesh Extraction

Once the level-set advection is completed and the signed distance field is re-initialized, the surface mesh extraction can be performed. It consists of two parts: (1) the contour extraction providing an explicit definition of the IB from the level-set and (2) the meshing of the surface (i.e., using GMSH). The first part is performed in the level-set module (geometry evolution solver) while the second part is performed when preprocessing the IB.

A surface discretization can be obtained by performing the extraction of the contour φ = 0. Note that the contour extraction is not performed before the re-initialization because of our general method which combines a Body-Fitted surface for the clean geometry and an Immersed Boundary for the ice shape. It is desired to extract a contour which matches the BF surface where there is no ice while the IB is normally extracted for the iced zones. To do so, a check based on φ is performed at the BF surface during the extraction process. If the BF surface is in the solid zone and far from the contour φ = 0 (e.g., leading edge in Fig. 5a), the signed distance field can be so distorted that the contour may be falsely detected at the BF surface. The re-initialization of the level-set solves this issue and thus helps the contour extraction.

For the contour extraction, well-known methods such as the marching cubes [START_REF] Lorensen | Marching Cubes: A High Resolution 3D Surface Construction Algorithm[END_REF], marching tetrahedra [START_REF] Treece | Regularised marching tetrahedra: improved iso-surface extraction[END_REF] and marching squares (a 2D equivalent of the marching cubes) can be used. In this paper, 2D unstructured meshes made of triangles are used which allows for a simple contour extraction method. We consider four possible configurations for the contour intersection with a triangular cell: edge to edge, edge to vertex, vertex to vertex or vertex only (Fig. 6).

Edge to edge Vertex

Edge to vertex These cases are all handled automatically by performing an edge-based interpolation assuming a single intersection point per edge. The process marches from cell to cell and adds consecutive intersection points to a linked list, forming a surface discretization. Tested edges are tagged along the way to avoid adding duplicates to the list. Once the marching process can no longer find any intersection on untested edges, the contour is completed. An edge is intersected by the contour if there is a sign change in φ between its two vertices. If φ at one vertex is below a specified threshold, the intersection is assumed to occur at the vertex and no interpolation is made. In this case, all the edges sharing the vertex are tagged as tested. This approach retains the discretization of the body-fitted surface where there is no ice accretion (φ ≈ 0) and perform a more classical contour extraction for the IB. It also directly provides an ordered list of points (surface mesh) for each body when dealing with a multi-element configuration. Note that in IGLOO2D, φ is reconstructed at the vertices from a weighted least square interpolation using the cell-center solution. An example of the marching process illustrated in Fig. 7 where the vertices are identified as positive, negative or zero and the extracted contour is illustrated in red.

Fig. 7 Example for the contour extraction marching process

The extraction process usually produces an irregular discretization where nodes can be very close to each other when the edge intersection is detected near a vertex. To help retrieve a more uniform and smoother surface mesh, nodes are merged if they are too close and inserted if they are too far apart. The merge and insertion process is performed by an arithmetic average followed by a correction to bring the node back on the level-set. Following the idea presented in [START_REF] Enright | A Hybrid Particle Level Set Method for Improved Interface Capturing[END_REF], the correction takes the form:

x corrected = x merged/inserted + (φ -φ target )ψn φ (23)
where φ target = 0 and ψ is a relaxation parameter set to unity but that can be reduced to avoid erroneous corrections (e.g., point near the wrong contour if multiple contours are involved). An example of contour extraction is provided in Fig. 8a near the leading edge of an iced NACA23012. The effect of node merging and insertion is illustrated in Figs.

8b-8c. Note that this node correction process is not mandatory as the surface is later re-meshed using GMSH. However it helps in avoiding surface wiggles due to an irregular initial node distribution when fitting a spline.

III. Ice Accretion Results

In this section, the new ice accretion framework using the IB and level-set methods is assessed. The objective is to reproduce the ice accretion results obtained with a classical BF approach while improving the robustness of the numerical tool (e.g., no negative volumes in the mesh, no overlaps in the surface mesh). In order to demonstrate the benefits of using the level-set approach, ice accretion over a manufactured ice shape is first performed using the level-set method and compared to the Lagrangian node displacement method. Then, rime ice case 241 and glaze ice case 242 from the 1 st AIAA Ice Prediction Workshop (IPW, [START_REF] Broeren | 1st AIAA Ice Prediction Workshop[END_REF]) are tested. These cases are respectively run ED1977 and run ED1978 from [START_REF] Lee | Implementation and Validation of 3-D Ice Accretion Measurement Methodology[END_REF], with slightly corrected icing conditions. Additional ice accretion cases from [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF] (cases 001, 003 and 004) are also tested to further demonstrate the behavior of the IBM. Finally, the new framework is tested on the multi-element McDonnell-Douglas LB606b Airfoil (MDA) [START_REF] Petrosino | Ice Accretion Model on Multi-Element Airfoil[END_REF] to illustrate the flexibility of the method on complex high-lift systems. The simulation parameters are summarized in Table 1. In this section, two methods are available for the representation of the ice shape (IBM or BF) and two for the geometry update (node displacement or Level-Set). This makes four possible combinations. When simply referring to the Immersed Boundary Method, the use of the level-set method is implied. Similarly, when referring to the Body-Fitted method, the use of the Lagrangian node displacement approach is implied (the standard approach in the icing community). In addition, the calculations are carried out with IGLOO2D. The default options described in [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF] were used for the body-fitted approach, in particular for the MESSINGER2D solver and the boundary-layer solver SIM2D. Regarding the meshes, unstructured grids generated by GMSH were systematically used. The wall mesh size is in the range of 1e-3 to 5e-3 chords (with a refinement in the range of 5e-4 chords for blunt trailing edges). These mesh sizes are fairly representative of default mesh sizes used in IGLOO2D. They generally allow obtaining a good trade-off between solution accuracy (predicted ice shape) and computational time.

For all the calculations, the wall mesh size is kept constant near the leading edge and extended over 0.75 chords as shown in figure Fig. 9. This is required when using the IBM in order to avoid the re-meshing during the multi-step process while maintaining an equivalent wall cell size compared to the BF method (with re-meshing).

A multi-step approach is adopted, using 2 to 10 steps. For the IBM approach, the calculations are performed by changing the calculation strategy for the volume solvers EULER2D and TRAJE2D (penalization) and for the ice shape transportation (level-set), all other parameters remaining the same.

-0.1 0 0.1 0.2 0.3

Fig. 9 Example mesh around a NACA23012 with an extended refinement zone near the leading edge

In the following sections, ice shapes are compared using the criteria presented in [START_REF] Wright | Validation Methods and Results for a Two-Dimensional Ice Accretion Code[END_REF]. Namely, two ice shapes are considered in good agreement if they have similar icing limits, ice thickness distribution, horn angle (or maximum thickness angle), maximum ice thickness and overall ice area. Contrary to [START_REF] Wright | Validation Methods and Results for a Two-Dimensional Ice Accretion Code[END_REF], only a qualitative comparison is made.

A. Manufactured Ice Shape

In order to clearly show the behavior of the level-set method against the usual Lagrangian node displacement approach, a manufactured ice shape is used with a fixed ice accretion rate (thickness and time). The ice accretion time is set to 400s and the ice thickness is enforced to 0.02m for every surface node with coordinate x < 0.05m.

This manufactured geometry was first presented in [START_REF] Lavoie | Comparison of thermodynamic models for ice accretion on airfoils[END_REF] and is generated from a NACA0012 airfoil with added artificial ice near the leading edge. The three-horn configuration was selected to obtain multiple flow recirculation zones and create a difficult situation for the ice growth solver because of the presence of highly concave and convex features.

On Fig. 10, the ice accreted three-horn geometry is illustrated with the enforced ice accretion thickness. The level-set solution is represented in blue, indicating the zone with φ < 0. The ice shape generated by the Lagrangian node displacement is shown as a solid black line, where the geometry overlaps can be seen near concave regions of the geometry. The contour φ = 0 is extracted by our edge marching method and represented by the red line with markers.

As observed in Fig. 10, the level-set method automatically handles the geometry overlaps and the extracted contour provides an explicit surface mesh discretization that can be used in the multi-step ice accretion process. For the results presented in this paper, only the ice shape is used as an Immersed Boundary and the clean geometry is still treated using a Body-Fitted approach. As an illustration, the aerodynamic and droplet fields for the rime ice case are shown in Fig. 11 where the IB (the ice shape) is represented by the red line and the solid body is white. In this section, two-step and 10-step ice accretion simulations are performed for the rime ice case 241. The results are compared between the Body-Fitted and IB methods. The wall mesh size is about 2e-3 chords with a refinement to 5e-4 chords at the trailing edge.

In Fig. 12, the pressure coefficients (Cp) and collection efficiency (β) are compared for the BF and IB methods on the 1 st ice layer of a two-step simulation. Ideally, the IB method should reproduce the results obtained with a BF approach. Fig. 12a illustrates a slight mismatch in Cp near the point of maximum suction. Nonetheless, the collection efficiency is very close between the two methods (Fig. 12b). As rime ice accretion is mostly governed by the collection efficiency, it generates very similar ice shapes for the BF and IB methods despite the difference in pressure coefficients (Fig. 13a). The ice shapes are also in good agreement with the experimental results. The experimental ice shape is the so-called MCCS (Maximum-Combined-Cross-Section) [START_REF] Broeren | Ice-Accretion Test Results for Three Large-Scale Swept-Wing Models in the NASA Icing Research Tunnel[END_REF] derived by the experimentalists from the ice scans (it is more or less the envelope of the ice shape). Since the simulation starts from a BF mesh, the 1 st step is not affected by the IBM and thus, the 1 st ice layer should be the same for both methods. However, a difference might be introduced by the geometry evolution solver which can use either a level-set with contour extraction or the Lagrangian node displacement approach. On Fig. 13a, a two-step ice accretion prediction is made where the 1 st ice layer is illustrated with a dashed line. As the generation of the 1 st ice layer is not influenced by the IBM, the Lagrangian and Eulerian (level-set) geometry updates can be compared, showing negligible difference. Thus, discrepancies observed in Fig. 12 for the surface data can be attributed to the IB method

and not to the level-set approach.

x

Fig. 13 Rime ice case 241 -multi-step ice shape predictions

A close-up of the leading edge Cp distribution is given in Fig. 14a to highlight the discrepancy observed near

x/c = 0.06. The difference between the two solutions can be explained by the interaction between the body-fitted wall and the IBM. Figure 14b provides a view of the airfoil's upper surface near the location of interest. For the rime ice case 241, the ice (illustrated in red) is getting thinner as we approach the impingement limits and eventually it merges with the body-fitted surface. The issue comes from the volume penalization method which, in our case, use a 1 st order implementation of the boundary condition. The cells are penalized if φ < 0 at the cell center and the usual physical equations are solved otherwise. The boundary conditions are applied at the cell center, regardless of the location of the IB within the cell. Thus, when the ice shape is thin enough (Fig. 14b), the cells are no longer penalized. In turn, a slip velocity is applied relative to the body-fitted wall instead of a velocity tangential to the ice shape. This premature switch to a BF wall boundary condition explains the difference in Cp near x/c = 0.06 in Fig. 14a. The issue could be solved by implementing a 2 nd order discretization of the penalization terms or by performing a local mesh refinement to obtain a better representation of the ice shape near the impingement limits. However, the current implementation is still able to provide a good prediction of the ice shape in comparison with the BF results (Fig. 13).

For a 10-layer ice shape prediction (Fig. 13b), the solution is still in good agreement with the experimental data for both methods. Increasing the number of steps reduces the thickness of each ice layer and this might affect the behavior of the penalization method for the same reason described earlier (1 st order discretization of the penalization terms). For instance, the penalization method might effectively see the same geometry for 2 consecutive ice layers even though the ice shape has actually moved. This typically occurs if the ice layer is too thin relative to the mesh cells. For the 10-step simulation presented here, the mesh cell size is about the same as the thickness of a single layer, providing good results. For the glaze ice case, the mesh characteristics are the same as for the rime ice case 241. The wall mesh size is about 2e-3 chords with a refinement to 5e-4 chords at the trailing edge.

For this case, there is again a slight mismatch on the Cp distribution (Fig. 15a), but a very similar collection efficiency for both methods (Fig. 15b). As glaze ice accretion is sensitive to the heat transfer coefficient which is in turn driven by the aerodynamics, the mismatch in Cp might explain the slight difference observed on the ice shape Fig. 16a.

The effect of using the Lagrangian node displacement vs. the level-set approach can again be estimated by analyzing the 1 st ice layer on Fig. 16a, where a negligible difference is observed. It suggests that the difference in Cp and β is due to the IBM. By observing Fig. 15a, the Cp distribution corresponds well between the BF and IB methods for x/c < 0.02 and

x/c > 0.06. The zone where the discrepancy occurs is located near the ice accretion limits where the ice shape stops sharply. For the IBM, this results in a detached flow with a recirculation zone (Fig. 16b) while it is not the case for the BF method, explaining the difference. Here, the comparison is made between the two methods with equivalent mesh size. However, this result suggests that the penalization method require a finer mesh near curved features to be equivalent to the BF approach. Although the ice shape prediction for the 2 nd layer is similar for both methods, it does not reproduce the experimental measurements (Fig. 16a). Note that icing experiments carry large uncertainties as well as spanwise variations [START_REF] Wright | Validation Methods and Results for a Two-Dimensional Ice Accretion Code[END_REF][START_REF] Ruff | Quantification of ice accretions for icing scaling evaluations[END_REF].

However, the divergence from the experimental ice shape seems too large to be attributed only to these uncertainties.

Some tests were performed by refining the mesh and manually increasing the wall roughness (by a factor 2), without significant improvement. Here, the use of a droplet size distribution might help in obtaining a prediction towards the experimental ice shape. This option was, however, not tested as it is not yet available for use with our penalization method.

By increasing the number of ice layers to 10 (Fig. 17a), the ice shape prediction is still far from the experimental results. Moreover, an ice horn is created but not in the same location. When comparing the BF and IB methods, the ice shape is similar for the most part, but with a larger difference near the ice horn (where the effect of the aerodynamics becomes more dominant on the ice accretion). The difference in ice shape is due to the combined effect of the penalization and level-set methods compared to the BF and Lagrangian approach (standard approach). In Fig. 17b, all four combination of methods are shown for the 10 th ice layer only. The figure illustrates that the use of the level-set method has only a limited impact while the IB methods have a larger effect on the difference in ice shape. This is similar to the observation made for the two-step ice accretion simulation, where the penalization method requires a finer mesh near curved features to be equivalent to the BF solution. Using an equivalent cell size, the IBM however provides a good approximation, with only a slightly longer ice horn and slightly different thickness distribution on the upper side.

D. Additional Cases on a NACA0012

In this section, 10-step ice accretion calculations are performed on cases 001, 003 and 004 from [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF] to further assess the behavior of the IBM. The calculations are performed on a coarser mesh (wall mesh size of 5e-3), but it is still representative of typical ice accretion simulations with IGLOO2D.

When comparing the ice shape prediction obtained from the IB and the BF methods, a good match is observed for the rime ice case 001 (Fig. 18a), but a larger difference is seen for the glaze ice cases 003 and 004 (Figs. 18b-18c). This is in line with the observation made in the previous sections. Glaze ice shape are more sensitive to the airflow solution and a perfect correspondence is not obtained for the wall data between the two methods (e.g., Cp distribution, Fig. 15a). Numerical experiments performed outside the scope of this paper showed that, in general, a mesh refinement can help the IBM in retrieving the BF solution. Nonetheless, using the IB and the level-set methods still provides a good estimation of the ice shapes when compared to the experimental data (Fig. 18) and body-fitted solution, even on coarser meshes. In this section, ice accretion is performed on the three-element McDonnell-Douglas Airfoil (MDA, Fig. 19) using the icing conditions provided in [START_REF] Petrosino | Ice Accretion Model on Multi-Element Airfoil[END_REF]. This test case is selected to show the flexibility of the IB and level-set methods. A two-step ice accretion simulation is performed with both the BF and IB methods. Here the objective is to see if the IB method can reproduce the BF solution on this more challenging configuration. To do so, a mesh refinement was performed to obtain similar pressure coefficients and collection efficiency on the 1 st ice layer (2 nd time step), as shown in Fig. 20. The wall cell sizes of the resulting mesh are summarized in Table 2. A finer mesh is required on this test case due to the flow separation downwind of the flap. Because a Euler flow solver is used (inviscid), this flow detachment is very sensitive to the mesh size. The current mesh set-up allowed the IB and BF to behave in a similar way (e.g., similar onset of the flow detachment). Although a better ice shape prediction can be achieved using a RANS solver and a droplet size distribution, it was

Table 2 Wall mesh characteristics

shown in [START_REF] Petrosino | Ice Accretion Model on Multi-Element Airfoil[END_REF] that a fair estimation of the ice shape can be achieved using a Euler flow solver and a single size of droplets for this specific test case. The ice accretion results obtained with IGLOO2D for both methods are shown for the flap, slat and main element in Fig. 21. The predicted ice shapes are not so far from the experiment for the flap and slat but are quite different from the expected solution for the main element. Perhaps a simulation involving more ice layers would improve the results. For instance, six steps are used in [START_REF] Petrosino | Ice Accretion Model on Multi-Element Airfoil[END_REF]. In this paper we are concerned about reproducing the ice shapes from the BF method with the IBM and in this regard, the ice shapes (Fig. 21) are in fact similar for both methods. The comparison still exhibits the usual discrepancy due to the accuracy of the IBM for the aerodynamics.

x/c y/c -0.12 -0.

IV. Conclusion

This paper investigates the potential of using an Immersed Boundary Method (IBM) combined with a level-set to alleviate the issues related to the mesh and geometry update in the numerical prediction of ice accretion. It also aims at assessing the accuracy of such a methodology in comparison with the usual body-fitted approach.

A penalization method is applied to the aerodynamics and droplet trajectories. The surface data is extracted using a weighted least square approach in order to use the boundary layer and ice accretion modules. The geometry (the ice shape) is updated using either a Lagrangian or an Eulerian (level-set) approach. A contour extraction process is also described for 2D meshes made of triangles in order to retrieve the explicit definition of the ice-air interface.

Using a manufactured test case, the level-set is shown to automatically handle geometry folding during the ice shape update while the Lagrangian approach fails at providing a useable surface discretization unless a correction is added.

For ice shape predictions, a BF mesh is used for the clean geometry and only the ice shape is treated as an IB. Following this approach, rime and glaze ice test cases from the Ice Prediction Workshop are performed using up to 10 ice layers.

The IBM predicted an ice shape equivalent to the body-fitted approach on the rime ice case. For the glaze ice case, the predicted ice shape is close to the body-fitted solution but exhibits a larger difference where ice accretion is most dependent on the aerodynamics (e.g., near ice horns). The difference is mostly attributed to the accuracy of the IBM and not to the use of the level-set. Additional rime and glaze ice cases on a NACA0012 showed that the current approach (IBM + level-set) provides a fair estimation of the ice shape when compared to both the BF method and the experimental results, even on coarser meshes. Moreover, a 2-step ice shape prediction on the McDonnell-Douglas multi-element airfoil showed that with proper mesh refinement, the IBM combined with the level-set method can reproduce the BF solution on a more challenging configuration.

Although some improvements can be made in terms of efficiency and accuracy, this paper shows the potential of the proposed methodology for automatic multi-step ice shape predictions. Also, the extension to 3D ice accretion is, in theory, straightforward except for the contour extraction process which will require some adaptation to deal with a 2D surface mesh.

Fig. 1

 1 Fig. 1 Sequential call to modules in multi-step icing simulations.

Fig. 2

 2 Fig. 2 Sequential call to modules in a multi-step icing simulation using an IBM.

Fig. 4

 4 Fig. 4 Example of a propagated icing velocity field for a clean NACA23012 airfoil, body-fitted surface. Coordinates and velocity non-dimensionalized by the chord (c).

Fig. 5

 5 Fig. 5 Example of level-set advection from a clean NACA23012 airfoil, single step ice accretion on a BF mesh

Fig. 6

 6 Fig. 6 Example of contour intersection with a triangle cell: edge and vertex cases

Fig. 8

 8 Fig. 8 Example of a level-set extraction on an iced NACA23012 airfoil. red: surface extraction; blue: ice

Fig. 10

 10 Fig. 10 Comparison between Lagrangian node displacement and Level-Set approach with contour extraction

Fig. 11

 11 Fig. 11 Rime ice case 241: solution from the volume solvers around the 1 st ice layer of a two-layer simulation.

Fig. 12

 12 Fig. 12 Rime ice case 241 -comparison of wall surface data on the 1 st ice layer of a two-layer simulation for the BF and IB methods.

  Fig. 14 Rime ice case 241 -illustration of the Cp discrepancy near x/c = 0.06, two-step ice accretion

Fig. 15

 15 Fig. 15 Glaze ice case 242 -comparison of wall surface data on the 1 st ice layer for the BF and IB simulations.

Fig. 16

 16 Fig. 16 Glaze ice case 242 -two-step ice accretion

Fig. 17

 17 Fig. 17 Glaze ice case 242 -10-step ice accretion

Fig. 18

 18 Fig. 18 Additional cases on a NACA0012 airfoil, 10-step ice accretion

Fig. 19

 19 Fig. 19 Global view of the McDonnell-Douglas LB606b airfoil and its experimental ice shape

Fig. 20

 20 Fig. 20 Wall data on the 1 st ice layer (2 nd step) for the McDonnell-Douglas multi-element airfoil (MDA LB606b)

Fig. 21

 21 Fig. 21 Two-step ice accretion on the McDonnell-Douglas multi-element airfoil (MDA LB606b)

Table 1 Simulation Parameters

 1

		Rime 241	Glaze 242	Case 001	Case 003	Case 004	Multi-Element
	Geometry	NACA23012 NACA23012 NACA0012 NACA0012 NACA0012	MDA
	Chord [m]	0.4572	0.4572	0.5334	0.5334	0.5334	0.9144
	AoA [deg]	2.0	2.0	4.0	4.0	4.0	8.0
	Mach	0.325	0.315	0.325	0.317	0.317	0.27
	P static [kPa]	92.528	92.941	101.325	101.325	101.325	101.325
	T static [K]	250.15	266.05	250.7	262.3	262.3	268.2
	LWC [g/m 3 ]	0.42	0.81	0.55	1.0	0.6	0.6
	MVD [µm]	30.0	15.0	20.0	20.0	15.0	20.0
	Icing Time [s]	300	300	420	231	384	360
	Roughness (ks) [mm]	0.4572	0.4572	0.5334	0.5334	0.5334	0.9144

Lagrangian Layer 1 BF Lagrangian Layer 2 IB Level-Set Layer 1 IB Level-Set Layer 2 Clean Geometry (b) Main Element

									0.08
	0
		Experimental Ice
		BF Lagrangian layer 1					0.06
	-0.02	BF Lagrangian Layer 2 IB Level-Set layer 1
		IB Level-Set layer 2					0.04
		Clean Geometry
	-0.04
									0.02
									Experimental Ice
	-0.12 -0.1 -0.08 -0.06	1	-0.08	-0.06 (a) Slat y/c -0.15 -0.04 -0.1 -0.05 0 0.05	-0.02 0.9	0	0.95	0.02	x/c 0.1 BF x/c y/c 0.02 0.04 0.06 0.08 0.12 -0.08 -0.06 -0.04 -0.02 0 1 1.05 1.1	0.14	0.16	0.18

Experimental Ice BF Lagrangian Layer 1 BF Lagrangian Layer 2 IB Level-Set Layer 1 IB Level-set Layer 2 Clean Geometry (c) Flap

Acknowledgments

This work is financed by the Natural Sciences and Engineering Research Council of Canada (NSERC), the Canada Research Chair program and the ONERA. The authors would like to thank Francesco Petrosino for kindly providing the MDA geometry and the experimental ice accretion data.