Numerical investigation of the influences of ridge ice parameters on lift and drag coefficients of airfoils through design of experiments

The complexity of ice accretion shapes on aircraft brings challenges in establishing the physical correlations between the geometric parameters of ice and the aerodynamic performances of an airfoil. The present study applied a quarter-circle simplified ridge ice model and designed four parameters, including the lower ice horn, the height of the upper ice horn, the angle of the upper ice horn, and the location of the upper ice horn. Design of Experiments (DOE) was adopted to design the simulation process and analyze the obtained data. Numerical investigations were carried out to study the effects of these four design parameters on the lift coefficients (CL) and drag coefficients (CD) of three airfoils at various angles of attack. The results revealed that the four design parameters had influences on CL and CD in all simulation cases. The influence of the height of the upper ice horn was always dominant, followed by the angle and location of the upper ice horn. The influence degrees of the lower ice horn varied with the angles of attack and airfoils. The change trends of CL and CD were opposite while changing the value of any parameter in most simulation cases. In addition, the leading-edge radii of these three airfoils, as well as the shapes of the upper and lower surfaces, could alter the influences of design parameters on CL and CD.


Introduction
Icing may occur on the windward side of the aircraft and engine, including the wing, windshield, rotor, tail rotor, engine intake device, instrument sensor, etc. [1][2][3][4] A large number of tests have proved that a small amount of icing on the key parts of an aircraft can cause a considerable decrease in lift, an increase in drag, and a deterioration in maneuverability and stability. 5Besides, the ice accretion in the intake port and power unit of the aircraft engine may induce fatal consequences such as engine damage and shutdown.At the same time, ice shedding might lead to internal engine damage. 6In 1994, the famous Roselawn air crash in Indiana of USA revealed the existence of supercooled large droplet (SLD) icing condition. 7,8This accident brought the ridge ice to the attentions of the researchers.When the heated leading-edge ice-protected surface was not running at 100% evaporation, water flowed back on the surface from the heated section to freeze downstream on the cooler unheated surface, forming the ridge ice.Based on the ice geometries, the shapes of the ice can be divided into ridge ice, 9 horn ice 10 (usually produced in glaze ice conditions and had some horn-shaped protrusions), streamwise ice 11 (usually formed in rime icing conditions and appeared to be an extension of the airfoils), and so on.The ridge ice was found to cause a more striking aerodynamic deterioration than other ice shapes (e.g., horn ice or streamwise ice).
Considering the ice-tolerant design of the control surfaces of aircraft is a useful way to reduce the adverse effects of icing.Generally, a robust aerodynamic design of airfoil under ice accretion scenarios needs a better insight into the correlations between performance losses and icing conditions.There are many advantages of using simplified ice shapes to study the aerodynamic performance of airfoils.For example, it simplifies geometric modeling of the ice shape and lowers aircraft development costs.Therefore, it is of great research value and practical significance to develop a simplified ice shape that has similar influences on the airfoil as the original ice shape.
So far, a series of studies on the aerodynamic performances of ridge ice have been reported, which include flow field investigations and geometric simplifications.Miller and Addy 12,13 recorded the icing photos of the ridge ice on the model with a pneumatic deicing system in the SLD icing test in the icing tunnel.The front area of the ridge ice was found to be relatively clean, with a certain inhomogeneity in the ridge direction.Lee and Bragg 14 studied the aerodynamic effects of simulated SLD ice accretion on a modified NACA 23012 airfoil with forward-facing quarter-round simulations.The height-to-chord ratio of the airfoil ranged from 0.0083 to 0.0139.They found that the simulated ice shape produced a large separation bubble as the real ice did.Moreover, Lee and Bragg 15 found that the variations in the simulated ice shape geometry had a considerable effect on the airfoil aerodynamics, with the half-round shape exhibiting a significantly higher lift coefficient.Dunn and Loth 16 predicted sectional aerodynamic characteristics for the quarter-round ice shape.Significant reductions in the lift were obtained for those relatively small protuberances, which was consistent with experimental findings.The results also showed a good predictive performance on drag, pitchingmoment, and hinge-moment variations.Bragg and Loth 17 studied the different effects of ridge ice accretion on airfoil aerodynamics between airfoils NACA 23012m and NLF 0414.The results showed that the protuberance of ice affected the airfoil NLF 0414 very differently from the airfoil NACA 23012m, and the lift coefficient of the half-round shape was closer to that of the real ice.Lee, 18 Kumar and Loth 19 studied the airfoil lift, drag, pitching moment, and hinge moment with different chordwise locations, sizes, and shapes of the ridge ice simulations.By utilizing the Particle Image Velocimetry (PIV) technique, Gurbacki 20 found that the ridge ice and horn ice had a similar flow pattern.Lee and Bragg 21 found that increasing the iceshape height generally resulted in more severe performance degradation.The exception was for the ice located at the leading edge of the airfoil.Meanwhile, the ice height no longer significantly degraded the aerodynamic performance once a critical height was reached.
In terms of numerical simulations, Pan and Loth 22 used the detached eddy simulation to improve the predictive ability of the stall behavior for an iced airfoil.The validation results were in good agreement with the experimental and previous numerical results concerning time-averaged aerodynamic forces as well as the unsteady flow pattern.Broeren et al. 23 presented the results of recent investigations on the aerodynamics of simulated runback ice accretion on airfoils and recorded the detailed geometrical shape and aerodynamic parameters of this ridge ice.In addition, Broeren et al. 24 reported the results of experimental studies which established a high-fidelity, full-scale, iced-airfoil aerodynamic performance database.Based on experimental data and previous studies, Li 25 numerically modeled the ridge ice of airfoil NACA23012.By using a forward quarter circle model, the aerodynamic changes caused by various ridge ice were studied.Pouryoussefi et al. 26 numerically investigated the aerodynamic behavior of three ice shapes proposed by Broeren et al. 23 They found that ridge ice had the most significant effects on the aerodynamics of an airfoil.Zhang et al. 27 numerically studied the relationship between the lift coefficient and flow field parameters.The results demonstrated that the zonal detached-eddy simulation model exhibited additional advantages compared to the Spalart-Allmaras turbulence model.Costes and Moens 28 further found that both the Spalart-Allmaras model and only the RANS model could effectively converge to a steady state, regardless of the flow conditions.Li et al. 29 developed a threeequation turbulence model, which proved to apply to a wide range of ice shapes.Mortensen et al. 30 compared different turbulence models by implementing the leading-edge icing calculation.The results showed that the S-A model worked best in coping with the aerodynamic force prediction.
The previous studies mainly focused on simplifying ridge ice into a forward quarter circle shape because the aerodynamics of this ridge ice was closer to the real ice.8][29][30] However, previous studies mainly considered the effects of the height and the location of the ice horn, [14][15][16][17][18][19][20][21][22][23][24][25][26] but ignored the lower ice horn.Moreover, few studies considered the effects of the same ridge ice shape on different airfoils.Thus, the present study numerically investigated the lift coefficients (C L ) and drag coefficients (C D ) of iced airfoils based on the forward quarter circle ridge ice with the Spalart-Allmaras model.In addition, the Design of the experiments (DOE) was adopted for the first time to design the simulation cases and analyze the obtained C L and C D .The influences of four design parameters (i.e., the height of the upper ice horn, the angle of the upper ice horn, the location of the upper ice horn, and the lower ice horn) of a ridge ice model on three airfoils (i.e., NACA23012, NACA64-215, and RAE2822) were carefully analyzed.The present study aimed to improve our understanding of the influences of ridge ice parameters on airfoils and propose some guidelines on ice shape simplification.

Geometric model of the ridge ice
Based on previous studies, 14,23,30 we investigated the quarter-circle simplified ridge ice model and defined four design parameters, namely the lower ice horn, the height of the upper ice horn (h), the angle of the upper ice horn (b), and the location of the upper ice horn.The 90º angle of the upper ice horn meant that the straight edge of the model was normal to the upper surface of the airfoil.The location of the upper ice horn was defined as the distance from the ice horn to the leading edge along the chordal direction, which was usually in units of chord length (c).The shapes of the lower ice horn remained unchanged in different cases.Besides, there were only two options for the lower ice horn: appearance and disappearance.

Numerical method
In this study, numerical investigations were carried out by the CFD tool, ANSYS-Fluent.The Spalart-Allmaras model was chosen according to the study of Mortensen et al. 30 In order to verify the accuracy of the numerical method, Broeren's EG1159 ridge ice shape 23 on the airfoil NACA 23012 was simulated in this section.The chord length was 3.6 m and the diameter of the calculation domain was set at 20 times the chord length.Structured grids were adopted in this study.The thickness of the first-layer grids near the wall was set as 6 3 10 À6 m, which ensured the y + approaching 1.The calculation information is shown in Table 1.

Validation of numerical investigation
Figure 2 shows the comparisons between the simulation and experimental results 23 with respect to the C L with different numbers of grids at AOA = 6°.The parameters of the simulation were consistent with those in the experiment described in Table 1.As the refinement of the grids, the discrepancy between the simulation and the experimental results converged to the order of 1 3 10 À4 .Considering the computational resources, the grid number of 1:8 3 10 5 was selected for the present numerical investigations.The grids around the whole airfoil and the upper ice horn were shown in Figure 3.
Figure 4 shows the comparisons of the present simulations and the previous experimental results 23 at various angles of attack.It can be seen that the simulation curves had a similar regular pattern as the experiment curves, and the differences between these values were very small.

Data analysis method
Design of experiments (DOE) is a structured and organized method to determine relationships between the input and output of a process. 31The main objective of DOE is to obtain the maximum amount of information while limiting the number of experiments needed. 32The DOE method specifically used is factorial design, which is a multi-factor and cross-classification design. 33actorial design enables an experimenter to study the effect of the design parameters on the objective coefficient.
The values for three design parameters are shown in Table 2.It should be noted that the parameter of the lower ice horn did not have a specific value.Instead, it only had two options: appearance and disappearance.
The analysis process with respect to geometric parameters of the ridge ice and airfoils was organized as follows.Firstly, a permutation and combination table of four parameters was designed (Table 3).The influences of four design parameters on C L and C D of a specific airfoil at a specific angle of attack required 20 simulations.The present study investigated the lift and drag coefficients of airfoil NACA23012 at three angles of attack (i.e., AOA = 2°, 5°, and 8°), as well as the airfoils NACA64-215 and RAE2822 at AOA = 2°.It should be noted that AOA = 2°, 5°, and 8°represented the cruise condition, the stall condition, and the poststall condition, respectively.The total number of simulation cases was 100.
Secondly, the airfoil models and grids were generated respectively based on the permutation and combination table.Thirdly, the C L and C D of each model were simulated.Finally, the reference value and the standardized effects of each parameter were figured out.In the present study, the reference value was the 0.95 quantile of a t-distribution with degrees of freedom equal to the degrees of freedom for the error term.The standardized effects were calculated by some formula. 31aking 2 3 2 factorial design 31,[34][35][36] as an example, Table 4 shows the relationship between the lift coefficients of the airfoil and two ice shape parameters.Each parameter had two values, which were denoted as A1, A2, B1, and B2.Due to the influences of these two parameters, the lift coefficients of the airfoil were A1B1, A1B2, A2B1, and A2B2.The standardized effect of the width of the upper ice horn can be calculated by equation (1).
The standardized effect of the height of the upper ice horn can be calculated by equation (2).
If the standardized effect of a parameter was larger than the reference value, it meant that this parameter was statistically significant at 0.1 level with the current model terms.In other words, this parameter had an influence on C L and C D .A larger standardized effect indicated that the influence degree of this parameter was greater.Meanwhile, the influence trends of the C L and C D with a certain ridge ice parameter were analyzed.The normalized effect (n) was defined as while s is the standardized effect and r is the reference value in the same case.Through this method, the differences in standardized effects in different cases were narrowed, and the simulation results of different cases could be plotted in a single figure.Results and discussions

C L and C D of airfoil NACA 23012 at different angles of attack
In this study, three angles of attack (i.e., AOA = 2°, 5°, and 8°) were selected for the analysis.It should be noted that AOA = 2°, 5°, and 8°represented the cruise, stall, and post-stall conditions, respectively.Figure 5 shows the influence degrees of the four design parameters on the C L of the airfoil NACA23012 at different angles of attack.The X coordinate represents the parameters, while the Y coordinate represents the normalized effect.
The normalized effects of these design parameters were all greater than 1, which meant that they all had an influence on C L .At AOA = 2°and 8°, the influence degrees of design parameters on C L were in descending order: the lower ice horn, the height of the upper ice horn, the location of the upper ice horn, and the angle of the upper ice horn.At AOA = 5°, the influence degrees of design parameters on C L were in descending order: the height of the upper ice horn, the location of the upper ice horn, the angle of the upper ice horn, and the lower ice horn.It can be found that the influence degrees of the height of the upper ice horn on C L were always dominant.The influence degrees of the angle and location of the upper ice horn on C L were weaker than that of the height of the upper ice horn.The influence degrees of the lower ice horn on C L were significant at 2°and 8°, but relatively weak at 5°.It should be noted that there was no relationship between the normalized effects of these three angles of attack.For the convenience of illustration, the data of these three angles of attack were plotted in one figure hereafter.Figure 6 shows the influence degrees of the four design parameters on the C D of the airfoil NACA23012 at different angles of attack.It can be seen that the normalized effects of four design parameters were all greater than 1, which meant that they had an influence on C D .At AOA = 2°, the influence degrees of design parameters on C D were in descending order: the lower ice horn, the height of the upper ice horn, the angle of the upper ice horn, and the location of the upper ice horn (Figure 6).At AOA = 5°and 8°, the influence degrees of design parameters on C D were in descending order: the height of the upper ice horn, the lower ice horn, the location of the upper ice horn, and the angle of the upper ice horn (Figure 6).It can be found that the influence degrees of the height of the upper ice horn on C D were always dominant.The influence degrees of the angle and location of the upper ice horn on C D were weaker than that of the height of the upper ice horn.These findings were consistent with those in Figure 5.It meant that these three design parameters had similar influence degrees on C L and C D .In addition, the lower ice horn had a significant influence on C D at these three angles of attack.
The influence degrees of four design parameters on C L and C D were obtained through the previous analysis.Further analysis was carried out in this study to get the influence trends of design parameters on C L and C D . Figure 7 shows the influence trends of the four design parameters on Cl and Cd of airfoil NACA23012 at AOA = 2°, 5°, and 8°.Considering the large difference between C L and C D at these three angles of attack, the values of C L and C D with specified parameters (i.e., 80°of the angle of the upper ice horn, 0.0855 c of the location of the upper ice horn, 0.0117 c of the height of the upper ice horn, and the existence of the lower ice horn) were selected as the benchmark.In order to better describe the influence trend, a term, namely the changes of coefficients (D), was defined as  It should be noted that the detailed values of C L and C D were not indicated in Figure 7. Table 5 shows these values as a supplement.The slopes of lines in Figure 7 can be calculated from the data listed in Table 5.These slopes can be used to derive C L and C D for other values of parameters.For example, when the angle of the upper ice horn is 95°and the angle of attack is 2°, the C L of airfoil NACA23012 can be estimated as 0.25253.

C L and C D of three airfoils at AOA = 2°I
n this section, the C L and C D of three airfoils (i.e., NACA 23012, NACA 64-215, and RAE 2822) at AOA = 2°were simulated and compared in detail.Figure 8 shows the influence degrees of the four design parameters on C L of these three airfoils.The normalized effects of design parameters were all greater than 1, which meant that they all had an influence on the C L .The influence degrees of design parameters on C L of airfoil NACA23012 were in descending order: the lower ice horn, the height of the upper ice horn, the location of the upper ice horn, and the angle of the upper ice horn.The influence degrees of design parameters on C L of airfoil NACA64-215 were in descending order: the lower ice horn, the height of the upper ice horn, the angle of the upper ice horn, and the location of the upper ice horn.The influence degrees of design parameters on C L of airfoil RAE2822 were in descending order: the height of the upper ice horn, the lower ice horn, the angle of the upper ice horn, and the location of the upper ice horn.It can be inferred that the influence degrees of the height of the upper ice horn and the lower ice horn on C L were dominant, while the other two parameters had relatively less influence.The influence degrees of the same parameter were different for these three airfoils.This might be caused by the leading-edge geometries of different airfoils.It should be noted that there was no relationship between the normalized effect of these three airfoils.For the convenience of illustration, the data of these three airfoils were plotted in one figure hereafter.
Figure 9 shows the influence degrees of the four design parameters on the C D of three airfoils.It can be seen that the normalized effects of four design parameters were all greater than 1.Each design parameter had an influence on the C D .The influence degrees of four design parameters on C D were in descending order: the lower ice horn, the height of the upper ice horn, the angle of the upper ice horn, and the location of the upper ice horn.Besides, the influence degrees of the height of the upper ice horn and the lower ice horn on C D were always dominant, while the other two parameters had less influence.These findings were consistent with those in Figure 8, which meant that these four design parameters had similar influence degrees on C L and C D .
Figure 10 shows the influence trends of four design parameters on three airfoils at AOA = 2°.It can be seen that C L increased and C D decreased as the angle of the upper ice horn increased or the height of the upper ice horn decreased (Figure 10    Some phenomena can be explained by the geometric structure of the airfoils.The leading-edge radius of the NACA23012 airfoil was relatively big in comparison to the other airfoils, with an evident protrusion on the upper surface of the airfoil.Thus, its suction peak was close to the leading edge of the airfoil.The presence of the upper ice horn caused a significant decline in the pressure distribution curve on the upper surface of the airfoil.And this decline was delayed as the upper ice horn moved towards the trailing edge of the airfoil.This could increase the C L of the airfoil, as shown in Figure 10(b).The upper surface of RAE2822 was relatively flat and the maximum thickness point was closer to the mid-section of the airfoil.As a result, the C L of airfoil RAE2822 were less impacted by the location of the upper ice horn throughout the parameter ranges of this study.The increased C L of the airfoil RAE2822 after removing the lower ice horn in Figure 10(d) might be due to the convex lower surface of the airfoil, while the lower surfaces of the other two airfoils were flatter.In addition, the subtle differences between the three airfoils in Figure 10(a) and (c) might be due to the different leading-edge radii of the airfoils.
Table 6 also lists detailed values of C L and C D of these three airfoils as a supplement.

Conclusion
The present work mainly focused on the simplified ridge ice parameters and investigated the qualitative and quantitative changes in C L and C D of airfoils caused by changes in these parameters.To the authors' best knowledge, this study adopted the DOE method for the first time to analyze the influences of ridge ice parameters on the lift and drag coefficients.Four design parameters (i.e., the height of the upper ice horn, the angle of the upper ice horn, the location of the upper ice horn, and the lower ice horn) were selected for analysis.Numerical investigations were carried out to study the influences of four design parameters on the C L and C D of three airfoils (i.e., airfoil NACA23012, airfoil NACA64-215, and RAE2822) at three angles of attack.Since there were no apparent interactions between the design parameters within the parameter ranges of this study, the present work mainly focused on the individual effects of design parameters.The following conclusions were obtained: 1.For airfoil NACA23012, the influence degrees of the height of the upper ice horn on C L and C D were always dominant.The influence degrees of the angle and location of the upper ice horn were weaker than that of the height of the upper ice horn.The C L increased and C D decreased when the angle of the upper ice horn increased, the location of the upper ice horn moved toward the trailing edge of the airfoil, or the height of the upper ice horn decreased.These conclusions were applicable to the cases of three angles of attack.2. For airfoil NACA23012, the influence of the lower ice horn on C L and C D was more complex.As far as the influence degrees were concerned, the lower ice horn had less influence on C L at AOA = 5°, but had significant influences in other cases.As far as the influence trends were concerned, removing the lower ice horn increased C L and decreased C D at AOA = 5°and 8°, but decreased C L and C D simultaneously at AOA = 2°.The present work provides some valuable insights into the ice shape simplification methods, which could not only improve our understanding of the simplification methods of ridge ice, but also provide a reference for researchers to carry out further experimental and numerical analyses.A larger range of design parameters is recommended to be investigated in future studies to establish more comprehensive relationships among design parameters, airfoils, and angles of attack.

Authors' contributions
Chengyi Zheng contributed to the investigation, methodology, simulations, data analysis, and writing-original draft;

Figure 3 .
Figure 3. Grids.(a) The whole airfoil and (b) The upper ice horn.

Figure 4 .
Figure 4. Comparisons between the simulation results and the experiment data.(a) Comparisons of the C L and (b) Comparisons of the C D.

Figure 6 .
Figure 6.Influence degrees on the C D of airfoil NACA23012.

Figure 5 .
Figure 5. Influence degrees on the C L of airfoil NACA23012.

Figure 7 .
Figure 7. Influence trends on the C L and C D of airfoil NACA23012.(a) The angle of the upper ice horn, (b) The location of the upper ice horn, (c) The height of the upper ice horn, and (d) The lower ice horn.
(a) and (c)).The C L of these three airfoils and C D of airfoil RAE2822 increased

Figure 8 .
Figure 8. Influence degrees on C L of three airfoils at AOA = 2°.Figure 9. Influence degrees on C D of three airfoils at AOA = 2°.

Figure 9 .
Figure 8. Influence degrees on C L of three airfoils at AOA = 2°.Figure 9. Influence degrees on C D of three airfoils at AOA = 2°.

Figure 10 .
Figure 10.Influence trends on C L and C D of three airfoils at AOA = 2°.(a) The angle of the upper ice horn, (b) The location of the upper ice horn, (c) The height of the upper ice horn, and (d) The lower ice horn.

3 .
For different airfoils at AOA = 2°, the height of the upper ice horn and the lower ice horn always had significant influences on C L and C D , while the other two parameters had relatively less influence.In most simulation cases, changing the values of parameters could increase C L and decrease C D .However, C L and C D increased simultaneously when the location of the upper ice horn moved towards the tailing edge of airfoil RAE2822.In addition, C L and C D decreased simultaneously when the lower ice horn of airfoils NACA23012 and NACA64-215 were removed.4. Changes in the height and angle of the upper ice horn exhibited comparable effects on the C L and C D of various airfoils.The location of the upper ice horn had less influence on the airfoil with a flat upper surface (RAE2822).The presence of a lower ice horn increased the C L and C D of airfoils with a flat lower surface (NACA23012 and NACA64-215).The presence of the lower ice horn decreased the C L of the airfoil RAE2822 since the lower surface of the airfoil RAE2822 had a protrusion near to the leading edge.

Table 1 .
Information on ridge ice.

Table 3 .
Permutation and combination table of four design parameters.

Table 2 .
Design parameters and their values.

Table 5 .
C L and C D of NACA23012 at three angles of attack.

Table 6 .
C L and C D of these three airfoils at AOA = 2°.Qiaotian Dong was involved in methodology; Zhigang Yang was involved in project administration; Zheyan Jin contributed to supervision, methodology, writing, review, and editing, project administration, funding acquisition.