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

Authors: Hao Zhang; Rui Zhao; Chih-Yung Wen
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

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.


---

## Body

## 1. 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 in-flight 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–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–21], studies also focused on broadening their engineering applications [22–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 Section2. 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 in-flight 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.


## 2. 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.Figure 1
Interactive loop of the proposed numerical model.


### 2.1. 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(1)∂ρ∂t+∇⋅ρU=0,Momentum equation(2)∂ρU∂t+∇⋅ρUU=‐∇⋅P+23μ∇⋅U+∇⋅μ∇U+∇UT,Energy equation(3)∂ρe∂t+∇⋅ρeU+12∂ρU2∂t+12∇⋅ρUU2+∂P∂t+∇⋅PU−P∇⋅U=−∇⋅23μ∇⋅UU+∇⋅μ∇U+∇UT⋅U+∇⋅λ∇T,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 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
(4)DDtρν~=Cb1ρS~ν~1−ft2−Cw1fw−Cb1κ2ft2ρν~2d2+1σνt∇⋅μ+ν~∇ν~+Cb2ρ∇ν~2+Sν~,where d=dmin+Ks signifies the nearest distance from walls (dmin) associated with the equivalent sand-grain roughness (Ks). The term fw is given as follows: fw=g1+Cw36/g2+Cw361/6, g=r+Cw2r6−r, r≡ν~/S~κ2d2. And ft2=Ct3exp−Ct4χ2, where Ct3=1.2, Ct4=0.5, and χ=ν~/ν+Cr1ks/d, Cr1=0.5. Cb1=0.1355, Cw1=Cb1/κ2+1+Cb2/σvt, κ=0.41, σνt=2/3, Cb2=0.622, Cv1=7.1, Cw3=2, and Cw2=0.3 are constants. The additional definition of S~ is relevant to the vorticity S as S~=S+ν~/κ2d2fv2, where fv2=1−ν~/ν+ν~fv1 and S=2ΩijΩij. The time-averaged rotation tensor is evaluated by Ωij=1/2∂ui/∂xj−∂uj/∂xi.The turbulent viscosity is modeled asμt=ρν~fv1, wherefv1=χ3/χ3+Cv13, Cv1=7.1. The wall boundary condition is ascertained via ∂ν~/∂n=ν~/d.The temporal and spatial discretization of the solution domain is implemented via the finite volume method. The pressure-based PIMPLE algorithm (combination of Pressure Implicit Split Operator (PISO) and Semi Implicit Methods Pressure Linked Equations (SIMPLE) algorithm) is adopted for pressure-velocity coupling [33]. The second order linear upwind scheme is employed to discretize the convective and diffusion terms.


### 2.2. 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 cost-effective 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:(5)∂α∂t+∇⋅αUp=0,(6)∂αUp∂t+∇⋅αUpUp=FU−Up,where α denotes the nondimensional local volume fraction of droplets to air phase and Up is the droplet velocity vector. The source term of drag components induced by the ambient airflow is formulated byF=αfRer/τpU−Up, where Rer=U−Upρpdp/μ characterizes the Reynolds number of droplets, dp is the droplet diameter and refers to the mean volume diameter (MVD) in this research, the momentum response timeτp=ρpdp2/18μ is an inertial parameter to describe variations of droplet velocities, and fRer is a nonlinear empirical formula of drag coefficient, estimated by fRer=1+0.15Rer0.687+0.0175Rer/1+45000Rer−1.16 [34].Taking formality of the Eulerian averaged approach, the nondimensional local collection efficiency is expressed as(7)β=ṀimpṀ∞=−αρpUp⋅nLWC×Up∞,where Up∞ 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.


### 2.3. Ice Accretion Module

#### 2.3.1. Three-Dimensional Ice Accretion Thermodynamic 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:
(8)Ṁimp+Ṁin=Ṁice+Ṁevap+Ṁout,(9)Q̇ice+Q̇imp+Q̇in=Q̇out+Q̇evap+Q̇conv.Figure 2
Mass (a) and energy (b) balance of a control volume during in-flight ice accretion.
(a)(b)The incoming mass flux is a summation of the entering run-back water mass flux from the upstream control volumes (Ṁin) and the impinging droplet mass flux (Ṁimp). The outgoing mass flux includes the leaving run-back water mass flux to the downstream control volumes (Ṁout), freezing mass flux (Ṁice), and evaporating/sublimating mass flux (Ṁevap). 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 ofṀ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
(10)Ṁout,3=Ṁout⋅l3τ⋅n3l3τ⋅n3+l4τ⋅n4,(11)Ṁout,4=Ṁout⋅l4τ⋅n4l3τ⋅n3+l4τ⋅n4,where l3 and l4 denote the lengths of edges 3 and 4, respectively, τ is the shear stress obtained from the airflow solver, and n3 and n4 are the edge normal vectors of edges 3 and 4, respectively.Figure 3
Schematic of flow directions in an individual control volume.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:
(12)Q̇imp=ṀimpCp,waterT∞−T0+Up22,(13)Q̇in=∑nṀin,nCp,waterTin,n−T0,where T∞ is the freestream temperature, T0 is the freezing temperature, and Tin,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 volumes (Q̇out), convective heat transfer with the surrounding airflow (Q̇conv), and heat loss by evaporation or sublimation (Q̇evap). Resorting to the near-wall flow calculation, Q̇conv refers to
(14)Q̇conv=hcTs−T∞=λ∂T∂nwallTs−T∞Tw,ini−T∞,where hc is the convective coefficient, ∂T/∂n is the normal gradient of temperature distribution to the aerodynamic surface, and Tw,ini is the initial temperature of the wall.The above conservative formulas include three unknowns: the mass flux of accreted iceṀice, the flow-out run-back water mass flux Ṁout, and the equilibrium temperature of the aerodynamic surface Ts. The freezing fraction f is therefore introduced to circumvent the inadequacy due to these unknowns and is defined as f≡freezingmass/incomingmass=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 downstream control volumes. The compatibility relation can be given as follows:
(15)f≤0,Ṁice=0,Ṁout≠0,Ts>0;Q̇ice=0,Q̇out=ṀoutCp,waterTs−T0,Q̇evap=ṀevapLevap+Cp,waterTs−T0.If the equilibrium temperature is the freezing point, the ice and water are coexistent in the current control volume, and the compatibility relation becomes(16)0<f<1,Ṁice≠0,Ṁout≠0,Ts=T0;Q̇ice=ṀiceLf,Q̇out≈0,Q̇evap=ṀevapLevap+Cp,waterTs−T0.If the equilibrium temperature is below the freezing point, all captured droplets freeze with no water remaining in the current control volume, and the relation is changed as given below:(17)f≥1,Ṁice≠0,Ṁout=0,Ts<0;Q̇ice=ṀiceLf+Cp,iceT0−Ts,Q̇out=0,Q̇evap=ṀevapLsub+Cp,iceTs−T0.Here,Cp,ice is the specific heat of ice, Lf is the fusion latent heat of water, and Levap and Lsub are the latent heat with respect to the freezing point for vaporization and sublimation, respectively.


#### 2.3.2. 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.Figure 4
Flow chart depicting the solving strategy of the developed ice accretion solver.After the thermodynamic iteration, the ice thickness is finalized by(18)Hice=Ṁice,n⋅Δtρice,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.


## 2.1. 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(1)∂ρ∂t+∇⋅ρU=0,Momentum equation(2)∂ρU∂t+∇⋅ρUU=‐∇⋅P+23μ∇⋅U+∇⋅μ∇U+∇UT,Energy equation(3)∂ρe∂t+∇⋅ρeU+12∂ρU2∂t+12∇⋅ρUU2+∂P∂t+∇⋅PU−P∇⋅U=−∇⋅23μ∇⋅UU+∇⋅μ∇U+∇UT⋅U+∇⋅λ∇T,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 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
(4)DDtρν~=Cb1ρS~ν~1−ft2−Cw1fw−Cb1κ2ft2ρν~2d2+1σνt∇⋅μ+ν~∇ν~+Cb2ρ∇ν~2+Sν~,where d=dmin+Ks signifies the nearest distance from walls (dmin) associated with the equivalent sand-grain roughness (Ks). The term fw is given as follows: fw=g1+Cw36/g2+Cw361/6, g=r+Cw2r6−r, r≡ν~/S~κ2d2. And ft2=Ct3exp−Ct4χ2, where Ct3=1.2, Ct4=0.5, and χ=ν~/ν+Cr1ks/d, Cr1=0.5. Cb1=0.1355, Cw1=Cb1/κ2+1+Cb2/σvt, κ=0.41, σνt=2/3, Cb2=0.622, Cv1=7.1, Cw3=2, and Cw2=0.3 are constants. The additional definition of S~ is relevant to the vorticity S as S~=S+ν~/κ2d2fv2, where fv2=1−ν~/ν+ν~fv1 and S=2ΩijΩij. The time-averaged rotation tensor is evaluated by Ωij=1/2∂ui/∂xj−∂uj/∂xi.The turbulent viscosity is modeled asμt=ρν~fv1, wherefv1=χ3/χ3+Cv13, Cv1=7.1. The wall boundary condition is ascertained via ∂ν~/∂n=ν~/d.The temporal and spatial discretization of the solution domain is implemented via the finite volume method. The pressure-based PIMPLE algorithm (combination of Pressure Implicit Split Operator (PISO) and Semi Implicit Methods Pressure Linked Equations (SIMPLE) algorithm) is adopted for pressure-velocity coupling [33]. The second order linear upwind scheme is employed to discretize the convective and diffusion terms.


## 2.2. 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 cost-effective 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:(5)∂α∂t+∇⋅αUp=0,(6)∂αUp∂t+∇⋅αUpUp=FU−Up,where α denotes the nondimensional local volume fraction of droplets to air phase and Up is the droplet velocity vector. The source term of drag components induced by the ambient airflow is formulated byF=αfRer/τpU−Up, where Rer=U−Upρpdp/μ characterizes the Reynolds number of droplets, dp is the droplet diameter and refers to the mean volume diameter (MVD) in this research, the momentum response timeτp=ρpdp2/18μ is an inertial parameter to describe variations of droplet velocities, and fRer is a nonlinear empirical formula of drag coefficient, estimated by fRer=1+0.15Rer0.687+0.0175Rer/1+45000Rer−1.16 [34].Taking formality of the Eulerian averaged approach, the nondimensional local collection efficiency is expressed as(7)β=ṀimpṀ∞=−αρpUp⋅nLWC×Up∞,where Up∞ 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.


## 2.3. Ice Accretion Module

### 2.3.1. Three-Dimensional Ice Accretion Thermodynamic 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:
(8)Ṁimp+Ṁin=Ṁice+Ṁevap+Ṁout,(9)Q̇ice+Q̇imp+Q̇in=Q̇out+Q̇evap+Q̇conv.Figure 2
Mass (a) and energy (b) balance of a control volume during in-flight ice accretion.
(a)(b)The incoming mass flux is a summation of the entering run-back water mass flux from the upstream control volumes (Ṁin) and the impinging droplet mass flux (Ṁimp). The outgoing mass flux includes the leaving run-back water mass flux to the downstream control volumes (Ṁout), freezing mass flux (Ṁice), and evaporating/sublimating mass flux (Ṁevap). 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 ofṀ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
(10)Ṁout,3=Ṁout⋅l3τ⋅n3l3τ⋅n3+l4τ⋅n4,(11)Ṁout,4=Ṁout⋅l4τ⋅n4l3τ⋅n3+l4τ⋅n4,where l3 and l4 denote the lengths of edges 3 and 4, respectively, τ is the shear stress obtained from the airflow solver, and n3 and n4 are the edge normal vectors of edges 3 and 4, respectively.Figure 3
Schematic of flow directions in an individual control volume.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:
(12)Q̇imp=ṀimpCp,waterT∞−T0+Up22,(13)Q̇in=∑nṀin,nCp,waterTin,n−T0,where T∞ is the freestream temperature, T0 is the freezing temperature, and Tin,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 volumes (Q̇out), convective heat transfer with the surrounding airflow (Q̇conv), and heat loss by evaporation or sublimation (Q̇evap). Resorting to the near-wall flow calculation, Q̇conv refers to
(14)Q̇conv=hcTs−T∞=λ∂T∂nwallTs−T∞Tw,ini−T∞,where hc is the convective coefficient, ∂T/∂n is the normal gradient of temperature distribution to the aerodynamic surface, and Tw,ini is the initial temperature of the wall.The above conservative formulas include three unknowns: the mass flux of accreted iceṀice, the flow-out run-back water mass flux Ṁout, and the equilibrium temperature of the aerodynamic surface Ts. The freezing fraction f is therefore introduced to circumvent the inadequacy due to these unknowns and is defined as f≡freezingmass/incomingmass=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 downstream control volumes. The compatibility relation can be given as follows:
(15)f≤0,Ṁice=0,Ṁout≠0,Ts>0;Q̇ice=0,Q̇out=ṀoutCp,waterTs−T0,Q̇evap=ṀevapLevap+Cp,waterTs−T0.If the equilibrium temperature is the freezing point, the ice and water are coexistent in the current control volume, and the compatibility relation becomes(16)0<f<1,Ṁice≠0,Ṁout≠0,Ts=T0;Q̇ice=ṀiceLf,Q̇out≈0,Q̇evap=ṀevapLevap+Cp,waterTs−T0.If the equilibrium temperature is below the freezing point, all captured droplets freeze with no water remaining in the current control volume, and the relation is changed as given below:(17)f≥1,Ṁice≠0,Ṁout=0,Ts<0;Q̇ice=ṀiceLf+Cp,iceT0−Ts,Q̇out=0,Q̇evap=ṀevapLsub+Cp,iceTs−T0.Here,Cp,ice is the specific heat of ice, Lf is the fusion latent heat of water, and Levap and Lsub are the latent heat with respect to the freezing point for vaporization and sublimation, respectively.


### 2.3.2. 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.Figure 4
Flow chart depicting the solving strategy of the developed ice accretion solver.After the thermodynamic iteration, the ice thickness is finalized by(18)Hice=Ṁice,n⋅Δtρice,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.


## 2.3.1. Three-Dimensional Ice Accretion Thermodynamic 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:
(8)Ṁimp+Ṁin=Ṁice+Ṁevap+Ṁout,(9)Q̇ice+Q̇imp+Q̇in=Q̇out+Q̇evap+Q̇conv.Figure 2
Mass (a) and energy (b) balance of a control volume during in-flight ice accretion.
(a)(b)The incoming mass flux is a summation of the entering run-back water mass flux from the upstream control volumes (Ṁin) and the impinging droplet mass flux (Ṁimp). The outgoing mass flux includes the leaving run-back water mass flux to the downstream control volumes (Ṁout), freezing mass flux (Ṁice), and evaporating/sublimating mass flux (Ṁevap). 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 ofṀ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
(10)Ṁout,3=Ṁout⋅l3τ⋅n3l3τ⋅n3+l4τ⋅n4,(11)Ṁout,4=Ṁout⋅l4τ⋅n4l3τ⋅n3+l4τ⋅n4,where l3 and l4 denote the lengths of edges 3 and 4, respectively, τ is the shear stress obtained from the airflow solver, and n3 and n4 are the edge normal vectors of edges 3 and 4, respectively.Figure 3
Schematic of flow directions in an individual control volume.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:
(12)Q̇imp=ṀimpCp,waterT∞−T0+Up22,(13)Q̇in=∑nṀin,nCp,waterTin,n−T0,where T∞ is the freestream temperature, T0 is the freezing temperature, and Tin,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 volumes (Q̇out), convective heat transfer with the surrounding airflow (Q̇conv), and heat loss by evaporation or sublimation (Q̇evap). Resorting to the near-wall flow calculation, Q̇conv refers to
(14)Q̇conv=hcTs−T∞=λ∂T∂nwallTs−T∞Tw,ini−T∞,where hc is the convective coefficient, ∂T/∂n is the normal gradient of temperature distribution to the aerodynamic surface, and Tw,ini is the initial temperature of the wall.The above conservative formulas include three unknowns: the mass flux of accreted iceṀice, the flow-out run-back water mass flux Ṁout, and the equilibrium temperature of the aerodynamic surface Ts. The freezing fraction f is therefore introduced to circumvent the inadequacy due to these unknowns and is defined as f≡freezingmass/incomingmass=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 downstream control volumes. The compatibility relation can be given as follows:
(15)f≤0,Ṁice=0,Ṁout≠0,Ts>0;Q̇ice=0,Q̇out=ṀoutCp,waterTs−T0,Q̇evap=ṀevapLevap+Cp,waterTs−T0.If the equilibrium temperature is the freezing point, the ice and water are coexistent in the current control volume, and the compatibility relation becomes(16)0<f<1,Ṁice≠0,Ṁout≠0,Ts=T0;Q̇ice=ṀiceLf,Q̇out≈0,Q̇evap=ṀevapLevap+Cp,waterTs−T0.If the equilibrium temperature is below the freezing point, all captured droplets freeze with no water remaining in the current control volume, and the relation is changed as given below:(17)f≥1,Ṁice≠0,Ṁout=0,Ts<0;Q̇ice=ṀiceLf+Cp,iceT0−Ts,Q̇out=0,Q̇evap=ṀevapLsub+Cp,iceTs−T0.Here,Cp,ice is the specific heat of ice, Lf is the fusion latent heat of water, and Levap and Lsub are the latent heat with respect to the freezing point for vaporization and sublimation, respectively.


## 2.3.2. 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.Figure 4
Flow chart depicting the solving strategy of the developed ice accretion solver.After the thermodynamic iteration, the ice thickness is finalized by(18)Hice=Ṁice,n⋅Δtρice,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.


## 3. 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.


### 3.1. 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].Table 1
Setup and inlet conditions of the airflow validation case over a 2D NACA0012 airfoil.
GeometryNACA0012—1 mVelocity51 m/sAoA10°Pressure101,325 PaDensity1.18 kg/m3Temperature300 KMach number0.15Reynolds number3,500,000Computational domainO-type, diameter: 30 mMesh number52,357Figure 5
Comparison of simulation and experimental results in airflow fields.


### 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/m3. 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 2
Setup and inlet conditions of the droplet impingement validation case over a 2D NACA0012 airfoil.
GeometryNACA0012—0.5334 mVelocity102.80 m/sAoA4°Pressure100,920 PaDensity1.46 kg/m3Temperature262.04 KMVD20μmLWC0.55 g/m3Reynolds number5,100,000Computational domainO-type, diameter -16 mMesh number49,176Figure 6
Comparison of collection efficiency among the numerical data obtained by the current code, LEWICE, and FENSAP-ICE.


### 3.3. Validation of Ice Accretion Solver

The developed ice accretion solver is validated by performing two representative ice-accumulating benchmarks, i.e., LEWICE Run 403 and Run 404 [8, 39], corresponding to the rime ice and glaze ice phenomena, respectively.


#### 3.3.1. Validation of Rime Ice Accretion

The ambient temperature of the LEWICE Run 404 case is set as 256.49 K. The numerical arrangement follows that in the Icing Research Tunnel (IRT) experiment, as detailed in Table3.Table 3
Settings for the rime ice accretion case.
GeometryNACA0012—0.5334 mVelocity102.8 m/sIDLEWICE Run 404AoA4°Pressure100,920 PaDensity1.457 kg/m3Temperature256.49 KMVD20μmLWC0.55 g/m3Reynolds number5,100,000Computational domainC-type, diameter: 16 mMesh number49,176Ice accretion time7 minFigure7 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.Figure 7
Comparison of rime ice accretion for the current numerical, LEWICE, FENSAP-ICE, and IRT experimental data.


#### 3.3.2. Validation of Glaze Ice Accumulation

The glaze icing case corresponds to LEWICE Run 403. As presented in Table1, 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 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.Figure 8
Comparison of glaze ice accretion between the current numerical result with LEWICE, FENSAP-ICE, and IRT experimental data.


## 3.1. 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].Table 1
Setup and inlet conditions of the airflow validation case over a 2D NACA0012 airfoil.
GeometryNACA0012—1 mVelocity51 m/sAoA10°Pressure101,325 PaDensity1.18 kg/m3Temperature300 KMach number0.15Reynolds number3,500,000Computational domainO-type, diameter: 30 mMesh number52,357Figure 5
Comparison of simulation and experimental results in airflow fields.


## 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/m3. 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 2
Setup and inlet conditions of the droplet impingement validation case over a 2D NACA0012 airfoil.
GeometryNACA0012—0.5334 mVelocity102.80 m/sAoA4°Pressure100,920 PaDensity1.46 kg/m3Temperature262.04 KMVD20μmLWC0.55 g/m3Reynolds number5,100,000Computational domainO-type, diameter -16 mMesh number49,176Figure 6
Comparison of collection efficiency among the numerical data obtained by the current code, LEWICE, and FENSAP-ICE.


## 3.3. Validation of Ice Accretion Solver

The developed ice accretion solver is validated by performing two representative ice-accumulating benchmarks, i.e., LEWICE Run 403 and Run 404 [8, 39], corresponding to the rime ice and glaze ice phenomena, respectively.


### 3.3.1. Validation of Rime Ice Accretion

The ambient temperature of the LEWICE Run 404 case is set as 256.49 K. The numerical arrangement follows that in the Icing Research Tunnel (IRT) experiment, as detailed in Table3.Table 3
Settings for the rime ice accretion case.
GeometryNACA0012—0.5334 mVelocity102.8 m/sIDLEWICE Run 404AoA4°Pressure100,920 PaDensity1.457 kg/m3Temperature256.49 KMVD20μmLWC0.55 g/m3Reynolds number5,100,000Computational domainC-type, diameter: 16 mMesh number49,176Ice accretion time7 minFigure7 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.Figure 7
Comparison of rime ice accretion for the current numerical, LEWICE, FENSAP-ICE, and IRT experimental data.


### 3.3.2. Validation of Glaze Ice Accumulation

The glaze icing case corresponds to LEWICE Run 403. As presented in Table1, 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 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.Figure 8
Comparison of glaze ice accretion between the current numerical result with LEWICE, FENSAP-ICE, and IRT experimental data.


## 3.3.1. Validation of Rime Ice Accretion

The ambient temperature of the LEWICE Run 404 case is set as 256.49 K. The numerical arrangement follows that in the Icing Research Tunnel (IRT) experiment, as detailed in Table3.Table 3
Settings for the rime ice accretion case.
GeometryNACA0012—0.5334 mVelocity102.8 m/sIDLEWICE Run 404AoA4°Pressure100,920 PaDensity1.457 kg/m3Temperature256.49 KMVD20μmLWC0.55 g/m3Reynolds number5,100,000Computational domainC-type, diameter: 16 mMesh number49,176Ice accretion time7 minFigure7 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.Figure 7
Comparison of rime ice accretion for the current numerical, LEWICE, FENSAP-ICE, and IRT experimental data.


## 3.3.2. Validation of Glaze Ice Accumulation

The glaze icing case corresponds to LEWICE Run 403. As presented in Table1, 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 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.Figure 8
Comparison of glaze ice accretion between the current numerical result with LEWICE, FENSAP-ICE, and IRT experimental data.


## 4. 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.Figure 9
Simplified pitot-static tube configurations considered in this study.The computational domain isC-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.Figure 10
Computation domain and coordinate orientation.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 Table4. The pitot tube cases are classified into five groups. The AoA varies from 0° to 10° for 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 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.Table 4
Ice accretion cases for the pitot-static tube.
Group 1AoA (rime icingT=241 K)Case 1AoA=0°Case 2AoA=2°Case 3AoA=4°Case 4AoA=6°Case 5AoA=8°Case 6AoA=10°Group 2Inflow velocity (rime icingT=241 K)Case 7U=100 m/sCase 8U=125 m/sCase 9U=150 m/sCase 10U=175 m/sCase 11U=200 m/sCase 12U=225 m/sGroup 3AoA (glaze icingT=262 K)Case 13AoA=0°Case 14AoA=2°Case 15AoA=4°Case 16AoA=6°Case 17AoA=8°Case 18AoA=10°Group 4Inflow velocity (glaze icingT=262 K)Case 19U=100 m/sCase 20U=125 m/sCase 21U=150 m/sCase 22U=175 m/sCase 23U=200 m/sCase 24U=225 m/sGroup 5Inflow temperatureCase 25T=253 KCase 26T=256 KCase 27T=259 KCase 28T=265 KCase 29T=268 KCase 30T=271 KGroup 6Wall temperatureCase 31TW=264 KCase 32TW=266 KCase 33TW=268 KCase 34TW=270 KCase 35TW=272 KCase 36TW=274 KTable 5
Simulation settings of the pitot-static tube cases.
GeometryPitot-static tubeDensity1.457 kg/m3Velocity100 m/sLWC0.55 g/m3AoA0°Temperature262 KComputational domainC-type, diameter: 0.762 mWall temperature272 K


## 5. Results and Discussion

### 5.1. Mesh Resolution

The mesh independence studies are conducted in this section. Three levels of mesh refinements are used for the pitot-static tube—coarse, medium, and fine meshes; the corresponding nodes in the computational domain are 1,997,536; 3,699,655; and 5,152,743, respectively. 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.Figure 11
Collection efficiency distribution for the pitot tube with three different mesh resolutions. Arc length is defined as the ratio of curve distance from the original point and the total length of the tubes.Figure 12
Structured mesh distribution over the wall of the pitot-static tube.


### 5.2. Distribution of Droplet Collection

As a preliminary step, the properties of droplets impinging on the pitot-static tubes are discussed. Figure13 presents a comparison of the collection efficiency contours of the pitot-static tube under AoA of 0° (Case 1) and 10° (Case 6). It should be noted that the tube experiences the heaviest droplet impingement at the leading edge. Besides, the other forward-facing surfaces possess high water contents. Compared to the AoA of 0°, the peak values of the collection efficiency shift with the stagnation points under the AoA of 10°. Furthermore, cases of different inflow velocities (Cases 7–12 and 19–24), inflow temperatures (Cases 25–30), and wall temperatures (Cases 31–36) exhibit analogous impinging behaviors, as compared to the cases in Group 1 for water collection.Figure 13
Collection efficiency contour of the pitot-static tube under Case 1 (a) and Case 6 (b).
(a)(b)


### 5.3. 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 Figures14 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.Figure 14
Blocked ice accretion on the pitot-static tube under the rime ice condition at 241 K and different AoAs: (a)AoA=0° at a blocking time of 180 s (Case 1); (b) AoA=10° at a blocking time of 93 s (Case 6).
(a)(b)Figure 15
Blocked ice accretion on the pitot-static tube under the glaze ice condition at 262 K and different AoAs: (a)AoA=0° at a blocking time of 310 s (Case 13); (b) AoA=10° at a blocking time of 279 s (Case 18).
(a)(b)


### 5.4. Blocking Time

For a quantitative evaluation, the failure times are summarized in Figure16 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
Failure time due to ice accretion under (a) different AoAs, (b) different inflow temperatures, (c) different inflow velocities, and (d) different pitot tube wall temperatures.
(a)(b)(c)(d)Figure16(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 Figure16(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 Figure16(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].


### 5.5. 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:(19)tt0Type1=exp73.961TT0−75.179+0.115VV02−0.405VV0−0.223AoAAoA02+0.046AoAAoA0+1.644TWTW0−0.642forglazeicing0.072VV02−0.47VV0−0.048AoAAoA0+1forrimeicing,where t0=300 s, T0=262 K, V0=100 m/s, AoA0=10°, and Tw0=272 K are the selected baseline values for nondimensionalization. The regression analyses indicate good fits to the data, with R2=0.996 and R2=0.999 under the glaze icing and rime icing scenarios, respectively; R2 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 Figure16(d), the linear expression is adequate for identifying the tendencies under different wall temperatures, which is in agreement with actual icing processes.Figure17 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 temperature of 264–272 K), the failure time can be well predicted by using the proposed empirical equations.Figure 17
Comparison of numerical data and results calculated via the proposed equations for a pitot tube.


## 5.1. Mesh Resolution

The mesh independence studies are conducted in this section. Three levels of mesh refinements are used for the pitot-static tube—coarse, medium, and fine meshes; the corresponding nodes in the computational domain are 1,997,536; 3,699,655; and 5,152,743, respectively. 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.Figure 11
Collection efficiency distribution for the pitot tube with three different mesh resolutions. Arc length is defined as the ratio of curve distance from the original point and the total length of the tubes.Figure 12
Structured mesh distribution over the wall of the pitot-static tube.


## 5.2. Distribution of Droplet Collection

As a preliminary step, the properties of droplets impinging on the pitot-static tubes are discussed. Figure13 presents a comparison of the collection efficiency contours of the pitot-static tube under AoA of 0° (Case 1) and 10° (Case 6). It should be noted that the tube experiences the heaviest droplet impingement at the leading edge. Besides, the other forward-facing surfaces possess high water contents. Compared to the AoA of 0°, the peak values of the collection efficiency shift with the stagnation points under the AoA of 10°. Furthermore, cases of different inflow velocities (Cases 7–12 and 19–24), inflow temperatures (Cases 25–30), and wall temperatures (Cases 31–36) exhibit analogous impinging behaviors, as compared to the cases in Group 1 for water collection.Figure 13
Collection efficiency contour of the pitot-static tube under Case 1 (a) and Case 6 (b).
(a)(b)


## 5.3. 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 Figures14 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.Figure 14
Blocked ice accretion on the pitot-static tube under the rime ice condition at 241 K and different AoAs: (a)AoA=0° at a blocking time of 180 s (Case 1); (b) AoA=10° at a blocking time of 93 s (Case 6).
(a)(b)Figure 15
Blocked ice accretion on the pitot-static tube under the glaze ice condition at 262 K and different AoAs: (a)AoA=0° at a blocking time of 310 s (Case 13); (b) AoA=10° at a blocking time of 279 s (Case 18).
(a)(b)


## 5.4. Blocking Time

For a quantitative evaluation, the failure times are summarized in Figure16 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
Failure time due to ice accretion under (a) different AoAs, (b) different inflow temperatures, (c) different inflow velocities, and (d) different pitot tube wall temperatures.
(a)(b)(c)(d)Figure16(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 Figure16(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 Figure16(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].


## 5.5. 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:(19)tt0Type1=exp73.961TT0−75.179+0.115VV02−0.405VV0−0.223AoAAoA02+0.046AoAAoA0+1.644TWTW0−0.642forglazeicing0.072VV02−0.47VV0−0.048AoAAoA0+1forrimeicing,where t0=300 s, T0=262 K, V0=100 m/s, AoA0=10°, and Tw0=272 K are the selected baseline values for nondimensionalization. The regression analyses indicate good fits to the data, with R2=0.996 and R2=0.999 under the glaze icing and rime icing scenarios, respectively; R2 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 Figure16(d), the linear expression is adequate for identifying the tendencies under different wall temperatures, which is in agreement with actual icing processes.Figure17 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 temperature of 264–272 K), the failure time can be well predicted by using the proposed empirical equations.Figure 17
Comparison of numerical data and results calculated via the proposed equations for a pitot tube.


## 6. 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. Compared to the broadly linear reduction under rime icing, the failure time under glaze icing is more sensitive to changes in the AoA and exhibits larger differences.A nonlinear regression analysis is conducted for quick estimations of failure times for the pitot tube. The correlation equations employ two relevant parameters for rime ice accretion and four relevant parameters for glaze ice accretion, yielding convenient predictions of blocking time within the applicable range.



---

*Source: 5599116-2021-07-19.xml*