come

Numerical Simulation of Polyester Stretch/Blow

Molding Process

Fabrice Schmidt, Jean François Agassant, Michel Bellet, G. Denis To cite this version: Fabrice Schmidt, Jean François Agassant, Michel Bellet, G. Denis. Numerical Simulation of Polyester Stretch/Blow Molding Process. Numerical Methods in Industrial Forming Processes, 1992, sophiaantipolis, France. hal-02017895 Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org)


Keywords: C f friction coefficient f, g, l closure functions

Theory and Validation of a 2D Finite-Volume Integral Boundary Layer Method for Icing Applications Charlotte Bayeux 1 and Emmanuel Radenac 2 and Philippe Villedieu 3 ONERA / DMPE Université de Toulouse, F-31055 Toulouse, France A two-dimensional integral boundary layer method is developed to enable fast and economical computations of boundary layer flows. The ultimate goal is to provide some experience for the extension of this method in three dimensions. In this study, the unsteady momentum and kinetic energy integral equations are solved numerically, together with a set of closure relations based on assumed velocity profiles for laminar and turbulent flows. The robustness of the method is ensured by a Finite-Volume formulation based on an upwind scheme and a semi-implicit time discretization. Moreover a control method has been developed in order to avoid the Goldstein singularity. The accuracy of the numerical method in the vicinity of the stagnation point is strongly improved by introducing a consistent corrective source term in the right-hand side of the equation system. The chosen closure relations are validated with test cases of self-similar flows. Numerical results are also compared with those of a full Prandtl equations code for NACA0012, GLC305 and MS317 airfoils test cases to demonstrate the capabilities of the method. Finally, preliminary results are shown proving the ability of the method to deal with iced airfoils even for complex glaze ice shapes. Today a usual and general method to obtain viscous flow solutions is to numerically solve the full Navier-Stokes equations. However, this approach is time-consuming and computationally demanding. This is a serious issue for applications which require a large number of calculations.

In aeronautics, this is for example the case of in-flight icing simulation. Icing simulation numerical tools are used both for design purpose (to improve the safety of flights) and in the certification procedure as means of compliance. Since a huge number of computations must be performed, tools must be very low time consuming. To date, 2D codes are mature enough to be routinely used by aeronautic engineers but their 3D extension still requires some improvements and especially a significant reduction of the computational time.

A way to drastically reduce the cost of viscous flow computations is to use a coupling between an inviscid flow solver and an integral boundary layer solver.An integral boundary layer code is based on the integration of the boundary layer equations in the direction normal to the wall. This integration results in the loss of one space dimension and thus in significant saving in computational time. As far as icing simulations are concerned, the choice of such a method (instead of using a full Navier-Stokes solver) is justified by the fact that the aerodynamic velocity field in the boundary layer is not useful to compute the droplet collection efficiency and the ice accretion growth rate. Only the values of the skin friction and the heat exchange coefficients must be provided by the air flow solver to the accretion solver [START_REF] Gent | Aircraft icing[END_REF] . The ability of the boundary layer solver to accurately compute separated flows is not mandatory either. Indeed ice accretion generally occurs in the leading edge area (region of highest droplet collection efficiency) where the boundary layer is attached and, when separation occurs due to the ice deposit itself, the heat exchange coefficient in the separated zone is much lower than in the attached part and hence does not need to be accurately predicted.

For these reasons, the development of a general and robust 3D integral boundary layer Finite-aeronautic industry with efficient tools for icing applications. This opinion also relies on the fact that today, even for 2D computations, many existing icing tools (LEWICE [START_REF] Wright | User Manual for the Improved NASA Lewis Ice Accretion Code LEWICE 1.6[END_REF], ONICE [START_REF] Guffond | Validation du programme bidimensionnel de captation[END_REF], IGLOO [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF]) use simplified integral boundary layer models, such as the well-known method of Thwaites [START_REF] Thwaites | Approximate calculation of the laminar boundray layer[END_REF],

and nevertheless are able to provide satisfactory results for most icing conditions. However, the 3D extension of these simplified boundary layer models is not straightforward since the formulation on which they are based is only two-dimensional. More general models and numerical methods, which could be formulated both in 2D and in 3D, are necessary. This is the aim of the present paper to address this objective. This work is in continuity with many other studies performed during the last decades to extend the application domain of integral boundary layer models [START_REF] Cousteix | Couche limite laminaire[END_REF][START_REF] Cousteix | Turbulence et couche limite[END_REF][START_REF] Drela | Two-dimensional transonic aerodynamic design and analysis using the Euler equations[END_REF][START_REF] Tai | An Integral Prediction Method for Three-Dimensional Flow Separation[END_REF][START_REF] Mughal | Integral Method for Three-Dimensional Boundary-Layers[END_REF][START_REF] Nishida | Fully simultaneous coupling of the full potential equation and the integral boundary layer equations in three dimensions[END_REF][START_REF] Milewski | Three-dimensional viscous flow computations using the integral boundary layer equations simultaneously coupled with a low order panel method[END_REF][START_REF] Rafael | CFD and Boundary Layer Models with Laminar-Turbulent Transition around Airfoils and a Rough Cylinder: Results Validation[END_REF][START_REF] Van Es | Comparison and Application of Unsteady Integral Boundary Layer Methods using various numerical schemes[END_REF][START_REF] Dipaola | A rotating reference frame, integral boundary layer method[END_REF][START_REF] Lokatt | Robust viscous-inviscid interaction scheme for application on unstructured meshes[END_REF][START_REF] Zhang | A Non-parametric Discontinuous Galerkin Formulation of the Integral Boundary Layer Equations with Strong Viscous/Inviscid Coupling[END_REF].

Here we restrict our attention to the 2D case. The focus is put on the model formulation, the derivation of the corresponding closure relationships, the choice of the numerical discretization method and validation tests. This work is an important step towards the development of a 3D

Finite-Volume integral boundary layer solver for icing applications.

The paper is divided into 3 parts. In the first one, the chosen formulation of the integral boundary layer model is presented and the corresponding closure relationships are derived for both laminar and turbulent cases. The resulting equation system is shown to be hyperbolic in both cases.

Moreover a control method has been developed in order to avoid the Goldstein singularity [START_REF] Goldstein | On laminar boundary-layer flow near a position of separation[END_REF].

Section III is devoted to the Finite-Volume discretization method. Finally validation test cases are presented in section IV. The accuracy of the proposed methodology is assessed by comparing the results against reference analytical solutions (self-similar boundary layers) and results produced by a full boundary layer solver for 2D flows around different airfoils [START_REF] Aupoix | Couches Limites Bidimensionnelles Compressibles. Descriptif et mode d'emploi du code CLICET -Version 2015[END_REF]. The ability of the method to produce correct results in a fast and robust manner on iced geometries is also assessed.

II. Modeling equations

The purpose of this section is to establish an integral boundary layer model based on two integral equations and a set of closure relations for laminar and turbulent flows.

A Finite-Volume semi-implicit method on unstructured surface meshes is employed to ensure the robustness of the method. Moreover, the need to explicitly locate the stagnation point is avoided by solving an unsteady formulation of the integral equations. A good accuracy of the method is required in the vicinity of the stagnation point, where icing occurs. On the contrary, a very accurate solution is not required in regions where the boundary layer is detached even though the method must produce a realistic result. Indeed, icing occurs mainly in the leading edge zone, where boundary layer separations are rare. In addition, the laminar-turbulent transition is generally located upstream of the laminar separation point, making the boundary layer less prone to separation. When it is not the case, the heat transfer coefficient is expected to be rather low in a separated zone. Consequently, heat transfer at the ice/air interface is less significant in the energy balance within a separated zone than in an attached region. Thus, it is assumed in the present work that a non-accurate boundary layer model for separated zones will not significantly change the final ice shape.

The calculation method developped in this study is valid for a 2D unsteady, incompressible, laminar or turbulent boundary layer. Only the dynamic boundary layer is solved here. It is indeed assumed that the heat transfer coefficient is inferred from dynamic quantities as often done in ice accretion codes [START_REF] Gent | Aircraft icing[END_REF]. The external flow is assumed to be known from a steady irrotational inviscid calculation. The curvature effects are neglected. A 1-way approach is adopted because the inviscid flow is not supposed to be significantly affected by the boundary layer.

B. 2D Integral boundary layer equations

Under the assumptions given in section II A, the unsteady and incompressible Prandtl equations can be written as follows:

                   ∂u ∂x + ∂v ∂y = 0 ∂u ∂t + u ∂u ∂x + v ∂u ∂y = - 1 ρ ∂P ∂x + 1 ρ ∂τ xy ∂y ∂P ∂y = 0 (1a) (1b) (1c)
The pressure gradient term can be expressed in terms of velocity gradients thanks to the momentum equation in the x direction written in the external flow region. This external flow is supposed to be steady, incompressible and inviscid:

-

1 ρ ∂P ∂x = - 1 ρ ∂P e ∂x = u e ∂u e ∂x
The 2D integral boundary layer (IBL) equations are obtained by integrating the following 2D equations with respect to the coordinate y, normal to the wall, over the boundary layer thickness:

u e × continuity equation (1a) -momentum equation (1b) (2a) u e 2 × continuity equation (1a) -2u × momentum equation (1b) (2b)
These equations involve incompressible integral quantities which are defined as:

δ 1 = ∞ 0 1 - u u e dy θ = ∞ 0 u u e 1 - u u e dy δ 3 = ∞ 0 u u e 1 - u 2 u e 2 dy
δ 1 , θ and δ 3 correspond to the displacement, momentum and kinetic-energy thicknesses respectively.

These quantities represent the deficit of mass, momentum and kinetic energy within the boundary layer. Finally, the system of two integral equations for unsteady, incompressible, two-dimensional flows is obtained in conservation form:

∂U ∂t + ∂F (x, U ) ∂x = S(x, U ) (3)
where

U =     u e δ 1 u 2 e θ     =     U 1 U 2     , F =     u 2 e θ u 3 e δ 3 -u 3 e θ     , S =     -u e δ 1 ∂u e ∂x + 1 2 u 2 e C f (u 2 e δ 1 -u 2 e θ) ∂u e ∂x - 1 2 u 3 e C f + 2u 3 e C D    
C f and C D denote respectively the friction coefficient and the dissipation coefficient: Model (3) is not the only possible integral model for 2D boundary layers. Several other formulations can be found in the literature. Some authors prefer to use the entrainment equation instead of the kinetic energy equation [START_REF] Smith | An Integral Prediction Method for Three-Dimensional Compressible Turbulent Boundary Layers[END_REF][START_REF] Coenen | Viscous-Inviscid Interaction with the Quasi-Simultaneous Method for 2D and 3D Aerodynamic Flow[END_REF] and many of them prefer to use a steady formulation [START_REF] Thwaites | Approximate calculation of the laminar boundray layer[END_REF][START_REF] Mughal | A Calculation Method for the Three-Dimensional Boundary-Layer Equations in Integral Form[END_REF] which is simpler and cheaper from a numerical point of view. But, even if it is more complicated than other models, there are important advantages of using model ( 3), especially if one takes into account that the final objective is to extend the chosen formulation to three dimensional problems.

1 2 u e 2 C f = τ xy w ρ u e 3 C D = 1 ρ ∞ 0 τ
First, since model ( 3) is an unsteady model, it can be solved iteratively using a time discretization method until convergence to the steady solution, starting from an arbitrary initial condition. It is not necessary to use a space marching algorithm, from upstream to downstream, which is difficult to generalize to 3D unstructured meshes and arbitrary configurations with several stagnation points or lines. Second, model ( 3) is based on balance equations for both the displacement thickness and the momentum thickness respectively and these two integral length scales are the most commonly used quantities in semi-empirical closure models, in particular for the skin friction coefficient C f and the viscous energy dissipation coefficient C D which are often expressed as function of θ and

H = δ 1 /θ.
The integration process used to derive (3) results in a loss of information especially on the velocity field within the boundary layer. Thus, system (3) involves three unknowns (u

3 e δ 3 , C f , C D )
in addition to the two primary ones (u e δ 1 , u 2 e θ). Closure relations are often derived by making assumptions on the shape of the velocity profiles. It is worth mentioning that only the closure relations are affected by the laminar or turbulent boundary layer flow regime, whereas the two integral equations are valid for both cases.

C. Laminar and turbulent closure relations

The choice we made for the closure relations is presented in this section. It is common practice to close all terms of equation ( 3) by the use of relations between C f , C D and u e 3 δ 3 and the shape factor defined as H = δ 1 /θ. They are often based on self-similar solutions (Cousteix [6], Drela [START_REF] Drela | Two-dimensional transonic aerodynamic design and analysis using the Euler equations[END_REF]), especially for laminar flows. For turbulent flows, empirical relations are also frequently used.

However, the method presented here is also intended to be used for 3D and thermal boundary layer [START_REF] Bempedelis | A 3D Finite-Volume Integral Boundary Layer method for icing applications[END_REF][START_REF] Bayeux | Méthode intégrale pour la couche limite tridimensionnelle -Applications au givrage[END_REF]. Thus, closure relations employed here are based on assumed velocity profiles because they will be easier to implement for the additional equations. Once the velocity profile is known, the computation of u 3 e δ 3 , C f and C D is straightforward.

Laminar closure relations

Many profile families have been proposed in the literature. Several of them can be found in Rosenhead [START_REF] Rosenhead | Laminar Boundary Layers[END_REF]. The velocity profile proposed in our case has the same form as those proposed by

Mangler [START_REF] Mangler | Das Impulsverfahren zur näherungsweisen Berechnung der laminaren Reibungsschicht[END_REF]:

u u e = u = 1 -(1 + aη)(1 -η) p-1 (4)
where η = y δ and δ is a measure of the boundary layer thickness. More precisely, y = δ is the location from which u = u e . This condition strictly applies at y -→ ∞ but it is assumed that it can be transferred from infinity to y = δ without significant error. Boundary conditions ( 5), ( 6) and [START_REF] Cousteix | Turbulence et couche limite[END_REF] are naturally satisfied by the velocity profile (4):

η = 0 : u = 0 (5) η = 1 : u = 1 (6) η = 1 : ∀n ≤ p -2 ∂ n u ∂η n = 0 (7)
One can notice that the exponent p must be strictly greater than 2 by construction. It is worth mentioning that p needs not necessarily be an integer. In this case, the number of zero gradient conditions at η = 1 is given by the largest integer less than or equal to the exponent p -2.

Three other conditions (8), ( 9) and [START_REF] Mughal | Integral Method for Three-Dimensional Boundary-Layers[END_REF] are required in order to express a, p and δ as functions of the primary unknowns u e δ 1 and u e 2 θ. They read:

∂ u ∂η η=0 = u e δC f 2ν (8)
δ 1 = δ 1 0 (1 -u) dη (9) θ = δ 1 0 u(1 -u) dη (10)
where the friction coefficient in condition ( 8) is given by a closure relation proposed by Cousteix [START_REF] Cousteix | Couche limite laminaire[END_REF],

which was extracted from the self-similar solutions of the Falkner-Skan boundary layer equation. It reads:

C f 2 Re θ = g(H) with Re θ = u e θ ν (11)
where Fig. 1 Closure function g against the shape factor H Conditions ( 9) and ( 10) are essential to ensure the consistency of the two primary unkonwns u e δ 1 and u e 2 θ. However, condition ( 8) is not commonly used in the literature. Since the friction coefficient is a major output of the IBL code, this option has been preferred over Mangler choice to use a condition which ensures an accurate modeling of the effect of the pressure gradient. Thanks to the velocity profile (4) and the definition of the shape factor H = δ 1 θ , the set of conditions (8), ( 9) and (10) becomes:

g(H) = 2
p -1 -a = Hg(H) δ δ 1 (13) δ 1 = δ(p + 1 + a) p(p + 1) (14)
θ = δ p + 1 + a p(p + 1) - 2a 2 + 2a(2p + 1) + 2p(2p + 1) (2p -1)2p(2p + 1) (15)
This forms a system of 3 equations with 3 unknowns a, δ and p. From condition [START_REF] Rafael | CFD and Boundary Layer Models with Laminar-Turbulent Transition around Airfoils and a Rough Cylinder: Results Validation[END_REF] More details are given in the appendix on the derivation of the velocity profile and the function p(H). Finally the velocity profile can be written as a function of a single parameter, H, which is directly related to the solved variables u e δ 1 and u e 2 θ:

u(η) = 1 -(1 + a(H)η)(1 -η) p(H)-1 (17)
where a(H) = p 2 -p(p + 1)Hg(H) -1 and p(H) is given by equation [START_REF] Lokatt | Robust viscous-inviscid interaction scheme for application on unstructured meshes[END_REF]. g(H) is given by equation [START_REF] Milewski | Three-dimensional viscous flow computations using the integral boundary layer equations simultaneously coupled with a low order panel method[END_REF]. η = y δ with δ = δ 1 p(p + 1)

p + p 2 -p(p + 1)Hg(H) .

Velocity profiles produced by [START_REF] Zhang | A Non-parametric Discontinuous Galerkin Formulation of the Integral Boundary Layer Equations with Strong Viscous/Inviscid Coupling[END_REF] are compared with several self-similar Falkner-Skan solutions [START_REF] Cousteix | Couche limite laminaire[END_REF] in figure 3. These self-similar solutions are obtained for external velocities expressed as u e = kx m .

The pressure gradient is thus dependent on m or β = 2m 1 + m : it is positive for negative values of β. β = 0 corresponds to a flow over a flat plate (zero pressure gradient). Then the pressure gradient is negative for positive values of β. One may notice that β = 1 represents a 2D stagnation point flow.

Figure 3 shows a fair agreement with the Falkner-Skan profiles especially for β > 0 corresponding to accelerated flows in an attached boundary layer. Thanks to the velocity profile [START_REF] Zhang | A Non-parametric Discontinuous Galerkin Formulation of the Integral Boundary Layer Equations with Strong Viscous/Inviscid Coupling[END_REF], δ 3 and C D can now be easily determined from their definition (regarding the other source terms, C f is given by equation ( 11)):

δ 3 = δ 1 0 u(1 -u 2 ) dη C D = 1 δ 1 0 ν u e ∂ u ∂η 2 dη
It can be shown that the closure relations can be written in the following form:

δ 3 θ = f l (H) 2C D f l (H) Re θ = l(H)
where the functions f l (H) and l(H) are given in the appendix.

Turbulent closure relations

The turbulent boundary layer is composed of an inner region close to the wall where the viscosity prevails and an outer region where turbulence is predominant. Complete velocity profiles, valid over the entire boundary layer, can be found in the literature. One of the most commonly used profiles for IBL methods (Drela [8], Mughal [START_REF] Mughal | A Calculation Method for the Three-Dimensional Boundary-Layer Equations in Integral Form[END_REF], Nishida [START_REF] Nishida | Fully simultaneous coupling of the full potential equation and the integral boundary layer equations in three dimensions[END_REF], Milewski [START_REF] Milewski | Three-dimensional viscous flow computations using the integral boundary layer equations simultaneously coupled with a low order panel method[END_REF]) was developped by Swafford [START_REF] Swafford | Analytical approximation of two-dimensional separated turbulent boundary layer velocity profiles[END_REF] and Whitfield [START_REF] Whitfield | Analytical description of the complete turbulent boundary layer velocity profile[END_REF]. But this model is very cumbersome. The velocity profile chosen in our case is simpler and often used as well for boundary layers (Tai [9], Mughal [START_REF] Mughal | Integral Method for Three-Dimensional Boundary-Layers[END_REF], Kays and Crawford [START_REF] Kays | Convective heat and mass transfer[END_REF],

Cousteix [START_REF] Cousteix | Turbulence et couche limite[END_REF]). It reads:

u = η 1/n (18)
where η = y δ . This power-law with n = 7 is often used for boundary layers [START_REF] Mughal | Integral Method for Three-Dimensional Boundary-Layers[END_REF][START_REF] Kays | Convective heat and mass transfer[END_REF] and gives a good approximation of the turbulent velocity profile especially over a flat plate. To be more general, in this calculation method, two conditions are imposed in order to express n and δ as functions of the primary unknowns:

δ 1 = δ 1 0 (1 -u) dη = δ n + 1 (19) θ = δ 1 0 u(1 -u) dη = nδ (n + 1)(n + 2) (20)
As for the laminar case, conditions ( 19) and ( 20) ensure the consistency of the primary unknowns.

This implies:

δ = (n + 1)δ 1 (21)
and

n = 2 H -1 (22) u = η (H-1)/2 (23)
where H = δ 1 θ . It is worth mentioning that the velocity profile of equation ( 23) is valid for the outer region but not for the region very close to the wall. However this does not introduce a significant error on the value of integral thicknesses δ 1 , θ, δ 3 since the outer region represents the major part of the boundary layer thickness. But the coefficients C f and C D cannot be directly derived from the velocity profile. Thus they are given by relations based on empirical data and the study of equilibrium boundary layers. These relations are proposed respectively by White [START_REF] White | Viscous fluid flow[END_REF] and Drela [START_REF] Drela | Two-dimensional transonic aerodynamic design and analysis using the Euler equations[END_REF]:

C f = 0.3 e -1.33H (log Re θ ) 1.74+0.31H (24)
C D = H * 2 C f 6 4 H -1 + 0.03 H -1 H 3 (25)
with

H * = 1.505 + 4 Re θ + 0.165 - 1.6 √ Re θ (H 0 -H) 1.6 H if H < H 0 H * = 1.505 + 4 Re θ + (H -H 0 ) 2      0.04 H + 0.007 ln Re θ H -H 0 + 4 ln Re θ 2      if H > H 0 and H 0 = 4 if Re θ < 400
H 0 = 3 + 400 Re θ if Re θ ≥ 400
It is worth mentioning that relations ( 24) and ( 25) have been established in smooth-wall conditions.

The rough-wall models are beyond the scope of the present paper and will be addressed in subsequent papers. It is current practice though in icing codes to employ semi-empirical relations which make it possible to calculate the skin friction coefficient on rough walls from the smooth integral quantities [START_REF] Gent | Aircraft icing[END_REF]. Although not fully satisfying, the method consisting of post-processing the integral quantities produced by the present integral model, including the use of relations [START_REF] Bayeux | Méthode intégrale pour la couche limite tridimensionnelle -Applications au givrage[END_REF] and ( 25), would thus be in line with the state-of-the-art. For improved modelling of the physics, equations ( 24) and ( 25) should be changed to include the effect of roughness (with usual relations like the ones of Kays and Crawford [START_REF] Kays | Convective heat and mass transfer[END_REF] or upgraded relations as proposed in [START_REF] Radenac | An extended rough-wall model for an integral boundary layer model intended for ice accretion calculations[END_REF]), so that the integral model solves the rough-wall-condition boundary layer equations.

The last term to close is δ 3 which can be easily determined from its definition thanks to the velocity profile [START_REF] Bempedelis | A 3D Finite-Volume Integral Boundary Layer method for icing applications[END_REF]:

δ 3 = δ 1 0 u(1 -u 2 ) dη
Thus, the closure relation can be written in the following form:

δ 3 θ = f t (H) = 2(n(H) + 2) (n(H) + 3)
where n(H) is given by ( 22).

Transition criterion

The transition can be imposed at a fixed location or calculated using a criterion. For the current paper, the local criterion proposed by Drela [START_REF] Drela | MISES Implementation of Modified Abu-Ghannam/Shaw Transition Criterion[END_REF] and based on the model of Abu-Ghannam and Shaw was chosen:

Re θ T = 155 + 89 0.25 tanh( 10 H -1 -5.5) + 1 n 1.25 n = -8.43 -2.4 ln τ 100 τ = 2.7 tanh τ 2.7
where Re θ T is the Reynolds number of transition and τ is the turbulence rate expressed in %. This criterion has already been used for icing applications by Fuzaro Rafael et al. [START_REF] Rafael | CFD and Boundary Layer Models with Laminar-Turbulent Transition around Airfoils and a Rough Cylinder: Results Validation[END_REF]. However this criterion is not valid on rough surfaces. Rough-wall criteria will not be employed for the present paper because roughness models are not addressed. However, it must be noted that most authors use the criterion Re k > 600, where Re k is a Reynolds number based on the roughness height k, for icing applications (LEWICE [START_REF] Wright | User Manual for the Improved NASA Lewis Ice Accretion Code LEWICE 1.6[END_REF], IGLOO [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF]). The criterion Re k > 600 will thus be a first step towards realistic rough-ice surface modelling in subsequent papers. Kerho and Bragg explain that the criterion should also be changed for little local Reynolds numbers Re x (up to Re k > 1200).

Moreover, the extension of the transitional region is particularly wide on rough walls [START_REF] Kerho | Airfoil Boundary-Layer Development and Transition with Large Leading-Edge Roughness[END_REF]. A more evolved criterion, using an intermittency function as proposed by [START_REF] Rafael | CFD and Boundary Layer Models with Laminar-Turbulent Transition around Airfoils and a Rough Cylinder: Results Validation[END_REF] for instance, could therefore be considered for even more realistic computations.

D. Nature of the integral equation system

Including the closure relations, the integral boundary layer system (3) now reads:

∂U ∂t + ∂F (x, U ) ∂x = S(x, U ) (26)
with

U =     u e δ 1 u e 2 θ     =     U 1 U 2     F =     U 2 u e U 2 (f (H) -1)     =     F 1 F 2    
where f (H) = f l (H) for a laminar boundary layer and f (H) = f t (H) for a turbulent boundary layer. It is assumed here that the model is consistent, which means that the two conditions ( 9) and ( 10) are exactly met. Hence it is assumed that θ in the term U 2 , obtained by the resolution of the integral boundary layer equations, is exactly equal to θ calculated with equation [START_REF] Mughal | Integral Method for Three-Dimensional Boundary-Layers[END_REF]. Consequently, F 1 is strictly equal to U 2 . The integral system is hyperbolic if the Jacobian matrix ∇ U F (x, U ) of the flux vector has only real and distinct eigenvalues.

∇ U F (x, U ) =     0 1 u e 2 f u e (f -Hf -1)     (27)
where f = df dH .

The Jacobian matrix depends on the closure relation for δ 3 . Thus, the hyperbolicity of the integral system depends on the closure assumptions (here the closure relation f (H)). The eigenvalues λ i of ∇ U F (x, U ) are solutions of the following equation:

λ 2 i -λ i u e (f -Hf -1) -u e 2 f = 0 (28)
Figure 4 illustrates that the dimensionless discriminant of the characteristic polynomial (∆/u e 2 = (f -Hf -1) 2 + 4f ) is always strictly positive for a large range of shape factors. Eigenvalues are also plotted in figure 5 for laminar and turbulent cases. These four graphs show that the integral boundary layer system associated with the closure functions is hyperbolic for both laminar and turbulent flows. Both eigenvalues are expected to be positive for an attached boundary layer and to be of opposite sign for a detached boundary layer. Indeed, the information propagates only in the downstream direction in the first case and in both directions (upstream and downstream) in the second. For the turbulent case, eigenvalues are positive for all values of the shape factor because the velocity profile proposed in this method does not allow negative velocities. For the laminar case, the two eigenvalues are positive for the low values of the shape factor and one of them is negative for H > 4.43. However, one may notice that there is a range of shape factor values between H crit and 4.43 where the friction coefficient is negative (i.e. the boundary layer is detached) and the eigenvalues are both positive. According to Cousteix [START_REF] Cousteix | Turbulence et couche limite[END_REF], the size of this zone depends on the chosen closure relations. For a laminar flow, f = 0 for H = 4.43 and the matrix ∇ U F (x, U ) is singular (equation ( 27)). This is the so-called Goldstein singularity [START_REF] Goldstein | On laminar boundary-layer flow near a position of separation[END_REF]. In practice, a steady solution can not be reached in the presence of a laminar boundary layer separation. A good way to avoid this singularity is to couple the boundary layer solver with an inviscid flow solver thanks to a viscous-inviscid interaction method. However, the inviscid flow solver used in the ONERA icing tool is an Euler solver and a viscous-inviscid interaction method can be time consuming and difficult to implement in 3D. In addition, as already mentioned, it is not mandatory to calculate the separation zone accurately.

An approach proposed by Lokatt and Eller is to change the closure relations to prevent the matrix singularity [START_REF] Lokatt | Robust viscous-inviscid interaction scheme for application on unstructured meshes[END_REF]. Here, a "control method" has been developped independently in order to avoid the Goldstein singularity. This method consists of artificially adding a source term in the equations so as to be able to control the value of the determinant of the Jacobian matrix ∇ U F (x, U ) and thus prevent it from vanishing, with the same consequence as Lokatt and Eller's method. This term must obviously vanish in the attached boundary layer zones so that it does not modify the solution in the areas where an accurate solution is sought. However, the solution in the separation zones will not really correspond to the expected solution while still exhibiting negative skin friction. This is the main limit of this method but it is not too restrictive for icing applications. Indeed, icing occurs mainly in the leading edge zone, where boundary layer separations are rare. This is especially true on rough surfaces like ice because the boundary layer becomes turbulent early and it is thus more resistant to separation. Section IV F, however, shows characteristic rime-ice and glaze-ice shapes, for which the transition has been artificially delayed (because a smooth-wall transition criterion was used) and for which the correction happens to be activated. But even under these conditions, the activation of the correction is very local so that it has a low impact on the overall solution. The Goldstein singularity is thus often avoided in practice, in the steady state. The method proposed here is mainly useful for the transitional phase preceding the steady state.

The system of integral equations, including the control term S c (x, U ), reads:

∂U ∂t + ∂F (x, U ) ∂x = S(x, U ) + S c (x, U ) (29)
where

S c (x, U ) =     u e v θ u e 2 v δ1     =     α(H)u e ∂θ ∂x α(H)u e ∂δ1 ∂x    
and α is a function of H that will be determined later. The control term is treated as a source term but it is intended to modify the flux term so as to avoid the singularity in the equations. Let ∇ U F c (s, U ) be the new Jacobian matrix, based on the modified fluxes. It reads:

∇ U F c (s, U ) =     0 (1 -α) u e 2 (f -α) u e (f -Hf -1)    
The function α(H) can then be considered as a control parameter. Indeed, it allows to have an influence on the sign of the determinant of the Jacobian matrix which becomes:

det = -u e 2 (f -α)(1 -α)
α must vanish for an attached boundary layer (H < H crit ) so that the system remains unchanged when no singularity is expected. In this case, det = -u e 2 f is strictly positive because f l < 0. Thus the determinant must be strictly positive for all values of H. Furthermore, α is chosen strictly less than 1 so that the value of the control terms are not too high. This implies that α must be strictly greater than f l . To account for all these constraints, α is chosen as a function of the form:

α(H) = c 1 + tanh H -H crit 0.25
where c is a constant striclty less than 0.5. Its value is fixed in such a way that the system is never singular (i.e. α > f l ) and is always hyperbolic (i.e. (∆/u e 2 = (f -

Hf -1) 2 + 4(f -α)(1 -α)) > 0)
whatever the value of H. Figure 6 shows that c must meet the condition:

0.01994 < c < 0.02457
c is therefore fixed at 0.020.

Figure 7 shows the eigenvalues of the modified system for a laminar boundary layer. As expected, one can see that they are both positive regardless of the value of the shape factor. It is worth mentioning that their values have not been changed in the attached boundary layer zone. The integral boundary layer equations together with the closure relations are solved with a Finite-Volume method to make extension to 3D general surfaces easier. For the 2D structured approach exposed here, equation ( 26) is integrated over an edge i of length ∆x i :

d dt Ci U dx + F (x i+1/2 , U i+1/2 (t)) -F (x i-1/2 , U i-1/2 (t)) = Ci S(x, U ) dx + Ci S c (x, U ) dx
Let the numerical flux G i+1/2 be an approximation of the flux vector F at the node x i+1/2 of edge i. The above equation thus reads:

∆x i dU i dt + G i+1/2 (t) -G i-1/2 (t) = ∆x i (S i (t) + S ci (t))) (30)
The flux terms are discretized thanks to an upwind scheme. This scheme is valid only when the eigenvalues λ ± /u e are both positive. In practice, this condition is always respected due to the control method developed in section II E. The numerical fluxes are calculated with the upstream or downstream values according to the sign of u e :

If u ei+1/2 > 0 G n i+1/2 =     u e 2 i+1/2 θ i u e 3 i+1/2 (δ 3i -θ i )     If u ei+1/2 < 0 G n i+1/2 =     u e 2 i+1/2 θ i+1 u e 3 i+1/2 (δ 3i+1 -θ i+1 )    
where

δ 3i = f (H i )θ i
It is worth mentioning that the external velocity is an input of the calculation. That is why its value at the nodes is known and can be directly used for the flux calculation. In fact, it can be shown that it is necessary to use the value at the nodes in order to enable the stagnation point zone to be accurately calculated (see section III D). The extension to second order accuracy is possible by the use of a MUSCL approach.

B. Semi-Implicit scheme for the time discretization

An explicit Euler method has been employed for the discretization of the fluxes while the source terms are computed implicitly in an attempt to increase the stability of the method and to allow arbitrary initial conditions. In addition, the term of the control method, S ci , is computed explicitly to be treated similarly to the flux terms. For the same reason, the gradients involved in the term S ci are computed with an upwind scheme.

∆x i U n+1 i -U n i ∆t i + G n i+1/2 -G n i-1/2 = ∆x i S n+1 i + ∆x i S c n i (31) 20
C. Solution of the system

The set of equations ( 31) can be written in the following form:

U n+1 i -∆t i S n+1 i = U n i - ∆t i ∆x i G n i+1/2 -G n i-1/2 + ∆t i S c n i
The resolution is performed by a Newton method. The goal is only to obtain the steady-state solution, not the transient regime. Thus it is not necessary to reach convergence at each time step.

Only one iteration for the Newton method per time step is sufficient. The solution is given by:

U n+1 i = U n i + [I -∆t i ∇ U S n i ] -1 ∆t i Res n i
where

Res n i = - 1 ∆x i G n i+1/2 -G n i-1/2 + S n i + S c n i
and ∇ U S is the Jacobian matrix of the source terms.

D. Treatment of the stagnation point

The stagnation point over an airfoil is defined as the point where the external velocity u e vanishes. The stagnation point is often located at the boundary of the computational domain and it is fixed by a boundary condition. However, in this study, the stagnation point is not explicitly localized and can be located within the computational domain. In this region, the numerical errors can become very important since the calculated quantities (u e 2 θ, u e 3 δ 3 ) are very low and lead to some divisions by very small numbers. This section deals with the adjustment of the numerical scheme in order to prevent numerical errors.

Identification of the issue

Focusing on a small zone very close to the stagnation point, it is assumed that the flow behaves as a bidimensional stagnation point flow (i.e. Falkner-Skan flow with m = 1). In that case, the external velocity varies as u e = kx (where x is the distance to the stagnation point) and the quantities H, δ 1 , θ, δ 3 do not depend on x. The flux term

F (x, U ) =     F 1 F 2     =     u e 2 θ u e 3 (δ 3 -θ)    
is derived as follows:

∂F (x, U ) ∂x =     2θu e ∂u e ∂x 3 (δ 3 -θ) u e 2 ∂u e ∂x     =     2θk 2 x 3 (δ 3 -θ) k 3 x 2     (32)
From a discrete point of view, a bidimensional stagnation point flow should verify the following hypotheses:

u ei = kx i , θ i = θ i-1 , δ 3i = δ 3i-1
Let ∆x i be the length of cell i, the external velocity reads:

u ei-1 = k x i - ∆x i 2 - ∆x i-1 2 , u ei-1/2 = k x i - ∆x i 2 , u ei+1/2 = k x i + ∆x i 2
Using the upwind scheme described in section III A , the first element of the derivative of the flux vector ( 32) is discretized as follows:

u e 2 i+1/2 θ i -u e 2 i-1/2 θ i-1 ∆x i = 2θ i k 2 x i
This discrete result is consistent with the continuous one obtained in [START_REF] Drela | MISES Implementation of Modified Abu-Ghannam/Shaw Transition Criterion[END_REF]. But it is worth mentioning that the external velocity must be taken at the node (i + 1/2) and not at the cell center (i). Indeed, in the latter case, we would have:

u e 2 i θ i -u e 2 i-1 θ i-1 ∆x i = 2θ i k 2 x i 1 2 + ∆x i-1 2∆x i - ∆x i-1 4 - ∆x i 8 - ∆x 2 i-1 8∆x i
In the vicinity of the stagnation point, ∆x is not negligible with respect to x. In fact, ∆x can become greater than x in the cell where u e vanishes. Thus, the discrete derivative of the flux term would not be consistent with the continuous one.

The same reasoning is used for the second element of the derivative of the flux vector [START_REF] Drela | MISES Implementation of Modified Abu-Ghannam/Shaw Transition Criterion[END_REF].

Using the upwind scheme described in section III A , its discrete derivative is obtained as follows:

u e 3 i+1/2 (δ 3i -θ i ) -u e 3 i-1/2 (δ 3i-1 -θ i-1 ) ∆x i = 3 (δ 3i -θ i ) k 3 x 2 i + ∆x 2 i 12 (33)
So for the kinetic energy equation, the discrete form of the derivative of the flux term is not fully consistent with the continuous one obtained in [START_REF] Drela | MISES Implementation of Modified Abu-Ghannam/Shaw Transition Criterion[END_REF] even if u e is taken at the node.

This difference can be explained by the fact that the discrete form of the equality ∂u e 3 ∂x = 3u e 2 ∂u e ∂x is not verified with the upwind scheme used. Far from the stagnation point, ∆x 2 is negligible compared to x 2 and the discretized scheme (33) is a consistent scheme. But near the leading edge, this is no more the case and the numerical error is not negligible.

Finally, the study above showed that the system of discrete equations ( 31) is not fully consistent with the corresponding continuous equations in the vicinity of the stagnation point. As theoretically shown above with the assumption u e = kx, the discretization of the system of equations using the upwind scheme of section III A produces an error. This non-consistency was evaluated as being equal to the discretization of the following term:

Er =     Er 1 Er 2     =     ∂u e 2 ∂x -2u e ∂u e ∂x θ ∂u e 3 ∂x -3u e 2 ∂u e ∂x (δ 3 -θ)    
It is worth mentioning that the discretization of the first component of Er has been shown to be equal to zero for a bidimensional stagnation point flow. However, this term can take non-zero values over an airfoil because the hypothesis u e = kx used in this demonstration is not exactly verified for a stagnation point in a non-similar flow.

Correction of the errors

This non-consistency can be corrected by adding a corrective source term in the right hand side of equation [START_REF] Radenac | An extended rough-wall model for an integral boundary layer model intended for ice accretion calculations[END_REF] to offset the error which appears in the left hand side of equation [START_REF] Radenac | An extended rough-wall model for an integral boundary layer model intended for ice accretion calculations[END_REF]. The corrective source term must therefore be equal to the discretization of Er. The continuous form of the vector Er is obviously equal to zero but its discrete form allows to offset the error which appears when the derivative of the flux vector is discretized. Taking the correction term into account, the system of discrete equations ( 31) reads (S c vanishes in the vicinity of the stagnation point):

∆x i U n+1 i -U n i ∆t i + G n i+1/2 -G n i-1/2 = ∆x i S n+1 i + ∆x i Er n+1 i ( 34
)
where

∆x i Er n+1 i =     u e 2 i+1/2 -u e 2 i-1/2 -2u ei u ei+1/2 -u ei-1/2 U 2 n+1 i u e 2 i u e 3 i+1/2 -u e 3 i-1/2 -3u e 2 i u ei+1/2 -u ei-1/2 U 2 n+1 i f (H n+1 i ) -1 u e 2 i    
It is worth mentioning that the first component of Er has been shown to be useless (i.e. equal to zero) for a bidimensional stagnation point flow. However, this term can take non-zero values over an airfoil because the hypothesis u e = kx used in this demonstration is not exactly verified for a stagnation point in a non-similar flow.

In addition, recasting the term Er 2i as Er 2i = k 2i ∂ue ∂x i (f (H i ) -1)U 2i , it can be shown that k 2i can take huge values very close to the stagnation point, which makes the system unstable. Thus, in order to ensure the stability of the corrected numerical scheme, an empirical CFL condition has been introduced. This condition is based on the characteristic time scale of the correction source term 1/ ∂ue ∂x . Moreover, the time step must be small for large values of k 2 . For these reasons, a good condition on the time step to make the numerical scheme with the correction source term stable is:

∆t i < 2 k 2i ∂u e ∂x i
Section IV B 2 shows how the correction source term allows to strongly increase the accuracy of the numerical method in the vicinity of the stagnation point.

E. Pressure gradient limitation

When the adverse pressure gradient is very important, the velocity gradient takes strongly negative values. In this case, the source terms, which can be simplified as follows:

S =     -u e δ 1 ∂u e ∂x + 1 2 u 2 e C f (u 2 e δ 1 -u 2 e θ) ∂u e ∂x - 1 2 u 3 e C f + 2u 3 e C D         -u e δ 1 ∂u e ∂x (u 2 e δ 1 -u 2 e θ) ∂u e ∂x     = - ∂u e ∂x     1 0 -u e 1         U 1 U 2     (35)
become too large and the calculation blows up. A simplified scalar equation representative of problem (3) becomes indeed:

∂v ∂t + a ∂v ∂x = bv
This equation is discretized as follows:

v n+1 i -v n i ∆t + a v n i -v n i-1 ∆x = bv n i
Thus it can be shown, through a matrix stability analysis, that |1 -a∆t/∆x + b∆t| < 1 is a necessary and sufficient condition so that the numerical scheme can converge to a steady solution.

This requires the following condition on the mesh size:

∆x < a b with b > 0
This condition has been adapted for the integral boundary layer equations system by taking a = λ -, the smallest eigenvalue of the flux jacobian matrix, and b = -∂ue ∂x , the maximum eigenvalue of the simplified source matrix (35). The condition on the mesh size becomes:

∆x < λ - -∂ue ∂x if ∂u e ∂x < 0
In practice, it is complicated to impose a condition on the mesh size. Thus this relation has been turned into a condition on the velocity gradient (i.e. the pressure gradient). It reads:

if ∂u e ∂x < 0 and - ∂u e ∂x > C λ - ∆x then - ∂u e ∂x = C λ - ∆x ( 36
)
where C is a constant which can be set by the user. It has been found empirically that C = 0.1 is a good value in most cases. It is worth mentioning that the solution for attached boundary layers is not modified because this condition applies only for adverse pressure gradients and only when they become very large. But, the solution is obviously incorrect in the cells where the limitation is activated. In practice, a strong adverse pressure gradient can occur mainly at the trailing edge of a clean airfoil or near the leading edge when the angle of attack is important. It can also occur downstream of the horn for airfoil with a horn-ice shape. However, as already mentioned, an accurate solution is not mandatory in these zones where ice accretion occurs rarely. Additionally, as is the case for the Goldstein singularity correction, section IV F shows that, for a significant buildup of ice, the pressure gradient correction can also be activated in the iced area, very locally and with a negligible impact on the solution.

Some results from validation tests are presented in this section. The first ones correspond to the simulation of some self-similar laminar boundary layer flows in order to evaluate the accuracy of the closure relations. It is worth mentioning that the turbulent model has been validated too but it will not be presented here because the turbulent closure relations come from the literature. Then the other results correspond to 2D airfoil computations with laminar-turbulent transition. The IBL method, which has been implemented in the so-called code BLIM2D, is then assessed against an ONERA in-house code (CLICET) based on the full Prandtl equations [START_REF] Aupoix | Couches Limites Bidimensionnelles Compressibles. Descriptif et mode d'emploi du code CLICET -Version 2015[END_REF].

A. Self-similar laminar boundary layer flows

The first step to validate the boundary layer calculation method is to compare the computed results with the Blasius and Falkner-Skan solutions, which are self-similar laminar solutions. For a self-similar flow, the velocity profiles only differ by a scale factor from each other. This type of flows can be obtained if the external velocity varies as: As can be seen from table 1, numerical results are in good agreement with the theoretical solutions. Figure 8 shows that for the 2D stagnation point flow (m = 1), the results are perfectly constant thanks to the method exposed in section III D. But for some other cases like m = 1/3, there is still a little error for small Re x . In practice, exactly constant results are necessary only for the case m = 1 in order to compute correctly the leading edge region over an airfoil profile.

u e =
It is worth mentioning that the values given in The L 2 error norm is given as follows:

Err X = N elmt i=1 (X BLIM 2D i -X CLICET i ) 2 N elmt
The errors obtained for the shape factor and the momentum thickness depending on the grid refinement are reported in table 3. The convergence is not very clear because this convergence study is somewhat limited by the fact that the solution given by CLICET is not the exact solution of the model. In addition this study represents a global convergence study for the whole aerodynamic solver because the inviscid flow must be computed with the same grid as the boundary layer. Moreover, it is worth mentioning that CLICET is a compressible solver whereas BLIM2D is an incompressible solver. Thus errors may converge towards an asymptote which represents the difference between the exact solution of CLICET and the exact solution of the model proposed in this paper.

However, it appears clearly in figures 9 and 10 that the solution is converged for meshes 3 and 4 because the errors no longer decrease. However one can see that the solution obtained with the mesh 2 is very satisfactory too. In this section, for NACA0012 cases, calculations performed on a 1024 grid will be presented.

Finally, non-uniform grids composed of 128 elements and refined in the vicinity of the stagnation 11 that the solution obtained with such a mesh is as converged as the solution obtained with the mesh 2. Such a mesh will also be used in the next sections dedicated to MS317 and GLC305 airfoils and iced airfoils.

Influence of the corrective source term on the stagnation point solution

This section deals with the study of the impact of the correction source term (see section III D) on the solution very close to the stagnation point. To that end, solutions obtained with and without the correction source term are compared in figure 12. The data which correspond to the test case used for this comparison are detailed in table 2.

It appears in figure 12 that the solution obtained without correction is incorrect in the vicinity of the stagnation point (about 6 cells). This is particularly apparent for the shape factor. However it is worth mentioning that the relative errors are broadly the same for the displacement and momentum thicknesses and the friction coefficient. In fact, the relative errors of the four quantities are between Table 4 Test case conditions 5% and 10%. In contrast, when the correction source term is activated, the solution in the vicinity of the stagnation point is very close to the one obtained with CLICET.

Finally, the correction developed allows to strongly improve the accuracy of the method in the vicinity of the stagnation point. This is essential to enable to simulate properly the ice accretion over an airfoil. Indeed, ice forms mainly on the leading edge where the droplet collection efficiency and the thermal exchanges are the most important. That is why an inaccurate calculation of the boundary layer in this zone would lead to significant errors on the final ice shape.

Influence of external velocity and angle of attack

Several calculations were performed with different values of external velocity and angle of attack in order to assess the capabilities of the method. One of them, with relatively high Mach number and angle of attack, is presented in this section. The data of the test case are detailed in table 4.

For a better comparison with CLICET, the transition criterion used for the calculations obtained with CLICET is the same as the one used by BLIM2D (see section II C 3).

Figure 13 shows that the results are in good agreement with the reference solution. The solution in the turbulent region is correct and the laminar region is very accurately computed by the code BLIM2D, including in the vicinity of the stagnation point. As already mentioned, especially for icing applications, it is essential to have good results in this region where thermal exchanges are very important and a little error can change significantly the final ice shape.

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 s/c 1.

C. MS317 cases

In this section, two calculations performed over a MS317 airfoil are presented, the first one with an angle of attack of 0 • and the second one with an angle of attack of 8 • . The data of the test cases are detailed in table 5. The accuracy of the solution is similar to the one observed over the NACA0012 airfoil. In particular, the solutions in the leading edge for laminar boundary layers are still excellent (see figures 14 and 15). Nevertheless it may be noted a slightly larger difference with the CLICET solution in the turbulent zone over the MS317 airfoil. The main reason is that the laminar-turbulent transition does not occur exactly at the same location for CLICET and BLIM2D, which induces a delay in the turbulent zone. In addition, the delay in this case is more visible than over the NACA0012 airfoil because the space step is larger than for the previous airfoil.

D. GLC305 cases

The solvers CLICET and BLIM2D have also been compared over a GLC305 airfoil. The data of the test cases presented in this section are detailed in table 6. The conclusions are broadly the same as for the previous cases (see figures 16 and 17). However, for the case with an angle of attack Fig. 14 Comparison between results of BLIM2D and CLICET over a MS317 airfoil for a chord of 0.914m, a Mach number of 0.2420 and an angle of attack of 0 • of 4.5 • , no comparison can be made on the upper surface (i.e. s/c > 0). Indeed, CLICET predicts a boundary layer separation just after the stagnation point and the solver no longer provides a solution because it can not solve the detached boundary layers. However, according to the solution of BLIM2D, the boundary layer is not detached even if it is strongly decelerated. The reason is that the adverse pressure gradient, which is too large, has been limited in this zone. The solution is therefore inaccurate on a little zone where the boundary layer should be detached. However this shows the robustness of the solver which allows to obtain a solution over the whole airfoil, knowing that an accurate solution is not mandatory for strong adverse pressure gradients.

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 s/c 1. First simulations over iced airfoils are presented in this section in order to show the capabilities of the method to perform calculations over complex geometries. The icing suite used to perform the ice accretion calculations is IGLOO2D [START_REF] Trontin | Description and assessment of the new ONERA 2D icing suite IGLOO2D[END_REF]. Figure 18 illustrates the predictor ice shape obtained for the two sets of external conditions detailed in table 7. The first case corresponds to a rime ice case and the second one to a glaze ice case. A boundary layer calculation has been performed with BLIM2D and CLICET over these two ice shapes. The results for the rime ice case are illustrated in figure 19. The solution obtained by BLIM2D is consistent with the one obtained by CLICET. Generally, calculations over rime ice are not much more difficult than those over a clean profile because the shapes are broadly similar (see figure 18). However, one can see that the leading edge has been thickened by the ice deposit and a strong adverse pressure gradient occurs close to the leading edge (at s/c ≈ ±0.1). The boundary layer is therefore strongly decelerated and the shape factor rises quickly. The difference which appears in this zone between CLICET and BLIM2D is due to the fact that the laminar-turbulent transition is not predicted exactly at the same location by the two codes. It is also worth mentioning that a boundary layer separation is predicted by BLIM2D at s/c ≈ -0.07. The treatment of the Goldstein singularity allows to produce a solution for which the skin friction is locally negative (on one point only). The agreement with CLICET remains very good everywhere else. Besides, the non-uniform mesh used for this calculation is not symmetrical. This explains why the results are not exactly symmetrical.

In contrast, for the glaze ice case, the ice shape is much more complex and it significantly affects the inviscid flow. In particular, for horn-ice shapes, the inviscid code usually predicts a large separation zone just downstream of the horns, which is a numerical issue that has to be addressed with BLIM2D. Indeed, this large separation zone results in backflow edge velocities (see figure 18), which may be a numerical challenge for an integral boundary layer code. That is why CLICET provides a solution only in the vicinity of the stagnation point (see figure 20).

However, BLIM2D is able to provide a steady solution over the entire profile. Again, a laminar boundary layer separation is produced by BLIM2D at s/c ≈ 0.024 where the skin friction becomes locally negative (on one point again). CLICET predicts a laminar separation roughly at the same location. It must be noticed that at s/c ≈ -0.044, BLIM2D predicts a laminar-turbulent transition whereas CLICET predicts a laminar separation. Moreover, one can notice in figure 18 that the external velocity vanishes at two different locations on the lower and upper surfaces, in addition to the stagnation point. The points s/c ≈ 0.1 and s/c ≈ -0.064 correspond to the onset of the large separation zone. The adverse pressure gradient is very strong. This is consistent with the fact that the integral quantities are very large and the friction coefficient becomes negative (see figure 20). The points s/c ≈ 0. well defined in a large separation zone. The skin friction given by BLIM2D here is globally lower than in the attached region. It means that the diffusive transfers to the walls are lower than in the upstream region, making their importance lower for the ice accretion computations. However, there is no evidence that the order of magnitude is correct. That is the reason why Navier-Stokes computations were run on the very same configuration, with the ONERA solver elsA.

For the present elsA computation, the grid of figure 21 was generated. It is composed of 117504 quadrangular cells. The grid was extended 50 chords around the airfoil to ensure proper farfield boundary conditions. The cell refinement is correct for a low-Reynolds approach, since the cell size is everywhere lower than y + = 1 in wall coordinates. In order to make a fine comparison between elsA and BLIM2D on the horn ice shape, the same wall conditions were considered for both codes.

Smooth and adiabatic wall conditions were used and transition was enforced at the same position as for BLIM2D, that is close to the top of the horns (s/c -0.04 and s/c 0.08).

Regarding turbulence, the turbulence rate was considered as 1% and a low value of turbulent viscosity was employed in the far field (a ratio 0.1 between the turbulent viscosity and the molecular viscosity was used). A k-ω Menter SST model was used. A steady local time-stepping approach was used with a backward Euler implicit time scheme (LUssor method, CFL=10). The spatial scheme was a Roe scheme with a minmod limiter and with a Harten entropic correction coefficient set to a very low value, 10 -6 , for the transport of turbulent variables. 200000 iterations were performed to achieve steady-state (a decrease in the residuals of around 9 orders of magnitude was obtained).

Figure 21 shows the flow separation obtained with elsA. The inviscid approach is not expected to catch the flow separation correctly (figure 18). But it is also well known that the flow separation is not easy to capture with RANS approaches either, especially its extent. Here, the streamlines of figures 18 and 21 show that elsA predicts the flow as detached slightly further than IGLOO2D. The Fig. 20 Comparison between results of BLIM2D and CLICET over a GLC305 with glaze ice shape agreement between the two codes on pressure and edge velocity is very good in the attached region (figure 22). At the top of horns (s/c -0.02 and s/c 0.09), there is a strong acceleration and a decrease in pressure. IGLOO2D overpredicts this acceleration, especially for the horn located on suction side (s/c > 0). This may be due to the fact that the boundary layer thickens the airfoil and smoothes the horn shape. In the separated zone, the Navier-Stokes computation captures the pressure plateau, whereas the evolution predicted by the inviscid solver of IGLOO2D is more chaotic.

The pressure levels finally get similar again in the downstream part of the airfoil.

Figure 23 shows that the agreement between the two codes is very good on skin friction in the attached region. The main discrepancy between the two codes is obtained at the top of the upper-side horn (s/c 0.09). This is due to the overestimation of the velocity by IGLOO2D.

BLIM2D logically also overestimated the skin friction at this point. Further downstream, the order of magnitude of skin friction produced by BLIM2D is rather correct as compared to elsA (as well as the location of the separation region, where the skin friction is negative). Skin friction is overall underpredicted by BLIM2D. The under-prediction of skin friction by BLIM2D is a sign that diffusive transfers to the wall are underestimated. The role of the convective transfer on ice accretion should therefore be reduced as well in the large separation area. This is safer than an over-prediction of the diffusive transfers to the wall, given the low level of accuracy of the solution in this area. Figure 24 shows where the Goldstein singularity correction and the pressure gradient limitation are activated on the the two cases studied. In both cases, the Goldstein correction is activated at convergence very occasionally. During the iterative process, the activation area is a little larger (cells where the correction was activated at least once during the iterations). Similarly, the limitation of the pressure gradient can be active until convergence in the iced area. For the rime case, the pressure gradient limitation is activated near the laminar-turbulent transition. For the glaze case, the limitation is particularly used in the flow separation area, mainly after the laminar-turbulent transition. But the limitation is also active in a small area near the stagnation point.

Figure 25 shows that the activation of these corrections has a limited effect on wall friction.

For the rime case, the difference between BLIM2D and CLICET is mainly due to the fact that the laminar-turbulent transition is a little delayed for BLIM2D (s/c = ±0.07 instead of ±0.064). In the areas where both codes produce laminar regime (s/c ∈ [-0.06, -0.055] and s/c ∈ [0.06, 0.062]), the pressure gradient limitation generates a rather small deviation on C f (although it is a little larger than elsewhere). The same observation holds in areas where both codes produce turbulent regime (s/c = 0.074). For the glaze case, the differences observed between BLIM2D and elsA in the flow separation zones have already been commented earlier. The Goldstein singularity correction by BLIM2D. But on both sides of this point, the level of friction is the same for both codes. At s/c 0.04, the pressure gradient limitation is activated and there is no significant difference between the skin frictions produced by the two codes.

No significant errors are therefore created by the activation of corrections, even though it should be noted that the two studied cases are particularly difficult. Indeed, the walls are considered as smooth, which means that the laminar-turbulent transition is delayed compared to a transition on a rough wall. Let us examine for example the case of an advanced transition at s/c = 0.03 instead of s/c 0.07 for the rime case.

It is interesting to note that the Goldstein correction no longer needs to be activated (figure 26(a)), which was expected. But the limitation of the pressure gradient is also almost no longer necessary (except in the immediate vicinity of the trailing edge and locally during the iterative process). This is due to the fact that the limitation is now activated in turbulent regime and the eigenvalue used in equation ( 36) is larger in turbulent regime than in laminar regime. There are two reasons for this. First, the eigenvalues are slightly larger in turbulent regime than in laminar regime (figure 5). Second,the eigenvalues are a decreasing function of the shape factor H and H is lower in turbulent regime than in laminar regime. Figure 26(b) shows that the agreement on the skin friction is overall very good between the two codes.

V. Conclusion

In this paper, an integral boundary layer method developed specifically for icing applications has been presented. In particular, it is known that ice accretion occurs mainly at the leading edge of an airfoil. For this reason, special attention has been given to the resolution of the boundary layer in this zone. Firstly, a closure of the integral model based on a laminar velocity profile has been developed in order to achieve a high precision for the calculation of the integral quantities and the friction coefficient at the leading edge. Then, a specific numerical procedure has been developed to improve significantly the accuracy of the solution in the vicinity of the stagnation point. In contrast, ice deposit occurs rarely in strong adverse pressure gradient and separation zones. Thus a condition to limit adverse pressure gradients has been developed to make the solver robust and to ensure that a steady solution is always provided, even if it is inaccurate in non-critical areas as far as icing is concerned. Moreover, the discretization method is based on a Finite-Volume formulation, adapted for unstructured meshes. This allows to perform calculations over airfoils with ice shapes which are often complex geometries. Finally, the model formulation and the discretization method have been chosen to ensure the possible extension of the method to the three dimensional case [START_REF] Bempedelis | A 3D Finite-Volume Integral Boundary Layer method for icing applications[END_REF]. Many validation test cases have been performed, both on academic and more realistic configurations. The results show the high accuracy of the method at the leading edge and its robustness in strong adverse pressure gradient zones. Calculations over airfoils with complex ice shapes are in progress but the first results exposed in this article are encouraging, both in terms of accuracy in the attached zones and robustness.

The ultimate goal of this work is to implement the IBL method in a 3D icing suite. To that end, the extension of the method to 3D is on-going [START_REF] Bempedelis | A 3D Finite-Volume Integral Boundary Layer method for icing applications[END_REF][START_REF] Bayeux | Méthode intégrale pour la couche limite tridimensionnelle -Applications au givrage[END_REF][START_REF] Blanchard | 3D ice accretion modeling using an integral boundary layer method[END_REF]. Besides, the solver must be able to calculate the friction coefficient and the heat transfer coefficient over an ice shape. So a follow-up to the present work is to solve the thermal boundary layer and to take into account the effects of surface roughness. The heat transfer coefficient can be inferred from dynamic quantities thanks to a Reynolds-like analogy. This approach is commonly used in 2D ice accretion codes. A more accurate method consists of solving a model for the thermal boundary layer by adding the integral energy equation to the system. The rough quantities can either be inferred from the smooth ones thanks to semi-empirical algebraic laws or be directly computed by modifying the closure relations to account for roughness effects. This work is currently on-going and promising preliminary results based on a

Reynolds-like analogy have already been obtained [START_REF] Bempedelis | A 3D Finite-Volume Integral Boundary Layer method for icing applications[END_REF][START_REF] Bayeux | Méthode intégrale pour la couche limite tridimensionnelle -Applications au givrage[END_REF].

One may notice that there are real and physical solutions of the above equation only for Hg(H) < (continuity of p(H))

This system of equations is solved using a least-squares method to approximate the function p(H).

The solution is given by equation ( 16).

  Re () Reynolds number Re () = ue() νe S source term vector U vector of primary variables u, v velocity components (m/s)

xy ∂u ∂y dy where τ xy = µ ∂u ∂y in laminar flow

  regime and τ xy = µ ∂u ∂y -ρ < u v > in turbulent flow regime. The first equation is an unsteady version of the well-known von Kármán equation. The second one is the kinetic energy integral equation. The conservative formulation is mandatory for the Finite Volume resolution.

Fig. 2

 2 Fig. 2 Degree of the analytic velocity profiles p against the shape factor H

16 Fig. 3

 163 Fig. 3 Polynomial velocity profiles compared with those obtained from the Falkner-Skan solutions for different pressure gradient parameters

Fig. 4 Fig. 5

 45 Fig. 4 Dimensionless discriminant of the characteristic polynomial against the shape factor H

Fig. 6

 6 Fig. 6 Regularity (det > 0) and hyperbolicity (∆ > 0) of the system with respect to c

kx m where m = x u e du e dx = β 2

 2 -β This edge velocity results from a wedge flow with an angle βπ. Under these conditions, the quantities H, Re δ1 / √ Re x , Re θ / √ Re x and C f √ Re x /2 are constant for each value of m. Four calculations were performed in order to assess the accuracy of laminar closure relations. They correspond to a stagnation point flow (m = 1), an accelerated boundary layer (m = 1/3), a flow over a flate plate (zero pressure gradient, m = 0) and a decelerated boundary layer (m = -0.07).

Fig. 9

 9 Fig. 9 Shape factor against the dimensionless curvilinear abscissa for different grid refinements

Fig. 11 Fig. 12

 1112 Fig. 11 Use of a non-uniform grid (zoom on the stagnation point)

Fig. 13

 13 Fig. 13 Comparison between results of BLIM2D and CLICET over a NACA0012 airfoil for a chord of 0.5m, a Mach number of 0.3 and an angle of attack of 4 •

Fig. 15

 15 Fig. 15 Comparison between results of BLIM2D and CLICET over a MS317 airfoil for a chord of 0.914m, a Mach number of 0.2420 and an angle of attack of 8 •

Fig. 16

 16 Fig. 16 Comparison between results of BLIM2D and CLICET over a GLC305 airfoil for a chord of 0.9144m, a Mach number of 0.2730 and an angle of attack of 4.5 •

Fig. 17

 17 Fig. 17 Comparison between results of BLIM2D and CLICET over a GLC305 airfoil for a chord of 0.9144m, a Mach number of 0.3940 and an angle of attack of 1.5 •

( a )

 a Fig. 18 Studied ice shapes

Fig. 19

 19 Fig. 19 Comparison between results of BLIM2D and CLICET over a NACA0012 with rime ice shape

  Fig. 21 Navier-Stokes elsA computations of glaze-ice GLC305 test-case

Fig. 22 Fig. 23

 2223 Fig. 22 Profiles produced by elsA and IGLOO2D around the glaze-ice GLC305 airfoil

Fig. 24

 24 Fig. 24 Locations where Goldstein singularity correction (Goldstein corr.) and pressure gradient limitation (PG Lim.) were activated during the BLIM2D computations of the cases of table 7

Fig. 25

 25 Fig. 25 Skin friction predicted by BLIM2D and the reference code and locations where Goldstein singularity correction and pressure gradient limitation were activated during the BLIM2D computations of the test-cases of table 7. : Goldstein singularity correction activated during the iterative process, : Goldstein singularity correction activated at convergence, : pressure gradient limitation activated during the iterative process, : pressure gradient limitation activated at convergence

Fig. 26

 26 Fig. 26 Skin friction predicted by BLIM2D and the reference code and locations where Goldstein singularity correction (Goldstein corr.) and pressure gradient limitation (PG Lim.) were activated during the BLIM2D computations of the rime-ice case of table 7. Laminar-turbulent transition imposed at s/c = 0.03. : Goldstein singularity correction activated during the iterative process, : Goldstein singularity correction activated at convergence, : pressure gradient limitation activated during the iterative process, : pressure gradient limitation activated at convergence.

  Function g is plotted against the shape factor H in figure1. It is apparent from this figure that g(H) (i.e. C f ) is negative for shape factors greater than H crit = 4.02923. This corresponds to areas where the boundary layer is detached from the body. In contrast, areas where the boundary layer is attached are characterized by a shape factor lower than H crit .

	.99259	1 H	-	1 8.05846	1.7	-	1 8.05846	1.7	if H ≤ H crit	(12a)
	g(H) = 0.20644 -90.30936	1 4.02923	1.3	-	H 1.3 1	3.35661
	+ (H -1) -0.06815 + 46.34236	1 (4.02923) 2 -	1 H 2	2.338238	if H > H crit	(12b)
					0.4
				g(H) = Re	0.2
					0
					0		2		4		H	6	8	10

θ Cf/2

table 1

 1

	are those obtained from a sufficiently large

Table 1

 1 Comparison of numerical results with analytical solutions for different Falkner-Skan

	flows
		1.05		1.05
		1.04		1.04
	H/H theo	1.02 1.03	H/H	1.03 1.02
		1.01		1.01
		1		1
		1e+03 1e+04 1e+05 1e+06 1e+07 1e+08		1e+03 1e+04 1e+05 1e+06 1e+07 1e+08
		Re x		Re x
	(a) Accelerated boundary layer (m = 1/3)
	B. NACA0012 case
	In this section, some results obtained over a NACA0012 airfoil are presented. The integral
	boundary layer code BLIM2D is compared with the code CLICET which is considered as the refer-
	ence.
	The fluxes are evaluated by a first order scheme and an arbitrary constant inital condition is
	imposed at the beginning of the calculation.

theo (b) 2D stagnation point flow (m = 1) Fig. 8 Ratio between H and H theo against the local Reynolds number Rex

Table 2

 2 Test case conditions 1. Mesh convergence A preliminary mesh convergence study was carried out. Four uniform meshes were constructed with a constant grid refinement ratio. They are composed of 128, 256, 512 and 1024 elements. The data of the test case used for this study are detailed in table 2. The mesh convergence study is limited to the leading edge of the airfoil, where the boundary layer is laminar.

Table 3

 3 Errors for the shape factor and the momentum thickness depending on the grid

	refinement
			CLICET 1024 elmts
	2.7		BLIM2D 128 elmts BLIM2D 256 elmts
			BLIM2D 512 elmts		2.34
	2.6		BLIM2D 1024 elmts		2.32
	2.5						2.3
	H					H	2.28
	2.4
							2.26
	2.3						2.24
							2.22
	2.2	-0.2	-0.1	0	0.1	0.2	-0.04 -0.02	0	0.02	0.04
				s/c				s/c
			(a) Global view			(b) Zoom on the stagnation point

Table 6

 6 Test cases conditions E. Calculation over airfoils with ice shapes

  29 and s/c ≈ -0.35 correspond to the end of the large separation zone. The flow behaves like a stagnation point: H ≈ 2.21 and the friction coefficient is large.The BLIM2D solution in the separation region is meaningless because the boundary layer is not

	Case	Rime ice Glaze ice
	Profile	NACA0012 GLC305
	Chord (m)	0.533	0.914
	AOA ( • )	0	0
	M∞	0.2807	0.3660
	T∞ (K)	252.44	260.78
	P∞ (P a)	93000	77000
	LWC (g/m 3 )	0.51	0.58
	MVD (µm)	40	44
	t (s)	600	876

Table 7

 7 Test cases conditions

  Finally, p can be obtained as follows. Knowing that H = δ 1 /θ, equation[START_REF] Dipaola | A rotating reference frame, integral boundary layer method[END_REF] becomes:Thanks to (37) and (38), the above equation depends only on p and H. This equation can then be solved numerically to obtain p(H) (see figure2, curve "target values"). Let's consider a piecewise function to fit this curve of target values (because g(H) is a piecewise function):p(H) = x 1 + x 2 (H -1.9538) x3 if H ≤ H crit p(H) = 2 + x 4 (H -x 5 ) x6 if H > H critwhere H crit = 4.02923. The expression above is justified by the fact that p must tend to +∞ when H tends to 1.9538 by a higher value. In addition, p must always be greater than 2, particularly when H tends to +∞. Moreover, p(H) and ∂p ∂H must be continuous in H crit hence:x 5 = H crit -x 4 x 6 (H crit -1.9538) x3+1x 2 x 3

							1
	(i.e. H > 1.9538) and p >	Hg(H) 1 -Hg(H)	. Under these conditions, the solutions for	δ δ 1	are:
				δ δ 1	=	p ± p 2 -p(p + 1)Hg(H) Hg(H)	=	p(p + 1) p ∓ p 2 -p(p + 1)Hg(H)
	The ratio	δ δ 1	must be positive for all values of H. Therefore, the solution is:
							δ δ 1	=	p(p + 1) p + p 2 -p(p + 1)Hg(H)	(38)
					1 H	=	δ δ 1	p + 1 + a p(p + 1)	-	2a 2 + 2a(2p + 1) + 2p(2p + 1) (2p -1)2p(2p + 1)
							1
							x 6 + 1	(continuity of p'(H))
	x 1 = 2 +	x 4 (H crit -x 5 ) x6 -	x 2 (H crit -1.9538) x3

Appendix Derivation of the laminar velocity profile

The laminar velocity profile has the following form:

where a, δ, and p must be expressed as functions of the primary unknowns thanks to the set of conditions ( 13), ( 14) and [START_REF] Dipaola | A rotating reference frame, integral boundary layer method[END_REF] which is recalled below:

From condition [START_REF] Rafael | CFD and Boundary Layer Models with Laminar-Turbulent Transition around Airfoils and a Rough Cylinder: Results Validation[END_REF], one can easily deduce a:

Then, δ δ 1 can be expressed as a function of H and p. Thanks to (37), condition [START_REF] Van Es | Comparison and Application of Unsteady Integral Boundary Layer Methods using various numerical schemes[END_REF] becomes:

Derivation of the closure functions

The closure functions for δ 3 and C D are obtained by integrating the veloctiy profile [START_REF] Zhang | A Non-parametric Discontinuous Galerkin Formulation of the Integral Boundary Layer Equations with Strong Viscous/Inviscid Coupling[END_REF]. The definition of these quantities is recalled below:

The closure functions f l (H) and l(H) are defined as:

where Re θ = u e θ ν . Thus, they read:

where

-1 p(4p 2 -1)

where a(H) = p 2 -p(p + 1)Hg(H) -1 and p = p(H) is given by ( 16).