Mesh Impact Analysis of Eulerian Method for Droplet Impingement Characteristics Under Aircraft Icing Conditions

，，，*，

1.Laboratory of Fundamental Science on Ergonomics and Environmental Control，School of Aeronautic Science and Engineering，Beihang University，Beijing 100191，P.R.China；

2.Shenyang Key Laboratory of Aircraft Icing and Ice Protection，AVIC Aerodynamics Research Institute，Shenyang 110034，P.R.China

Abstract: The research of super-cooled water droplet impingement characteristics is the basis of the ice shape prediction and the aircraft anti-icing system performance analysis.The Eulerian method is frequently used to compute the droplet motion and impingement characteristics，so it is significant to analyze the accuracy of the Eulerian method when calculating the droplet impingement on complex surfaces.Taking a NACA 0012 airfoil，an icing wind tunnel，an S-shape duct and a multi-element airfoil as the research objects，the Eulerian method is used with different meshes to obtain the water droplet motion and collection efficiency.The results show that when the water droplet is not deflected or blocked by upstream components，the results obtained by the Eulerian method are slightly affected by the mesh.However，the results calculated by the Eulerian method are mesh-dependent with upstream trajectory deflections，and the liquid water content and collection efficiency calculated using different meshes are inconsistent.The impact of the mesh on the calculation of the Eulerian method needs to be considered when droplets are affected by upstream effects.The findings of this research are beneficial for the accuracy of aircraft icing simulation.

Key words：aircraft icing；droplet impingement characteristics；Eulerian method；mesh impact；numerical simulation

0 Introduction

Super-cooled water droplets in clouds hit the windward surfaces of the aircraft components，leading to the aircraft icing during flight.Aircraft icing can seriously affect the aerodynamic performance，stability and maneuverability of aircraft，and even cause flight accidents［1］.To ensure the flight safety，it is essential to investigate the aircraft icing process and design anti-icing systems［2］.Since the impingement of super-cooled water droplets on the windward surfaces is the initial cause of aircraft icing，the numerical simulation of the droplet impingement characteristics is the basis for the analysis of aircraft icing process［3］and the design of ice protection system［4］.The collection efficiency，which is the ratio of surface to free-stream water mass fluxes，is always used to represent the droplet impingement characteristics［5］.

The super-cooled water droplet motion and impingement process is an air-droplet two-phase flow problem，and the Eulerian method is one of the main methods to solve it［6］.The water droplets in the cloud are considered as a continuous phase in the Eulerian method，and in order to determine the spatial distributions of the droplet volume fraction and velocity，the continuity and momentum equations of the droplet phase are solved［7］.Then，the collection efficiency of each surface can be determined.Yet，due to the highly nonlinear feature of the governing equations for droplets，the numerical instability in the Eulerian method could be caused by the local singularity of the droplet concentration［8］.Thus，in order to prevent droplet oscillations and guarantee convergence，numerical diffusion［8］and artificial viscosity［9］should be introduced.The Eulerian method has so far been frequently used for the water droplet impingement characteristics of complex surfaces and rotating components［9-11］.

However，it has been found that the motion of supercooled water droplets can be deflected by the frontal body along with the airflow，and some droplet trajectories can even overlap or cross［12］.The characteristics of water droplets such as the spatial distribution and the distance between droplets differ from those in the free stream，contributing to the difference of droplet impingement characteristics on the downstream surface，which is known as the upstream effect.Additionally，it has been discovered that the particle trajectory crossing effect cannot be captured by the Eulerian method in the simulation of the granular flows［13-14］，resulting in the phenomenon of “delta shocks”［15］.As the Eulerian method cannot simulate the droplet trajectory crossings，the accuracy of the Eulerian method should be further investigated.

Meanwhile，the mesh plays a significant part in the numerical simulation of computational fluid dynamics（CFD），and the impact of mesh generation on the Eulerian method for the droplet impingement characteristics needs to be considered.With the absence of upstream effect，the Eulerian method is mesh-independent in the calculation of droplet impingement characteristics，and the simulation results obtained with different meshes are well consistent［6，16］.However，there is a lack of research on the mesh impact of the Eulerian method for droplet impingement characteristics with upstream effects.Thus，further research is needed to investigate the mesh impact on the calculation results obtained by the Eulerian method in this case.

To sum up，the Eulerian method has been extensively used in the numerical simulation of supercooled water droplet impingement characteristics for aircraft icing.However，some simulations of aircraft icing and granular flow with the Eulerian method have been found to be flawed.And there are few researches on the mesh impact of the Eulerian method for the droplet impingement characteristics with upstream effects.Therefore，in order to analyze the Eulerian method for the water droplet impingement characteristics of aircraft icing，this paper uses the typical Eulerian method software FENSAP-ICE to simulate the droplet motion and impingement characteristics.The goal of this research is to analyze the differences between the results of different meshes，to investigate the mesh impact on the Eulerian method for the droplet impingement characteristics，and to analyze the accuracy of the Eulerian method.Firstly，the model of the Eulerian method is established with assumptions to calculate the local collection efficiency.The Eulerian method is verified with a NACA 0012 airfoil by meshing with different densities，and the analysis of mesh impact without upstream effects is performed.Then，an icing wind tunnel is used to calculate the motion of water droplet with different mesh densities and analyze the mesh impact on the prediction results after water droplet deflection.Finally，the accuracy of the Eulerian method for the droplet impingement characteristics is analyzed by using an S-shape duct and a multielement airfoil.

1 Numerical Simulation Model and Method

1.1 Mathematical model

In the research of aircraft icing，air is treated as the carrier phase，while the super-cooled water droplet is treated as the dispersed phase.The diameter of water droplets in clouds is about 20 μm，the liquid water content（LWC）is about 1 g/m3，and the droplet volume fraction is on the order of 10-6，which is very low［17］.Therefore，the air-droplet twophase flow can be thought of as being one-way coupled［18］.This means that the airflow field affects the droplet motion through aerodynamic drag，but the droplet motion does not have an influence on the airflow field.So，the airflow field is solved independently and can be directly used in the droplet calculation.The following assumptions are delivered for the motion and impingement characteristics of water droplets［18］：（1）Droplets are spherical with no deformation or breaking；（2）droplets do not have mutual influence without collision，coalescence or splashing；（3）droplets do not exchange heat and mass with the airflow；（4）the effects of turbulence on the droplets are disregarded；（5）aerodynamic drag，gravity and buoyancy are treated as the only forces acting on the droplets.

Based on the above assumptions，the continuity and momentum droplets equations in the Eulerian method can be written as［18-19］

whereαis the droplet volume fraction，tthe time，uthe droplet velocity vector，uathe air velocity vector，μathe air viscosity，ρthe droplet density，dpthe droplet diameter，ρathe air density，CDthe drag coefficient of the droplet，Frthe Froude number andRethe relative Reynolds number given as

The spatial field of the droplet volume fraction and velocity can be obtained by solving the above equations.Then，the local collection efficiencyβcan be calculated by

whereunis the droplet velocity in the impingement control volume，u∞the freestream velocity，αnthe droplet volume fraction on the surface of the control volume，α∞the droplet volume fraction in the freestream，andnthe unit normal vector on the surface of the control volume.

1.2 Implementation method

As the air-droplet two-phase flow can be considered to be one-way coupled，the field of airflow and droplet flow are separately obtained using the CFD software ANSYS FLUENT 19.2［20］and the typical Eulerian method software FENSAP-ICE 19.2［19］.The airflow fields of the NACA 0012 airfoil，the icing wind tunnel and the S-shape duct are solved by Reynolds-averaged Navier-Stokes（RANS）equations with thek-ωSST turbulence model，whereas the airflow field of the multi-element airfoil is solved by RANS equations with the Spalart-Allmaras turbulence model.In FLUENT，the finite volume method is utilized to calculate the airflow filed，the semi-implicit method for pressurelinked equations（SIMPLE）method is used as the solution method，and the second order upwind scheme is chosen as the difference scheme.Then，the airflow field obtained by FLUENT is read into FENSAP-ICE，and the DROP3D module is used to solve the Eulerian continuity and momentum droplets equations.In DROP3D，the droplet velocity at far fields is set to be consistent with the air velocity，and the droplet diameter and LWC are set as environmental parameters，after which the droplet flow field can be solved.

2 Results and Discussion

2.1 Method validation

To verify the correctness of the present Eulerian method，and to analyze the mesh impact without upstream effects，a NACA 0012 airfoil having a unit chord is used to compute the collection efficiency，which is contrasted with the results in Ref.［5］.Since FENSAP-ICE can only solve 3D geometric models，the computational domain of the airfoil is stretched by one mesh in the normal direction，with the two sides set as symmetry boundary conditions.Then，a 2.5D geometric model is established for simulation.As there is always no normal velocity in the calculation of the 2.5D model，the results will be consistent with those obtained with the 2D geometric model.The freestream velocity is 0.4Mawith the angle of attack（AOA）of 5°.The inlet pressure and the ambient temperature are set as 1.01×105Pa and 300 K，respectively.The LWC is 1 g/m3with the droplet diameter of 16 μm.Two Ctype meshes with different densities are used to calculate the validation cases，and the number of nodes at mesh boundaries in Case 2 is 1.5 times that in Case 1.Each mesh has 17 700 and 25 200 elements.A schematic diagram of mesh for the NACA 0012 airfoil is presented in Fig.1.

Fig.1 Schematic diagram of mesh for NACA 0012 airfoil

The contour of LWC around the NACA 0012 airfoil is shown as Fig.2.Both the upper and lower airfoil surfaces have a region where LWC approaches zero，which is the water-shaded region.Since the mass and inertia of water droplets are much larger than those of air particles，the droplets deviate from the air streamline and partly hit the airfoil surface when the airflow changes sharply near the airfoil.And the water droplets without impingement continue to move with the airflow after deflection，forming the airfoil blocking of water droplet motion.Additionally，the airflow change has a larger influence on the deflection of the droplets near the airfoil surface，causing the LWC to be typically higher outside the shadow region.

Fig.2 Contour of LWC for NACA 0012 airfoil

Fig.3 compares the pressure coefficients of the NACA 0012 airfoil obtained by the two meshes，and the two results show good agreement.The pressure coefficient distribution of the NACA 0012 airfoil is independent of the number of mesh elements，which proves the mesh-independence of airflow field.

Fig.3 Comparison of pressure coefficients obtained by two meshes for NACA 0012 airfoil

Fig.4 compares the collection efficiency of the two meshes with that in Ref.［5］.The maximum value of the local collection efficiency appears at the stagnation point of the lower airfoil surface.And the collection efficiency gradually decreases as it proceeds backwards along the airfoil，and drops to zero at the droplet impingement limit.As shown in Fig.4，a good agreement is found between the collection efficiencies obtained by the two meshes and Ref.［5］，validating the present Eulerian method.Furthermore，the collection efficiencies calculated by the two meshes are consistent.It indicates that when the droplet motion is not affected by the upstream effect，the droplet motion and impingement characteristics obtained by the Eulerian method are mesh-independent，verifying the accuracy of the Eulerian method.

Fig.4 Comparison of collection efficiencies with results in Ref.[5] for NACA 0012 airfoil

2.2 Comparison with an icing wind tunnel

According to the analysis in Introduction of the paper，when the water droplet is affected by the upstream components，the direction of movement will be deflected and even the trajectories will overlap or cross.Because the real droplet trajectory crossings cannot be simulated using the Eulerian method，the accuracy of the Eulerian method should be analyzed in these situations.An icing wind tunnel is used to simulate the airflow field and water droplet motion.The geometric model is shown in Fig.5，the air flows in the +xdirection，and the computation domain is 2 m long with the inlet height of the wind tunnel being 1 m.As mentioned above，the computational domain is stretched in the normal direction to establish a 2.5D simulation model.A velocity boundary condition with a constant velocity of 100 m/s is set at the inlet of the wind tunnel，where the uniformly-distributed droplets have the same speed.The inlet air is chosen as an incompressible fluid at normal temperature，and the right outlet is at a constant pressure of 1.01×105Pa.The LWC at the wind tunnel inlet is 1 g/m3with the droplet diameter of 50 μm.The cases are calculated using three meshes of different densities.The number of nodes at mesh boundaries in Case 2 and Case 3 are 1.5 times and 2 times that in Case 1，respectively.And there are 43 263，97 893 and 174 523 mesh elements in each case.A schematic diagram of mesh for the icing wind tunnel is shown as Fig.5.

Fig.5 Schematic diagram of mesh for icing wind tunnel

Fig.6 compares the contour of the LWC for the icing wind tunnel obtained by the three meshes.After flowing parallel to the wind tunnel inlet，the airflow is deflected by the contraction of the wind tunnel.Due to the deflection of the airflow and the inertia of the water droplets，some water droplets hit the contraction wall，and the water droplets without impingement enter the test section with airflow.Water droplets converge along with the airflow towards the centerline of the test section，and the maximum of the LWC appears there.As shown in Fig.6，the distributions of the LWC obtained by the three meshes are not consistent.The maximum LWC obtained in Case 1 is 44.4 g/m3，while in Case 2 and Case 3 it is 55.9 g/m3and 73.8 g/m3，resulting in a large difference between the three cases.

Fig.6 Comparison of LWC obtained by three meshes for icing wind tunnel

Fig.7 shows the LWC comparison of multiple cross sections in the icing wind tunnel obtained by the three meshes.Atx=0.4 m，water droplets moving with the air flow slightly deflect in the contraction section，and the LWC distributions of the three meshes are in agreement.Atx=1.2 m andx=1.6 m，where the cross sections are inside the test section，the droplets have already deflected sharply with the airflow，the distance between the water droplets becomes much closer，and the LWC values of the three meshes appear to be different.Atx=2 m，where the cross section is at the outlet of the wind tunnel，the flow of water droplets is steadier than that inside the wind tunnel，and the difference between the three results decreases.In the simulation of the Eulerian method，they-directional droplet velocity component eventually zeros out at the centerline of the test section［21］.As the water droplets move backward，the LWC gradually increases along the centerline and remains larger near it.According to the fluid mechanics theory，the streamlines of a steady flow field are consistent with the motion trajectory and will never cross.And the droplet streamline in the Eulerian method is also its trajectory line.Since the droplet velocity in each control volume is a single value［8］，the obtained water droplet streamlines would never interact with each other［12］，which conflicts with the no-collision and no-coalesce assumption.

Fig.7 Comparison of LWC of multiple cross sections in the icing wind tunnel obtained by three meshes

It can be seen that the Eulerian method cannot capture the motion and impingement characteristics of water droplets when their motion is deflected by the upstream wall［21］.The velocity component in theydirection is handled differently from reality using the Eulerian method，and there is a numerical diffusion effect，both of which result in the deviation of water droplet motion and LWC distribution.Furthermore，the calculation of droplet motion and LWC distribution is affected by meshing，which leads to inconsistent results of the Eulerian method when using different meshes.The calculation results obtained by the Eulerian method are mesh-dependent.In the following cases，the mesh impact analysis is also carried out by comparing the results of the water droplet impingement characteristics with different meshes.

2.3 Comparison with an S-shape duct

The internal structure of the engine is complex，and once super-cooled water droplets enter the engine，they are affected by upstream components［22］.Therefore，the accuracy of the droplet impingement characteristic on its internal surfaces using the Eulerian method needs to be further analyzed.Due to the long flow path of the engine intake and its large effect on droplet motion，a 2D section of an S-shape duct is used to investigate the mesh impact of the Eulerian method.The 2.5D geometry model is shown in Fig.8，with a duct deflector cone at the rear position.The freestream velocity is set as 0.2Main the +xdirection，and the ambient pressure and temperature are 1.01×105Pa and 263.15 K，respectively.Both the outlet at the outside of the engine lip and the outlet of the internal flow field are set as constant pressure boundaries of 1.01×105Pa.The LWC at the inlet is set as 1 g/m3with the droplet diameter of 50 μm.Meshes with the same number of mesh elements and different boundary layer refinement methods are used to calculate the two cases.The mesh near the cone is shown in Fig.9.The mesh height of the boundary layer on the cone in Case 1 is 0.03 mm，while in Case 2 is 0.01 mm.

Fig.8 Schematic diagram of the LWC contour for S-shape duct

Fig.9 Schematic diagram of mesh near the cone

Fig.8 is a schematic diagram of the LWC contour for the S-shape duct.After the air and water droplets enter the duct，the movement of water droplets will deflect along with the airflow deflection because of the wall diversion.Due to the inertia of water droplets，water-shaded regions appear on both the upper and lower surfaces.Moreover，as a result of the bending direction of the duct，there is a larger area of the water-shaded region on the lower surface.With the backward flow of water droplets along with the airflow，the LWC around the engine cone is directly affected，and the droplet impingement characteristics on the cone surface are affected as well.

The comparison of the collection efficiencies between these meshes for the engine cone is plotted in Fig.10，wheresis the arc length of the cone with the upper and lower outlet as the beginning and end，respectively.As shown in Fig.10，the local collection efficiency reaches its maximum near the leading edge of the cone.Water droplets can hit the entire upper surface，while the lower surface has a smaller impingement range.The difference is found between the droplet impingement characteristics obtained by the two meshes，which indicates that the mesh refinement method also affects the results of the Eulerian method.The water droplet impingement characteristics of the S-shape duct obtained by the Eulerian method are mesh-dependent.The accuracy of the Eulerian method under this circumstance needs to be taken into account.

Fig.10 Comparison of local collection efficiencies obtained by two meshes for the engine cone

2.4 Comparison with a multi-element airfoil

A multi-element airfoil is used for the simulation to further examine the mesh impact of the Eulerian method for droplet impingement characteristics.The multi-element airfoil，as shown in Fig.11，is made up of a front slat，a main wing，and a flap.The computational domain is stretched in the normal direction to establish a 2.5D simulation model as well.The freestream velocity is set as 0.23Mawith an AOA of 4°，and the ambient pressure and temperature are 95 630 Pa and 278 K，respectively.The LWC at the far field is set as 1 g/m3with the droplet diameter of 21 μm.The cases are calculated using three meshes with different densities，and the number of mesh elements are 94 592，130 575 and 214 530，respectively.The number of nodes at mesh boundaries in Case 2 and Case 3 are 1.2 times and 1.5 times that in Case 1，respectively.The meshes of the three cases around the airfoil are shown in Fig.11.

Fig.11 Comparisons of overall meshes of three cases around the airfoil

Fig.12 is a schematic diagram of the LWC contour for the multi-element airfoil.Since the slat ahead blocks both the movement of the airflow and water droplet，droplets are deflected along with the airflow.As a result，the LWC around the airfoil is directly affected，as well as the droplet impingement characteristics on the airfoil.

Fig.12 Schematic diagram of contour of LWC for multi-element airfoil

Fig.13 compares the pressure coefficientsCpof the multi-element airfoil obtained by the three meshes，and the results show good agreement.The pressure coefficient distribution of the multi-element airfoil is independent of the number of mesh elements，which indicates the mesh independence of the airflow field.

Fig.13 Comparison of pressure coefficients obtained by three meshes for multi-element airfoil

The collection efficiencies obtained by the three meshes using the Eulerian method are compared for the multi-element airfoil，as shown in Fig.14.There is no frontal body to deflect the droplet before hitting the slat，and the droplet impingement area and local collection efficiency of the three meshes are consistent.In contrast，after the droplet movement is blocked by the slat，there is a large difference between the droplet impingement area and the local collection efficiency on the main wing of the three meshes.It can be seen that when there is no upstream effect，the droplet motion and impingement characteristics obtained by the Eulerian method are mesh-independent，and simulation results can be considered accurate.When there is an upstream effect，the results of the airflow field are in good agreement，proving that the airflow field is mesh-independent.However，due to the influence of trajectory crossing and numerical diffusion after droplet deflection，mesh has an impact on the calculation results of droplet motion and impingement characteristics obtained by the Eulerian method.Meanwhile，the results do not develop towards mesh independence with the increase of mesh density.Therefore，it indicates that the simulation results of the Eulerian method are mesh-dependent.The mesh impact on the accuracy of simulation using the Eulerian method should be considered in these cases.

Fig.14 Comparisons of collection efficiency obtained by three meshes for the multi-element airfoil

3 Conclusions

The Eulerian method is established and investigated for the calculation of water droplet impingement characteristics using different meshes.The impact of mesh on the results of the Eulerian method is analyzed，and then the accuracy of the Eulerian method with upstream effects is discussed.To validate the established Eulerian method and the computation settings，a case of a NACA 0012 airfoil is applied.The mesh impact of the Eulerian method is analyzed with the simulation of water droplet motion in an icing wind tunnel.And the accuracy of the Eulerian method for water droplet impingement characteristics is investigated by adopting the cases of an S-shape duct and a multi-element airfoil.The main conclusions are as follows：

（1）When the droplet motion is not affected by the upstream effect，different meshes achieve the almost identical droplet impingement characteristics to those in the literature for the NACA 0012 airfoil，validating the mesh independence of the Eulerian method without upstream droplet deflections.

（2）Water droplets are deflected by the wall diversion in the contraction section of the icing wind tunnel，and the prediction results of the Eulerian method after droplet deflection are mesh-dependent.The water droplet streamlines in the Eulerian method cannot cross，leading to the deviation of the Eulerian method in calculating the droplet motion after deflection.Moreover，mesh has an impact on droplet motion calculations，which results in the difference between water droplet motion calculated by different meshes.

（3）When the droplet motion is deflected by upstream components，due to the effects of trajectory crossing and numerical diffusion after droplet deflection，there is a difference between the droplet impingement characteristics of various meshes.And the prediction results do not develop towards mesh independence with the increase of mesh density，while the airflow filed results have already been mesh-independent.With the upstream effect，the droplet impingement characteristics obtained by the Eulerian method are mesh-dependent.The impact of mesh on the accuracy of the Eulerian method should to be considered in future studies.