Performance Deterioration of Pitot Tubes Caused by In-Flight Ice Accretion: A Numerical Investigation

In-flight ice accretion on typical pitot-static systems is numerically investigated to reveal their performance deterioration under both rime and glaze icing. Coupled with the open source computational fluid dynamics (CFD) platform, OpenFOAM, the numerical strategy integrates the airflow determination by the Reynolds-averaged Navier-Stokes equations, droplet collection evaluation by Eulerian representation, and ice accumulation by mass and energy conservation. Under varying inflow conditions and wall temperatures, the calculated ice accretion performance indicates that the ambient temperature has the most significant effect on the icing-induced failure time, leading to an almost exponential growth. Meanwhile, the blocking time is found to be linearly proportional to the increase in wall temperature. With the increase in inflow velocity, the failure time follows a parabolic variation with glaze ice accretion while shows a monotonic reduction under rime icing conditions. In addition, when the angle of attack increases, failure accelerates under both the glaze and rime icing scenarios. These findings provide guidance for the protection design of pitot tubes. A nonlinear regression analysis is further conducted to estimate the failure performance. The predicated failure times show reliable consistency with numerical results, demonstrating the capability of the obtained empirical functions for convenient predictions of failure times within the applicable range.


Introduction
In-flight ice accretion is an extremely severe hazard to aviation safety. When passing through clouds containing supercooled droplets, aircraft/helicopters/unmanned aerial vehicles (UAVs) suffer from ice accumulation on exposed upwind parts, degrading flight performance and even leading to flight accidents and fatalities [1]. In 2009, Air France Flight 447 from Brazil to France crashed into the Atlantic Ocean, and all 228 people on board were killed. The final report on the incident announced that the tragedy was induced by ice accretion on the pitot tubes. The functioning of the pitot-system of the aircraft was hampered by ice crystals, leading to temporary inconsistencies between airspeed measurements. Therefore, the autopilot was disconnected, and the crew had incorrect responses; this eventually resulted in a nonrecoverable aerodynamic stall and control failure [2].
In general, a variety of complicated factors dominate inflight ice accumulation processes, such as meteorological conditions (temperature, liquid water content in the air, size and distribution of the super-cooled water droplets in clouds), navigation status of aircraft (flight velocity, angle of attack (AoA), and altitude), and aircraft instantaneous profiles (especially, shapes of wing leading edges) [3]. As a result of these factors (and differences therein), in-flight ice formation presents various appearances; however, it is generally distinguished into two typical types: rime icing and glaze (or clear) icing. When air temperature falls below −10°C, rime ice instantaneously forms when small supercooled droplets collide with aircraft surfaces. The attached ice appears opaque with a white texture and a rough coating. Glaze icing occurs under the combined conditions of a comparatively high temperature near the freezing point (−10°C to 0°C) and the presence of larger supercooled water droplets in clouds. Driven by aerodynamic forces, water droplets spread and form a thin liquid film (of the order of tens of microns), run back some distance along aircraft surfaces, and then gradually freeze downstream. Hence, the irregular "double horn" shapes (protrusions) can be visibly observed on the upper and lower sides of the leading edges of aircraft.
Thus far, in-flight ice accretion remains a severely pressing issue that needs to be thoroughly investigated and overcome. Intensive research, especially in the field of numerical simulations, is ongoing in this regard. Numerical approaches on ice accretion are capable of addressing a wide variety of combinations of test model configurations, flight conditions, and icing occurrence conditions. Currently, LEWICE 2D&3D (developed by NASA Lewis Research Center, the United States) and FENSAP-ICE (merged by ANSYS) are accepted as the most representative ice accretion simulating software. From the 1980s, NASA commenced the development of the LEWICE code. Along with continuous improvements and updates thus far [4][5][6][7][8], LEWICE has been proven to be robust and adequately accurate based on the comparison of its predicted results with a variety of experimental outcomes. LEWICE adopts the three-dimensional (3D) panel-method with potential flow equations to estimate airflow fields, utilizes the Lagrangian approach with droplet motion equations to track water trajectories, and the Messinger model to compute ice accretion [9]. Other codes based upon Lagrangian droplet motion models are very similar to LEWICE, such as ONERA [10] and CANICE [11]. Depending on the theory of Eulerian two-phase flow, Bourgault et al. [12,13] initially proposed an Eulerian approach to numerically model droplet impingement with different stabilization terms. On this basis, the commercial code FENSAP-ICE was put forth for 3D ice shape (rime or glaze) predictions. The code uses a shallow-water thermodynamic icing model (SWIM), developed by combining the classical Messinger thermodynamic model and a partial differential equation (PDE) system [14]. Building further on this foundation, FENSAP-ICE was subsequently enhanced with several advanced models, such as models for ice crystals [15] and complex geometries [16], deicing model [17], and supercooled large droplet (SLD) splashing and bouncing model [18].
In addition to improving numerical models pertaining to ice accretion [19][20][21], studies also focused on broadening their engineering applications [22][23][24]. In contrast to the numerous icing studies on airfoils, engine inlets, antennas, and other aircraft components [25,26], icing problems of pitot-static systems have been rarely investigated. The existing numerical studies on this subject either have very simplified geometries [27] or focus more on anti-/deicing issues [28,29]. Currently, pitot tubes are being increasingly implemented in light airplanes, and even UAVs; this poses a serious risk due to in-flight ice accretion.
Given this background, the objective of the current study is to systematically reveal the performance degradation of pitot-static systems under different in-flight icing conditions. The dedicated numerical model, for efficient and accurate prediction of in-flight ice accretion, is introduced in Section 2. The effectiveness of the model is validated in Section 3 by comparing the calculated results with previous experimental and numerical data from initial airflow fields to final rime and glaze ice shapes. The geometric configurations and physical models are described in Section 4. With the validated numerical approach, Section 5 addresses the issue of inflight ice accretion on the pitot-static system with particular emphasis on the influence of varying flight conditions on this system. A regression analysis is further conducted to reveal the performance degradation of the pitot tubes under rime and glaze ice accretion. Finally, conclusions from the study are detailed in Section 6.

Numerical Methods
The proposed numerical approach is developed on the C++ finite-volume open-source platform-OpenFOAM® [30]. To attain a balance between calculation efficiency and accuracy, the method undertaken yields the quasi-steady hypothesis and assumes the one-way interaction by only considering the effects of airflow on droplets, whereas neglecting the reverse impact of droplets on airflow. Figure 1 illustrates the modular iterative logic of the numerical strategy. The integrated submodules are solved independently, and only several correlative parameters are transferred among these submodels. The submodules of airflow field determination, droplet impingement, and ice accumulation are elaborated in detail in the following subsections.

Airflow Field.
Employing the rhoPimpleFOAM solver of OpenFOAM [30], the aerodynamic flow field is determined by solving the 3D compressible Reynolds-averaged Navier-Stokes (RANS) equations for the conservation of mass, momentum, and energy, as follows: Continuity equation Momentum equation Energy equation where ρ is the air density, U is the air velocity, P is the pressure, μ is the viscosity, e is the energy per unit volume, and λ is the thermal conductivity. One-equation Spalart-Allmaras (SA) turbulence model [31] is introduced as a supplementary to close the partial differential equation (PDE) system. Since the original SA turbulence model is derived for smooth surfaces, the convective heat flux (one of predominant factors determining ice growth) is dramatically underestimated. Therefore, roughness effects are implemented into the original SA model to accommodate rough frozen surfaces. This extension includes    3 International Journal of Aerospace Engineering the rough wall treatment as stated in the work of Aupoix and Spalart [32] and the equivalent sand-grain roughness model following the empirical representation summarized by NASA LEWCIE [9]. Therefore, the transport equation for the SA model is modified as where d = d min + Ks signifies the nearest distance from walls (d min ) associated with the equivalent sand-grain roughness (Ks). The term f w is given as follows: f w = g ½ð1 + C 6 w3 Þ/ðg 2 + C 6 w3 Þ 1/6 , g = r + C w2 ðr 6 − rÞ, r ≡ν/Sκ 2 d 2 .
The temporal and spatial discretization of the solution domain is implemented via the finite volume method. The

Droplet Impingement.
The Eulerian formulation is proposed to evaluate water droplet trajectories and impinging characteristics. Compared to the traditional Lagrangian approach, Eulerian method is more practical and costeffective for calculating 3D icing. The droplet-phase governing equations are numerically solved in the same set of grids as the airflow solver and derived as a continuity equation and a momentum equation, which read as follows: where α denotes the nondimensional local volume fraction of droplets to air phase and U p is the droplet velocity vector. The source term of drag components induced by the ambient airflow is formulated byF = αð f ðRe r Þ/τ p ÞðU − U p Þ, where Re r = jU − U p jρ p d p /μ characterizes the Reynolds number of droplets, d p is the droplet diameter and refers to the mean volume diameter (MVD) in this research, the momentum response timeτ p = ρ p d 2 p /18μ is an inertial parameter to describe variations of droplet velocities, and f ðRe r Þ is a nonlinear empirical formula of drag coefficient, estimated by f ðRe r Þ = 1 + 0:15 Re 0:687 r + ð0:0175 Re r /ð1 + 45000 Re −1:16 r ÞÞ [34]. Taking formality of the Eulerian averaged approach, the nondimensional local collection efficiency is expressed as where U p∞ is the velocity of freestream, n is the surface normal vector, LWC stands for liquid water content, and _ M imp and _ M ∞ are the impinging and freestream water mass flux, respectively. Accordingly, droplet impinging mass fluxes are determined by _ M imp = β ⋅ _ M ∞ for subsequent calculations regarding ice accumulation.  Model. The thermodynamic model for ice accretion is developed on the basis of the classical Messiger model [35]. As illustrated in Figure 2, depending on the concept of the conservation law, the mass and energy equilibrium formulations are established as follows: The incoming mass flux is a summation of the entering run-back water mass flux from the upstream control volumes . _ M evap is estimated as a correlation among the saturation vapor pressure, convective heat transfer coefficient, and equilibrium temperature, proposed by Macarthur et al. [36] and utilized by FENSAP-ICE code [14,37].
Determination of _ M in and _ M out depends on the air shear stress [26]. As diagrammatized in Figure 3, the corresponding mass fluxes through the flow-out edge 3 and edge 4 are ascertained via where l 3 and l 4 denote the lengths of edges 3 and 4, respectively, τ is the shear stress obtained from the airflow solver, and n 3 and n 4 are the edge normal vectors of edges 3 and 4, respectively.
The gained heat fluxes are contributed by the kinetic energy and the enthalpy carried by the impacting droplets ( _ Q imp ), enthalpy of run-back water from upstream control volumes ( _ Q in ), and latent heat released by accumulated ice ( _ Q ice ). _ Q imp and _ Q in are determined as follows: where T ∞ is the freestream temperature, T 0 is the freezing temperature, and T in,n is the temperature of the upstream control volume, which has a common edge n with the current control volume. The heat loss comprises of the enthalpy carried away by the leaving run-back water to the downstream control where h c is the convective coefficient, ∂T/∂n is the normal gradient of temperature distribution to the aerodynamic surface, and T w,ini is the initial temperature of the wall. The above conservative formulas include three unknowns: the mass flux of accreted ice _ M ice , the flow-out run-back water mass flux _ M out , and the equilibrium temperature of the aerodynamic surface T s . The freezing fraction f is therefore introduced to circumvent the inadequacy due to these unknowns and is defined as f ≡ ðfreezing mass/ incoming massÞ = ð _ M ice /ð _ M imp + _ M in ÞÞ: Depending on the equilibrium temperature, the compatibility relations are different. If the equilibrium temperature is above the freezing point, no ice formation occurs in the current control volume, and all of the liquid water flows out to the adjacent down-stream control volumes. The compatibility relation can be given as follows: If the equilibrium temperature is the freezing point, the ice and water are coexistent in the current control volume, and the compatibility relation becomes If the equilibrium temperature is below the freezing point, all captured droplets freeze with no water remaining   International Journal of Aerospace Engineering in the current control volume, and the relation is changed as given below: Here, C p,ice is the specific heat of ice, L f is the fusion latent heat of water, and L evap and L sub are the latent heat with respect to the freezing point for vaporization and sublimation, respectively.

Solving Strategy.
With the assumption of no flow-in run-back water in control volumes, _ M in = 0, the flowchart of the ice accretion solver is elaborated in Figure 4.
After the thermodynamic iteration, the ice thickness is finalized by whereΔt is the time step and ρ ice is the ice density. Each grid node is looped to achieve the procured ice thickness along the normal vector direction of ice-forming surfaces. An ideal position is predefined by the mesh deformation algorithm in OpenFOAM to reconstruct the ice shape. The grid system is regenerated when the mesh quality is too poor for the simulation to be continued.

Validation of the Numerical Approach
In this section, the proposed icing code is thoroughly validated part by part-the airflow fields, the droplet impingement, the rime ice accretion, and the glaze ice accretion. The predicted results are compared with the respective numerical results or the experimental data reported in the literature.

Validation of the Airflow Field Solver.
The airflow field is validated by comparing the current surface pressure coefficient around a two-dimensional (2D) NACA0012 airfoil with the experimental results by NASA [38]. Inlet conditions are listed in Table 1. As depicted in Figure 5, the current result is consistent with the experimental data obtained by Ladson et al. [38].

3.2.
Validation of the Droplet Impingement Solver. The 2D test case for the proposed droplet impinging solver is performed on a NACA0012 airfoil with a MVD of 20 μm and a LWC of 0.55 g/m 3 . Details of simulation conditions are presented in Table 2. Figure 6 compares the local collection efficiency calculated by the current code with that by LEWICE and FENSAP-ICE [39]. The predicted curve shape exhibits good consistency with that of LWWICE. The slight disparity in the current impinging rate with that observed in FENSAP-ICE may be caused by the differences in the predicated airflows.  Table 3.   Figure 7 presents a comparison of the predicted ice accumulation within 7 min based on the IRT experimental, FENASP-ICE, and LEWICE data. The current and FENSAP-ICE results are observed to be in good agreement. Compared to the experimental measurements, LEWICE well resolves the lower icing bump, while the current result exhibits a good agreement at the upper leading edge. Except for a few overestimations around the stagnation point, the predicted ice shape, in terms of ice thickness and freezing locations, is generally comparable to either the LEWICE or FENSAP-ICE profiles.

Validation of Glaze Ice
Accumulation. The glaze icing case corresponds to LEWICE Run 403. As presented in Table 1, the meteorological temperature is prescribed to be 262.04 K. At this temperature, the final ice shape is appreciably irregular, with the formation of two evident ice horns, characterized as typical glaze ice accretion. Figure 8 presents    11 International Journal of Aerospace Engineering a comparison between the current results with LEWICE, FENSAP-ICE, and experimental data. Compared to the IRT experiment, all three numerical codes have noticeable discrepancies on either the pressure or the suction side of the airfoil. It is accepted that a consensus is rarely achieved when comparing two ice accretion profiles [37,40,41]. The current model well predicts the two horns of glaze icing and the ice thickness around the leading edge. Although the simulated run-back behavior is less accurately captured on the suction side, the deviation of the predicted result is within a reasonable and acceptable range. To summarize, the developed icing code is concluded to be reliable for further investigations.

Problem Description and Modeling
Pitot-static systems are speed detection instruments used in the aviation industry. In this study, a typical configuration with an "L" shape is considered. Since the total pressure holes undergo complete failure, the effect of total pressure hole on the blockage of pitot-tube is investigated. Therefore, only the total pressure holes are considered in this research; static holes are neglected. Based on NASA's technical report for pitot-static systems [42], the geometry under consideration can be simplified as shown in Figure 9.
The computational domain is C-type, with a diameter 30 times larger than that of the pitot tube. The orientation of the coordinate axis is presented in Figure 10. The symmetry plane is located at Z = 0 of the XY coordinate plane.
Different ice accretion scenarios are investigated to explore the performance degradation of the pitot tubes induced by different in-flight conditions. The inflow conditions are summarized in Table 4. The pitot tube cases are classified into five groups. The AoA varies from 0°to 10°f or rime icing at T = 241 K (Group 1) and glaze icing at T = 262 K (Group 3). At the fixed AoA = 0°, the inlet temperature ranges from T = 253 K to 271 K, and it increases by 3 K for each case (Group 5). The inlet velocity covers a wide speed range from 100 m/s to 225 m/s (Group 2 for rime icing and Group 4 for glaze icing). For Group 6, the wall 13 International Journal of Aerospace Engineering temperature of the pitot tubes is increased from 264 K to 274 K, with increments of 2 K. Except for the concerned controlled parameters, other inflow conditions within the same group remain unchanged, as presented in Table 5. Based on the criteria of turbulent near-wall treatment, the normal spacing is assigned to be 10 -6 m from the body surface to ensure that y + has an order of magnitude of one. The predicted collection efficiency, using the three mesh densities along the symmetry plane, is compared in Figure 11. The plot exhibits larger distinctions between coarse meshes and the other two types of meshes. Moreover, differences between the medium and fine meshes are minimal. In summary, the medium mesh system is considered suitable for the icing simulations, as shown in Figure 12.

Performance Degradation.
Typically, the accumulated ice on pitot tubes is considered to have a limited influence on speed measurements until it completely blocks the total pressure holes at the leading edges. Therefore, the time at which the total pressure holes are blocked is investigated as a criterion for performance deterioration. For straightforward comparisons, the final blocked pitot tubes due to ice accretion under distinct flight conditions (with AoA of 0°and 10°) are presented in Figures 14 and 15 for the rime and glaze icing conditions, respectively. An increase in the AoA leads to a diminution of ice thickness on the upper side of the leading edge. Additionally, variations in the AoA induce different failure times under either rime icing or glaze icing. Based on the evolution of icing, it can also be deduced that the glaze icing scenario requires a greater amount of time for malfunctions to occur (i.e., the total pressure hole is completely blocked), as compared to the rime icing scenario. For glaze ice accretion, as shown in Figure 15(a), two modest horns can be observed on the upper and lower leading edges of the pitot-static tube. Under the rime icing condition, as shown in Figure 14, a thicker and larger lump of the ice layer emerges around the leading edge, as compared to that under the glaze icing condition in Figure 15, for both AoA.

Blocking Time.
For a quantitative evaluation, the failure times are summarized in Figure 16 and are also discussed in this section. As shown in Figure 16(a), the cases under rime icing have a significantly shorter failure time than those under glaze icing. In addition, the blocking time of the pitot probe has a similar response to the variation in the AoA, under both the rime and glaze icing conditions. With an increase in the AoA, the tube gains additional ice on the upper part of the total pressure hole, while the shape of the ice on the lower wall remains unchanged; this results in a short failure time. Figure 16(b) presents the accelerated increase in the failure time in response to the increase in the inflow temperature. When the inflow temperature approaches the freezing point, the blocking time increases significantly and exceeds 1300 s. It is demonstrated that ambient temperature has a dominant influence on ice accretion processes.
As illustrated in Figure 16(c), the blocking times under the glaze and rime icing conditions exhibit different trends in response to the inflow velocity. For glaze icing, an increase in the inflow velocity reduces the failure time at the low airflow speed range and then modestly enhances it at the higher speed range. An increase in the airflow velocity promotes convective heat transfer (negative value). In addition, the water droplet velocity increases with the airflow velocity, resulting in increments in the mass fluxes, enthalpy (negative value), and kinetic energy (positive value) brought by the collected water. Under low airflow speeds, the enthalpy of droplets and the convective heat flux are enlarged more than the kinetic energy, leading to an acceleration of the failure. Further increase in the inflow velocities causes greater amplifications in the kinetic energy than those in the enthalpy and convective heat flux, resulting in an improved functional performance. For rime icing, the failure time decreases monotonically with an increase in the inlet velocity. As droplets directly freeze into ice, heat exchange is not essential during rime icing processes. Therefore, rime ice thickness is only related to the incremental impinging droplet mass fluxes when the inflow velocity increases. The failure performance is observed to vary linearly with the wall temperature, as indicated in Figure 16(d). Wall temperature is directly associated with the calculation of convective heat transfer during the ice accretion process. As evidenced by the comparison, an increase in wall temperature effectively inhibits the performance deterioration of pitot tubes. Therefore, wall heating is concluded to be an efficient anti-icing approach for in-flight aircrafts and instruments [43].

Regression Analysis.
A regression analysis is conducted to further elucidate the correlations among the abovementioned relevant parameters and the failure time of the pitot tube. As rime ice is generated immediately after the droplets contact the leading surfaces, the accreted ice is closely related to the upwind areas and collected water mass fluxes. Therefore, rime icing processes are affected by variations in the AoA and inflow velocity, irrespective of the inflow and wall temperatures. For the glaze icing scenarios, nonlinear regression analyses deal with all these four parameters.
First, the influential parameters are assumed to be independent, without any interactive effects. Based on the numerical outcomes (Cases 1-36), nonlinear regression analyses involving all the influential parameters are performed. Once the parameters are transformed to dimensionless forms, the blocking time for the pitot tube can be determined as follows: where t 0 = 300 s, T 0 = 262 K, V 0 = 100 m/s, AoA 0 = 10°, and T w0 = 272 K are the selected baseline values for nondimensionalization. The regression analyses indicate good fits to the data, with R 2 = 0:996 and R 2 = 0:999 under the glaze icing and rime icing scenarios, respectively; R 2 is a statistical measure for representing how well the variance in the variables is explained by the regression model [44].
As discussed above, the inlet temperature has a significant influence on ice accretion processes and causes an almost exponential growth in the blocking time. For different inlet velocities, the quadratic function can well describe variation trends, as it is related to the impinging mass and energy fluxes. For cases with different AoAs, the glaze icing scenarios follow a quadratic curve, whereas the rime icing scenarios vary linearly. As shown in Figure 16(d), the linear expression is adequate for identifying the tendencies under different wall temperatures, which is in agreement with actual icing processes. Figure 17 compares the numerical data and the empirical prediction obtained via the proposed equations. It is evident that the results are in good agreement for a majority of the cases, with the exception of Case 29. The acquired equations are competent to predict the blocking time of the pitot tube; in other words, provided the initial conditions within the applicable ranges (i.e., inflow temperature of 243-271 K, inflow velocity of 100-225 m/s, AoA of 0-10°, and wall tem-perature of 264-272 K), the failure time can be well predicted by using the proposed empirical equations.

Conclusions
In this study, in-flight ice accretion on pitot-static systems is numerically investigated by means of a developed numerical approach. A complete set of case studies, involving different ambient and wall temperatures, AoA, and inflow velocities, is considered to reveal the performance deterioration of a pitot tube under rime and glaze ice accretion.
The obtained results indicate that ambient temperature has the most significant impact on the icing-induced failure time, leading to an almost exponential growth. The blocking time is found to be linearly proportional to increments in wall temperature. Therefore, wall heating was concluded to be an efficient anti-icing approach for in-flight aircrafts and instruments. For cases involving different AoAs and inflow velocities, the failure time under the rime icing scenario is considerably shorter than that under the glaze icing scenario. Meanwhile, with an increase in the inflow velocity, the failure time undergoes a parabolic variation under glaze ice accretion; however, it exhibits a monotonic reduction under rime icing conditions. With an increase in the AoA, both the glaze and rime icing scenarios have a tendency to accelerate failure.