Abstract

To avoid the use of a computationally intensive full unsteady simulation while providing an accurate solution, a quasi-steady simulation has been performed to study glaze and mixed ice accretion in the supercooled large droplet (SLD) regime in turbulent flow of atmospheric air. We have attempted to find the minimum time-step necessary to adequately simulate the icing process on an airfoil surface. Based on node displacement, a mesh morphing scheme has been adopted in the computations to account for the moving boundaries that are caused by the continuous ice buildup. We have modeled ice accretion on an airfoil surface for a time period of 232 s using several time steps. At each time-step we solved the steady-state conservation equations for the air and droplet phases and then used the results as initial conditions for the next time position. The magnitude of the time-step ranged from a least accurate two-shot simulation (where the time-step is 116 s) to a most accurate 2320-shot simulation (where the time-step is 0.1 s). In-between these two extreme time steps, we have performed a three-shot simulation (where the time-step is 77.33 s), a four-shot simulation (where the time-step is 58 s), a six-shot simulation (where the time-step is 38.67 s), a 46-shot simulation (where the time-step is 5 s), and a 232-shot simulation (where the time-step is 1 s). We have done so to find out the degree of accuracy (or inaccuracy) of the multishot simulation approach and to find out the appropriate time-step needed for a successful and valid quasi-steady simulation. A valid quasi-steady simulation needs to use a time-step that is small enough to reproduce the full time-dependent solution within a very small error band. We have found that both the 1 s and the 0.1 s time steps produce virtually identical results. This is the primary litmus test that proves the validity of the quasi-steady-state assumption. The results adopted in this paper are thus all based on the more conservative 0.1 s time-step. In the process of performing the simulations, remeshing was required in order to maintain the grid density in zones of high curvature to be able to capture the full physics in those zones. After successfully modeling glaze ice accretion over the airfoil surface using the 0.1 s quasi-steady simulation approach, the effects of supercooled large droplets (SLDs) impacting the surface have been examined and presented in terms of the variation of the local collection efficiency, the water film thickness, and the heat transfer rate. Examination of the variation of the angle that the ice horn makes with the airfoil chord line demonstrated a 20% improvement in angle prediction when the time-step is reduced from 116 s to 0.1 s. The analysis also reveals a 12% or 8.5% increase in the maximum collection efficiency, βmax, depending on whether the value of the liquid water content (LWC) has been doubled from 0.5g/m3 to 1g/m3, or the value of the freestream velocity has been doubled from 75m/s to 150m/s, respectively. Because of the need to monitor the local collection efficiency and convective heat fluxes at each shot (regardless of the number of shots employed), the approach adopted here was found to be effective in successively and successfully reproducing the curvature of the glaze ice horn.

1 Introduction

Any cloud containing liquid water can create a significant icing environment if the temperature is 0 °C or less. Cloud water may still not freeze at that temperature and instead exists as a supercooled liquid. When a flying surface such as an airplane wing or a wind turbine blade is exposed to these airborne droplets, atmospheric icing can occur, reducing the overall performance of the wing or the wind turbine blade. Ice is likely to form when the latent heat of supercooled droplets is removed. In terms of meteorological and wall conditions, ice can be classified into two main types: rime ice and glaze ice. The local icing rate is primarily governed by the ability to remove latent heat and hence the local convective heat transfer. When the ambient temperature is sufficiently low, supercooled droplets freeze quickly upon impact, resulting in rime ice. In glaze ice, which is identified with a strong clear ice structure, the heat transfer is insufficient to freeze all the impinging supercooled droplets instantaneously. Thus, a high proportion of water droplets can still remain liquid and emerge as a film of water running back down the wing surface. The ice that develops on an aircraft wing is normally a mix of rime and glaze ice and is referred to as mixed ice. Ice accretion, particularly the glazing and mixed types, are complex processes that necessitate a thorough understanding of phase transition as well as heat and mass transfer. The heat transfer/freezing processes are governed by several factors in addition to the temperature difference between the exposed surface and the impinging water droplets. These factors include the speed of the airflow, the liquid water content (LWC), the droplet diameter, and the surface roughness, among other factors.

The Messinger thermodynamic model, one of the earliest models in the field of ice accretion simulation, proposes that the icing rate can only be estimated by a mass and energy balance [1] for every single control volume. That model has been utilized in well-known ice accretion codes such as the NASA Lewis (Glenn) Ice Accretion Code (LEWICE) [2] and the Office National d'Etudes et de Recherches Aerospatiales (ONERA) code [3]. Since the thermal conductivity of ice is low, the original Messinger model overlooks the transitory nature of the icing process as well as the conduction heat transfer inside the ice layer. Both of those effects are important when the ice layer is thin as in the case of glaze ice. By eliminating the adiabatic condition of the surface onto which the ice forms, Myers [4] later addressed these concerns and incorporated a temperature gradient within the ice layer. Sherif et al. [5] presented a model that predicts ice accretion and local heat transfer rates when the equilibrium wall temperature is near freezing by accounting for the effects of freezing, sublimation, and evaporation. Naterer [6] developed an analytical solution to investigate the temperature gradient in the unfrozen liquid layer caused by the partial solidification of the incoming droplets. Bourgault et al. [7] introduced a novel thermodynamic model to tackle the shortcomings of the control volume approach by predicting ice accretion and droplet runback on the surface using a set of partial differential equations (PDEs). Zhu et al. [8] demonstrated that both the pressure gradient and the shear stress in the air layer are critical factors in the mass transition of the water film.

Additionally, meteorological examinations into icing-induced aviation accidents reveal the existence of supercooled large droplets (SLDs) that are larger than the median volumetric diameter of 40μm in clouds. Two distinct atmospheric conditions invoke the formation of SLDs: temperature inversion and a collision-coalescence process. These typically result in droplet diameters larger than 40μm. SLDs are thought to be the primary cause of a number of aviation accidents such as the American Eagle ATR-72 accident [9]. Ashenden et al. [10] reported that SLDs such as freezing drizzle may result in significant performance degradation of aircraft. Because of their size, SLDs cause abnormal ice accretion, such as faster ice accretion than normal droplets and striking at different locations. Accretions due to SLDs are harder to predict than accretions from normal size droplets.

Accuracy in predicting the local collection efficiency is a prerequisite for being able to perform accurate icing studies and accurately determining ice deposition shapes and rates on aircraft surfaces. Being able to do so typically requires taking into account aerodynamic droplet deformation and break up, as well as the bouncing and splashing of droplets on the surface. One of the earliest research about the effect of aerodynamic droplet breakup and the impingement of SLDs on an aircraft surface was reported by Honsek et al. [11]. Kim et al. [12] further improved the previous model to be more suitable to the Eulerian–Eulerian approach. Wang et al. [13] proposed a 2D splashing model to address the lack of investigations into the interaction between SLD dynamics and ice accretion. Raj et al. [14] studied the effect of SLD dynamics, specifically splashing and deformation, on a multi-element airfoil.

Therefore, precise modeling of ice shape and how it varies with environmental conditions is critical for being able to predict its occurrence and mitigating its negative impact. Numerical approaches are adequately accurate and cost-effective to provide intricate information on many features of fluid flow and the associated heat transfer mechanisms that may be difficult and/or costly to obtain through experiments. Stebbins et al. [15] carried out computational fluid dynamics (CFD) investigations that focused on studying the aerodynamic influence of in-flight icing and discussed the impact of wing shape, flow conditions, numerical schemes, and other factors. Many researchers have used CFD methods to replicate the ice buildup process on surfaces [16,17]. Using open-source software, Gori et al. [18] proposed an ice accretion model aimed at improving inconsistencies of the Myers model [4] by taking into account the local value of the air temperature outside the boundary layer rather than the constant freestream temperature value to solve the unsteady Stefan problem. Cao et al. [19] numerically investigated the natural melting behavior to study de-icing methods by instantly exposing the iced-up wing to a relatively higher surrounding temperature. Their results showed that the melting rate increased with the Mach number. Shad and Sherif [20] introduced a numerical methodology for simulating rime ice accretion on an airfoil surface in which all required modules such as airflow and droplet fields, icing growth, and the iced surface displacement were fully unsteadily coupled. Li and Paoli [21] used a compressible openfoam solver for simulating ice accretion over an aircraft wing in which their employed mesh-morphing model was capable of altering the form of the wing at each time-step.

However, a full unsteady 3D simulation in the turbulent flow regime is highly demanding in terms of computational time, even with today's super computers. Quasi-steady-state simulations where the steady-state governing equations are solved for sufficiently small time intervals (steps) and then advanced in time to the next time position while utilizing the steady-state solution at the previous time position as an initial condition to the current time position is a more cost-effective and highly accurate strategy. This approach is sometimes referred to as a multishot approach where the airflow and droplet fields are updated at the beginning of each time position and then kept constant until the next one. A “sufficiently small” time interval (step) is one that is small enough so as to produce a solution to the governing equations within a preselected very small enough error band relative to the full time-dependent solution. In that approach, the computational domain is updated to reflect the change in the ice-covered geometry. Studies that employ a version of that approach are several. For example, Ozcer et al. [22] used a multishot approach using the finite element Navier–Stokes analysis package (FENSAP-ICE) to perform a long-exposure 3D icing simulation where the results obtained for both the rime and glaze ice types were well supported by the experimental data, thus demonstrating the robustness of the method. Aliaga et al. [23] achieved a significant reduction in the overall computational time while maintaining numerical accuracy with respect to the maximum allowable time-step for a 2D rime ice case.

The contribution of the work described in this paper is thus multifold. Because of the complexity of the problem in terms of its highly time-dependent nature and because of the fact that the glaze and mixed ice accretion modes are the dominant modes, a full time-dependent solution of the governing equations will be costly as far as computational time is concerned. Added to the complexity is the fact that we're-also investigating the supercooled large droplet regime. Using a quasi-steady-state approach to solve the problem will help achieve virtually the same accuracy without having to endure the computational time penalty. As indicated earlier, a quasi-steady-state approach is one where the steady-state governing equations are solved for sufficiently small time intervals, and then the solution is advanced in time using the solution at the previous time position as the initial condition for the new time position. A valid quasi-steady simulation, therefore, needs to employ a time-step size that is small enough to be able to reproduce the full time-dependent solution within a very small error band. The terms “small enough” or “sufficiently small” necessarily mean that the quasi-steady-state solution should virtually appear no different than the full time-dependent solution. So, one of the contributions of this paper is to investigate several quasi-steady-state solutions (or several multishot solutions) to find out how small the time-step have to be to obtain a meaningfully accurate solution. The size of the time-step is thus an important variable that we needed to investigate. Another contribution is the combined use of the mesh morphing approach with a manual remeshing process that is capable of capturing the ice accretion physics when large deformations to the icing profile would have occurred (after the passing of enough time). This has been successfully applied in this paper and will have the advantage of efficiency without sacrificing the ability to capture the physics when large ice profile deformations have occurred (after the lapse of sufficient time). Said differently, the mesh morphing approach is fast and efficient to reconstruct the new mesh but it is incapable of producing an accurate mesh when large deformations have occurred. Remeshing is a more accurate way but its computational cost is high when it needs to be done for every time-step. Last but not least, the paper addresses important and complex ice accretion modes, which are glaze and mixed ice accretion. These complex accretion modes have been made more complex due to the fact that they are happening when supercooled large droplets are impacting the surface. Modeling supercooled large droplets add significantly more complexities to the physics of the problem.

2 Computational Approach

The NACA 0012 airfoil with a chord length c=53.34cm and a maximum thickness is 12% of the chord length, a typical CFD validation example, is used for simulations in this paper. The problem here is the transport of a liquid phase in a carrier air flow, where an Eulerian solution requires solving a coupled partial differential equation system for each phase. Several simplifying assumptions can generally be used to reduce the complexity of the computational solution for a typical set of conditions that apply to aircraft icing. These are discussed in more detail in this section.

2.1 Air Flow Modeling.

The flow field is considered a dilute air-droplet mixture because of the small values of the liquid water content. Therefore, only the effects of air on spherical droplets are considered and the two phases are solved separately. The flow field is turbulent as the Reynolds number based on the chord length is Re>5×105 for all velocity values. The Reynolds-averaged Navier–Stokes equations are employed to model incompressible turbulent airflow, wherein an instantaneous flow variable is decomposed into its mean and fluctuating components. This set of equations for the air phase (specified by the subscript a) can be written as follows:
·V¯a=0
(1)
V¯at+V¯a·V¯a=1ρaP¯+νa2V¯a·VaVa¯
(2)
T¯at+V¯a·T¯a=kρacp2T¯a·VaTa¯
(3)

where ρ, ν, cp, and k are the density, kinematic viscosity, specific heat, and thermal conductivity of the air phase. The quantities P¯,V¯,T¯ represent the mean pressure, velocity, and temperature, respectively, and V and T are the fluctuating velocity and temperature. The quantities VV¯ and VT¯ represent the specific Reynolds stress tensor and the specific turbulent heat fluxes. To calculate the components of the Reynolds stress tensor, an empirical closure known as a turbulence model is required. A very convenient model particularly useful in icing problems is the Spalart–Allmaras model, which is applied to all simulations [24].

2.2 Droplet Motion Modeling.

The governing equation of droplet trajectory can be derived through the integration of Newton's second law with the following assumptions: (i) mass and heat transfer between air and droplets are ignored, (ii) thermophysical properties of the droplets are constant, and (iii) coalescence and collision of the droplets is disregarded. The droplet continuity and momentum equations (specified by the subscript d) in the Eulerian framework can be expressed as follows:
αdt+·(αdVd)=0
(4)
(αdVd)t+(Vd·)(αdVd)=αdCdRed24K(VaVd)+(1ρaρd)αdgFr2
(5)
where K=ρad2/18μa is the droplet inertial parameter. The first term on the right side of Eq. (5) represents the drag force on the droplets, where Cd is the droplet drag coefficient and Red is the Reynolds number determined based on the relative velocity of air and the droplets. The second term on the right side reflects buoyancy and gravity forces, with Fr representing the Froude number based on the airfoil chord length and the freestream velocity. The volume fraction of the supercooled droplets can be defined as follows:
αd=ωρd
(6)
where ω is the liquid water content. The catch rate of the impinging droplets, which is normally called the local collection efficiency, β, is defined as the ratio of the surface collection rate of the water droplets to the freestream flow. Having the air flow velocity and droplet fields, the collection efficiency of the droplets onto the airfoil surface can be determined as follows:
β=αdVd·nV
(7)

where n is the surface normal vector.

When dealing with the SLD regime, additional dynamic effects such as droplet splashing, deformation, and break up, which influence the collection efficiency, should be considered. During the break up process, a large droplet is unstable and tends to break up into a number of smaller stable droplets due to the aerodynamic forces (e.g., shear force and pressure gradient). The critical Weber number (Wecr) for the onset of droplet breakup is as given below according to Pilch and Erdman [25]
Wecr=12(1+1.077(Oh)1.6)
(8)
where Oh is the dimensionless Ohnesorge number that relates the viscous forces to inertial and surface tension forces. A large droplet can also reach a critical condition at which its shape begins to deviate from spherical to an oblate disk. The droplet drag coefficient can be defined as follows according to Clift et al. [26]:
Cd=(1ϕ)Cd|spherical+ϕCd|disk
(9)
where ϕ is an eccentricity function depending on the value of the relative Weber number, We
ϕ=11(1+0.007Wer)6
(10)
The outcome of the impingement of droplets onto a surface comprises stick, spread, rebound, and splash processes or a combination of them. These processes may be characterized by the impact of Weber number according to Bai and Gosman [27]. The bouncing and splashing mechanisms are considered in our simulations because of their importance in so far as determining the local collection efficiency is concerned. Rebounding is a process, in which the impinging droplet rebounds off the surface after impact as a result of the air trapped between the droplet and the wetted surface. At high-impact energy, splashing of the droplets also occurs which creates numerous small ejected droplets, and ultimately reduces the amount of deposited water. The impact of bouncing-splashing on droplet momentum can be represented by a body force term, f, applied to the vicinity of the surface
f=VPVIδt
(11)
where VP is the postimpact velocity of the splashed droplets, VI is the impact velocity of the primary droplets, and δt is the collision contact time. The Mundo model is selected to determine the probability of splashing according to the following equation by Mundo et al. [28]:
(We)(Oh)2/5540R0.35
(12)

where R is the mean roughness height of the surface. Because the problem investigated is in the SLD regime, the droplets behave as rain drops falling with a terminal velocity. This terminal velocity should be accounted for in the modeling.

2.3 Glaze Ice Modeling.

Glaze icing zones constitute a thin film of water that grows above the iced surface and its velocity, uf, is presumed to be a function of the shear stress and gravity. From the momentum equation, the average film velocity can be expressed as follows:
u¯f(x,y)=hf2μfτwall(x,y)ρg3μfhf2
(13)
where x=(x1,x2) stands for the curvilinear coordinate system on the surface, y is normal to the surface, τwall is the shear stress from the air movement, and hf is the film height [7]. The resulting PDE formulation for mass conservation can be written as follows:
ρf(hft+·(u¯fhf))=m˙impm˙evapm˙ice
(14)
where m˙imp=ωβVd,,m˙evap, and m˙ice correspond to the mass flux rates due to impingement, evaporation, and solidification, respectively. The first term on the right side of Eq. (14) represents the impingement mass flux. The energy balance across the film can also be written as follows:
ρw((hfcpT̃f)t+·(u¯fhfcpT̃f))=ωβVd,(cpTd,+Vd22)+q˙iceq˙evapq˙convq˙rad+q˙cond
(15)

where q˙ice,q˙evap,q˙conv,q˙rad, and q˙cond are, respectively, the heat flux rates due to ice accretion, evaporation, convection, radiation, and conduction. The last two heat transfer mechanisms are zero in the current icing simulations. The quantity Vd2 is the square norm of the droplet impact velocity and T̃=TTref is the equilibrium temperature of the film, where Tref is a reference temperature whose value was set as the triple point of water in Kelvin. The three unknowns (hf,T̃f,m˙ice) can be computed with the aid of compatibility relations to enforce the following two conditions: (1) no water film forms when the equilibrium temperature is below the freezing point, (2) no ice accumulates if the film temperature exceeds 0 °C. It is worth of note that ice may yet form even if the local surface temperature is above freezing owing to the evaporation of the ice/water mixture. More details about the individual terms of Eqs. (14) and (15) can be found in Beaugendre [29].

2.4 Boundary Conditions.

To solve the governing equations for both the air and the droplets, Eqs. (1)(5), proper boundary conditions should be prescribed. A constant temperature boundary condition is adopted at the airfoil surface, which is also considered a no-slip wall. The far-field stream is considered as having constant velocity and temperature fields. When the flow (air-droplet) strikes the airfoil surface, the first layer formed on the surface is ice, and the unfrozen droplets produce a layer of water film on top of the ice. The water film is partially influenced by the shear forces and partially by gravitational forces. An angle of attack (α=4deg) has been prescribed for the flow field throughout the study as indicated in Fig. 1. When the flow (air-droplet) strikes the airfoil surface, the first layer formed on the surface is ice, and the unfrozen droplets produce a layer of water film on top of the ice as shown in Fig. 1. The water film is partially influenced by the shear forces and partially by gravitational forces. This combined effect is captured by Eq. (13). Additional equations pertaining to the mass and energy balances of the droplets are captured by Eqs. (14) and (15), respectively.

Fig. 1
(a) Flow around the airfoil showing the relevant parameters and the different heat fluxes for glaze ice accretion and (b) structured grid around the airfoil with a higher density in the vicinity of the airfoil surface
Fig. 1
(a) Flow around the airfoil showing the relevant parameters and the different heat fluxes for glaze ice accretion and (b) structured grid around the airfoil with a higher density in the vicinity of the airfoil surface
Close modal

2.5 Ice Roughness.

Turbulent boundary layers are highly sensitive to the ice roughness element because the law of the wall is modified if any roughness element protrudes through the viscous sublayer. Roughness affects the convective heat transfer rate, which consequently affects the ice accretion rate, especially under glaze icing conditions. Pais et al. [30] showed that surface roughness can also change the maximum liquid film thickness, the droplet diameter, and its departure characteristics. The empirical sand-grain correlation developed by Shin and Bond [31] is used in the current simulations, where the equivalent sand-grain roughness, ks, is a function of the freestream air temperature, the liquid water content, and the median volume diameter of the droplets as demonstrated by Eq. (16) below. The quantity c denotes the chord length
ks=8.05×104(0.5714+0.2457ω+1.2571ω2).(0.047T11.27).f(d)·c
(16)
where
f(d)={1d20μm1.6670.0333dd>20μm
(17)

2.6 Grid Convergence Study.

The entire computational domain is structurally gridded using the icemcfd module. To generate high-quality hexahedral grids, particularly those with acceptable element orthogonality, a Cartesian set of blocks is carved around the airfoil using a split block function available in the icemcfd module. Since the flow is viscous and turbulent, grid points have been clustered around the airfoil to better capture the boundary layer and surface movement as can be seen in Fig. 1. This figure also shows the flow field around the airfoil and the different fluxes of heat transfer affecting the airfoil surface. Conducting a mesh sensitivity study is always a good practice to ensure that the simulations provide a mathematically accurate solution. The process entails conducting successive levels of mesh refinement until two consecutive meshes generate nearly identical results for the given variable. According to Fig. 2, the percentage deviation in the drag coefficient between the first two coarse meshes is about 3.7%, while the percentage deviation can be seen to decrease as the number of elements increases. This has been done until the percentage deviation in the drag coefficient between two successive meshes has been reduced to within 1.4%. A grid convergence index (GCI) analysis is additionally used to measure the uncertainty of grid convergence and evaluate the discretization error as per Roache [32]. The GCI represents how much the solution would vary if the grid were further refined. The analysis must be carried out for three grid resolutions in order to compute two GCIs, for instance, from a fine to an intermediate mesh size GCI53 and from an intermediate to a coarse mesh size GCI31. For comparisons over three or more grids, the GCI can be defined as follows:
GCI=1.25|ϵ|(rp1)
(18)
where ϵ is the relative error for two subsequent values of the studied quantity. The quantity r1.3 is the ratio of grid refinement and p =2.26 denotes the order of convergence as given below:
p=ln(f3f2f2f1)/ln(r)
(19)
where f here is the investigated quantity (in this case Cd), and the subscripts represent the order of the element. Based on the provided GCI formulation, grid convergence indices of GCI31=12% and GCI53=7% can be computed. The following relationship should be examined to ensure that the asymptotic range of convergence is achieved:
GCI31rpGCI531
(20)
Fig. 2
Variation of the predicted drag coefficient against fivedifferent numbers of elements at V∞=67m/s, T∞=258 K, d=20μm, ω=1 g/m3, t=6 min, α=4deg
Fig. 2
Variation of the predicted drag coefficient against fivedifferent numbers of elements at V∞=67m/s, T∞=258 K, d=20μm, ω=1 g/m3, t=6 min, α=4deg
Close modal

To be in the asymptotic range, the left-hand side of the above equation should be close to 1. Results obtained for the left-hand side of the above equation in the current simulations were found to be 0.966, thus indicating good convergence.

3 Results and Discussion

3.1 Validation.

Numerical predictions from the current simulations were first compared with existing experimental data for the 2D NACA 0012 airfoil. All three sets of PDEs for airflow, droplet field, and water film were orderly and independently solved using the fluenticing module for the two-, three-, four-, and six-shot simulations, and were solved using FENSAP-ICE for quasi-steady simulations employing time steps of 5 s, 1 s, and 0.1 s. The total time duration for ice accretion has been selected to be 232 s for all simulations. The governing partial differential equations were linearized using the Newton–Galerkin method with the help of a second-order upwind formulation for numerical stability. The air/droplet flow was converged to a steady-state by an implicit time stepping scheme in order to enhance the stability of the linear system and reduce the total computational time. An arbitrary Lagrangian–Eulerian scheme has been used to account for the movement of the ice-accreted surface and mesh management. This mesh morphing method is comparably fast since the geometry's new shape is computed only by updating nodal position without changing the topology of the original grid.

In order to find a time-step that is sufficiently small to adequately simulate the time-dependent icing process in a quasi-steady-state context, we have divided the 232 s duration of simulation time into segments during which the steady-state solution of the conservation equations was carried out and the results of those simulations were used as initial conditions for advancing to the next time position. The magnitude of the time-step selected ranged from a least accurate two-shot simulation (where steady-state conditions were assumed to exist for 116 s at a time) to a most accurate 2320-shot simulation (where steady-state conditions were assumed to exist for only 0.1 s at a time). In-between these two extreme time steps, we have also performed a three-shot simulation (where steady-state conditions were assumed to exist for 77.33 s), a four-shot simulation (where steady-state conditions were assumed to exist for 58 s), a six-shot simulation (where steady-state conditions were assumed to exist for 38.67 s), a 46-shot simulation (where steady-state conditions were assumed to exist for 5 s), and a 232-shot simulation (where steady-state conditions were assumed to exist for 1 s). From Fig. 3, it can be seen that the ice profile reaches almost an identical shape when the time-step is equal to either 1 s or 0.1 s and that no discernible improvement in ice profile can be achieved with a larger number of shots. Table 1 shows the computational time versus the size of the time step, while Table 2 compares the angle (θ) between the ice horn and the chord line (with Point (0,0) being located at the tip of the leading edge) for different time steps. As can be seen from the table, the maximum percentage deviation in the value of the angle as compared to experimental data is reduced by about 20% as we reduce the value of the time-step from 116 s to either 1 s or 0.1 s since the 1 s and the 0.1 s simulation results are very close to each other. Said differently, a time-step of 0.1 s is in fact a sufficiently small time interval that should be able to reproduce the full time-dependent solution. Since the grid quality deteriorates after each shot, several numerical factors such as the Courant number and the number of iterations can be modified to improve the robustness and convergence of these simulations. Moreover, the remeshing strategy can be a useful alternative for the mesh morphing method to track high curvature zones around the ice horn and maintain optimal grid spacing that is usually influenced by large mesh deformations, particularly in the glaze ice accretion regime. Obviously, the main advantage of using fewer shots (or larger values for the time-step) is the reduction in the total computational time as summarized in Table 1. In this table, we can see the total computational time employing a Core i5 processor is drastically reduced from 98 h for the 2320-shot simulation to about half an hour for the six-shot simulation. As we will see later in the paper, the degree of detail in some of the other parameters that were computed was an important factor in selecting the 2320-shot simulation.

Fig. 3
Comparison of the final ice profile employing different time steps ranging from a two-shot (Δt=116 s) simulation to a 2320-shot (Δt=0.1 s) simulation at V∞=102.8 m/s, T∞=262 K, ω=1 g/m3, d=20μm, t=232 s, α=4deg. The light dotted line represents the experimental results reported by Ozcer et al. [33].
Fig. 3
Comparison of the final ice profile employing different time steps ranging from a two-shot (Δt=116 s) simulation to a 2320-shot (Δt=0.1 s) simulation at V∞=102.8 m/s, T∞=262 K, ω=1 g/m3, d=20μm, t=232 s, α=4deg. The light dotted line represents the experimental results reported by Ozcer et al. [33].
Close modal
Table 1

Computational time versus time step size

SchemeTime interval, sComputational time, hour
2 shots116 0.1
3 shots77 0.25
4 shots58 0.3
6 shots38 0.5
46 shots5 3.8
232 shots1 19
2320 shots0.1 98
SchemeTime interval, sComputational time, hour
2 shots116 0.1
3 shots77 0.25
4 shots58 0.3
6 shots38 0.5
46 shots5 3.8
232 shots1 19
2320 shots0.1 98
Table 2

Comparison of the angle that the Ice Horn makes with the chord length for different time steps as compared to experimental data (experimental angle θ59.2deg) of Ozcer et al. [33]

SchemeIce horn angle (θ)Percentage deviation
2 shots37.1 deg37%
3 shots37.5 deg36%
4 shots41.3 deg30%
6 shots42.2 deg29%
46 shots45.7 deg22%
232 shots49.1 deg17%
2320 shots49.5 deg16%
SchemeIce horn angle (θ)Percentage deviation
2 shots37.1 deg37%
3 shots37.5 deg36%
4 shots41.3 deg30%
6 shots42.2 deg29%
46 shots45.7 deg22%
232 shots49.1 deg17%
2320 shots49.5 deg16%

The film thickness, the local collection efficiency, and relevant ice morphology are depicted in Fig. 4, which shows that the ice accretion rate is directly dependent on the local collection efficiency. As can be seen, time steps of 58 s, 5 s, 1 s, and 0.1 s were chosen for the simulations for comparison among themselves as well as with the experimental data reported in Ozcer et al. [33], resulting in a reasonable prediction of the ice shape. The horn discrepancy observed in Fig. 4(c) may be attributed to the mesh morphing mechanism. This scheme inevitably enlarges the wall grids and their aspect ratio as ice accretes by virtue of fixed mesh topology between time positions. After a certain number of automatic shots, manual or automatic remeshing is therefore required to maintain proper mesh and solution quality. It is also worthwhile to note that the whole spectrum of droplets in a cloud is represented only by the median volume diameter of the droplets (d=20μm), which may result in some deviations. The glaze icing conditions, which have been considered are listed in Table 3.

Fig. 4
Variation of the water film thickness, local collection efficiency, and ice profile shape employing different timestepscompared with the experimental results reported by Ozcer et al. [33] at V∞=102.8 m/s, T∞=262 Km/s, ω=1 g/m3, d=20μm, t=232 s, α=4deg
Fig. 4
Variation of the water film thickness, local collection efficiency, and ice profile shape employing different timestepscompared with the experimental results reported by Ozcer et al. [33] at V∞=102.8 m/s, T∞=262 Km/s, ω=1 g/m3, d=20μm, t=232 s, α=4deg
Close modal
Table 3

Numerical values of the environmental parameters used for simulation of glaze ice accretion

V (m/s)T (K)ω (g/m3)d(μm)α
75, 102.8, 150262, 2680.5, 1, 220, 50, 100, 2004 deg
V (m/s)T (K)ω (g/m3)d(μm)α
75, 102.8, 150262, 2680.5, 1, 220, 50, 100, 2004 deg

3.2 Ice Accretion Versus Free Stream Velocity (V) and Liquid Water Content (or ω).

The thicknesses of the water film and ice layer strongly depend on the convective heat transfer rate. Results shown in Fig. 5 demonstrate that the film thickness is maximum at the stagnation point (y0.01m), where a vast number of water droplets impinge on the surface, and where the convective heat transfer rate is minimum. Thus, a substantial amount of water remains unfrozen and the ice thickness remains small. During the glaze icing regime, the water film on the upper surface is composed of impinging droplets, the remaining water, and the runback water. It is assumed that the runback water is shed by the combined impact of the aerodynamic and gravitational forces on the lower surface of the airfoil. In addition, higher LWC causes an expansion of the film/ice-covered area over the leading edge according to our data. It is worthy of note that once the LWC becomes equal to 2g/m3, the film thickness starts to develop another peak located just after the horn (around y0.02m), where the heat transfer rate is not large enough to convert the extra incoming water droplets to ice. This can be seen in Fig. 5(a). Finally, the key feature of glaze ice accretion is the simultaneous formation of liquid and ice on the airfoil surface. This can be clearly seen in Fig. 5.

Fig. 5
Variation of the water film thickness, the ice profile thickness, and the convective heat flux for a 2320-shot simulation (Δt=0.1 s) as a function of the LWC at V∞=102.8 m/s, T∞=262K, d=20μm, t=232 s, α=4deg
Fig. 5
Variation of the water film thickness, the ice profile thickness, and the convective heat flux for a 2320-shot simulation (Δt=0.1 s) as a function of the LWC at V∞=102.8 m/s, T∞=262K, d=20μm, t=232 s, α=4deg
Close modal

While temperature is a critical factor in influencing the shape of the final ice profile, the shape also depends on the freestream velocity (V) and the liquid water content (ω) as can be observed for different time steps in Fig. 6. This has also been confirmed by the experimental data of Shin and Bond [34]. The ice profile becomes virtually identical for all cases as we approach a true quasi-steady situation where a sufficiently small time interval (step) is used (e.g., the 1 s and 0.1 s simulations). Overall, any decrease in the local collection efficiency due to the lower liquid water content is compensated for by the increase in the freestream velocity at a constant angle of attack. The analysis also reveals a 12% or 8.5% increase in the maximum collection efficiency, βmax, depending on whether the value of the LWC has been doubled from 0.5 g/m3 to 1 g/m3, or the value of the freestream velocity has been doubled from 75 m/s to 150 m/s, respectively. Thus, the analysis reveals that the LWC has a greater influence on the overall ice buildup than does the freestream velocity. An ice horn may arise when the airfoil becomes exposed to conditions conducive to those associated with glaze ice accretion. The core aspects that determine flow separation are the thickness of the horn, its orientation angle, and its position on the airfoil surface. These are all highly influenced by the freestream velocity and/or the liquid water content. For example, at a freestream velocity of 75 m/s as shown in Figs. 6(a), 6(d), and 6(g), “streamwise ice” is totally converted to the horn-structured glaze ice (sometimes known as “horn ice”) when the LWC is doubled. The term “streamwise ice” is a common term in ice accretion terminology that refers to ice accreting along the curvilinear surface of the airfoil (i.e., parallel to the surface) [35]. Based on Figs. 6(d), 6(e), and 6(f) at a constant value for the LWC, the ice horn shrinks and its location seems to move downstream of the airfoil leading edge as the freestream velocity increases, forming a “streamwise ice” layer on the airfoil surface. This type of ice accretion typically occurs when all the impinging droplets are not evaporated, thus causing the formation of runback droplets. The freezing of even a small proportion of supercooled droplets may release enough latent heat for the remainder of the droplets to become warm, stay as liquid, and flow over the surface as runback, where they could be converted to ice farther downstream. Whereas this type of ice accretion is usually associated with supercooled large droplet conditions, it can exist for all droplet sizes resulting in larger impingement areas in the aft region of the airfoil.

Fig. 6
Variation of the ice profile horn shape and its location obtained for different time steps (Δt=58 s,  5 s,  1 s,  0.1 s) as a function of the freestream velocity and LWC at T∞=262K, d=20μm, t=232 s, α=4deg. Measured data are obtained at V∞=67m/s and t=360 s from Shin and Bond [34].
Fig. 6
Variation of the ice profile horn shape and its location obtained for different time steps (Δt=58 s,  5 s,  1 s,  0.1 s) as a function of the freestream velocity and LWC at T∞=262K, d=20μm, t=232 s, α=4deg. Measured data are obtained at V∞=67m/s and t=360 s from Shin and Bond [34].
Close modal

Liquid water content affects both the type of ice accreting and the rate at which it accretes. With all other variables relevant to ice accretion remaining constant, LWC has a direct relationship with the amount and number of impinging droplets. According to Figs. 6(a), 6(d), and 6(g), the LWC appears to have a significant influence on the horn shape and especially on the emergence of a lower horn. Larger values of the LWC increase the likelihood of glaze ice accretion. On the other hand, at lower LWC such as at ω=0.5g/m3, the ice shape resembles the original airfoil profile in general and can be categorized “streamwise ice.”

3.3 Ice Accretion Versus Droplet Diameter (d).

Given the fact that the droplet mass and aerodynamic forces are, respectively, in proportion to the cube and square of the droplet diameter, large droplets are less influenced by the airflow in the vicinity of the airfoil surface. This fact allows these large droplets to strike the airfoil farther at the aft region. Ice accretion from larger droplets is more prone to grow into shapes that can disrupt the flow of air over the airfoil. Variation of the ice structure with the freestream velocity and droplet diameter is shown in Fig. 7. Compared with the previous cases shown in Fig. 6, larger droplets are more likely to impinge on the airfoil because of their non-negligible terminal velocity, whereas tiny droplets are typically reflected. Thus, the ice is thicker and more irregular in structure in the SLD regime, with the ice horn forming even at lower freestream velocities (e.g., V=75m/s). As a result, significant ice accretion is more likely to occur in the SLD regime. As can be seen in Fig. 7, there is an upper horn that is more discernible when d=50,100,200μm due to an increasing trend in the local collection efficiency as the horn curvature increases. As the freestream velocity increases, a significant amount of runback and secondary impingement is produced as illustrated in Figs. 7(c), 7(f), and 7(i). Secondary impingement refers to the action of unfrozen water droplets that can no longer adhere to the accreted ice and thus are more likely to re-impinge on the surface farther downstream, thereby expanding the area upon which ice accretes. Examination of the shape and thickness of the ice profiles provided by Figs. 6 and 7 reveals that the secondary impingement is generally caused by either a large droplet diameter or a large LWC, with the former appearing to be more significant than the latter. Figure 8 shows that the convective heat flux decreases when the freestream velocity increases for all droplet sizes. As predicted in Fig. 9, the droplet size increases the rate of ice accretion by increasing the local value of the collection efficiency and by broadening the impingement area farther downstream of the airfoil. Because of the larger inertia (or larger mass) associated with larger droplets, it is more difficult for the airflow to change the droplet path of larger droplets.

Fig. 7
Variation of the ice profile shape in the SLD regime obtained for different time steps (Δt=58 s,  5 s,  1 s,  0.1 s) as a function of the freestream velocity and droplet size at T∞=262K, ω=1 g/m3, t=232 s, α=4deg
Fig. 7
Variation of the ice profile shape in the SLD regime obtained for different time steps (Δt=58 s,  5 s,  1 s,  0.1 s) as a function of the freestream velocity and droplet size at T∞=262K, ω=1 g/m3, t=232 s, α=4deg
Close modal
Fig. 8
Variation of the convective heat flux obtained by the 2320-shot simulation (Δt=0.1 s) as a function of droplet size and freestream velocity at T∞=262K, ω=1 g/m3, t=232 s, α=4deg
Fig. 8
Variation of the convective heat flux obtained by the 2320-shot simulation (Δt=0.1 s) as a function of droplet size and freestream velocity at T∞=262K, ω=1 g/m3, t=232 s, α=4deg
Close modal
Fig. 9
Distribution of the local collection efficiency obtained by the 2320-shot simulation (Δt=0.1 s) as a function of droplet size and freestream velocity at T∞=262K, ω=1 g/m3, t=232 s, α=4deg
Fig. 9
Distribution of the local collection efficiency obtained by the 2320-shot simulation (Δt=0.1 s) as a function of droplet size and freestream velocity at T∞=262K, ω=1 g/m3, t=232 s, α=4deg
Close modal

3.4 Ice Accretion Versus Free Stream Temperature (T).

Free stream temperature is an important variable that affects the final ice profile shape since it influences the heat transfer rate of droplets as they impinge on the airfoil surface. A higher freestream temperature causes a reduction in the amount of ice accreted on the airfoil and increases the likelihood of droplet runback [36]. The behavior of the runback droplets changes the shape of the ice profile, which in turn, affects the droplet impingement rate as well as the water film around the airfoil. Son et al. [37] observed that a secondary effect of the freestream temperature is the reduction of the mean value of the liquid water content when the freestream temperature decreases since the clouds are likely to contain more ice than liquid water at lower temperatures. In the present investigation, we have elected to study the effect of the freestream temperature by examining the behavior of the ice profile shape and the convective heat flux. Figure 10 displays the ice profile shape as a function of the freestream temperature (T) and droplet diameter (d). Examining the evolution of the ice profile shape in the figure reveals that the thicker ice layer at the lower temperature T=262K transforms to a thinner ice layer that is more spread out at the higher temperature T=268K where the ice profile shape extends on the top and lower sides of the airfoil surface. Examination of the distribution of the convective heat flux as a function of the freestream temperature and droplet diameter as shown in Figs. 10 and 11 reveal a profound effect of both variables. As we have seen earlier, a higher freestream temperature thins and spreads the ice profile shape. This in turn impacts the convective heat flux for the reasons discussed earlier in the paper.

Fig. 10
Variation of the ice profile shape in the SLD regime obtained for different time steps (Δt=58 s,  5 s,  1 s,  0.1 s) as a function of the freestream temperature and droplet size at V∞=75 m/s, ω=1 g/m3, t=232 s, α=4deg
Fig. 10
Variation of the ice profile shape in the SLD regime obtained for different time steps (Δt=58 s,  5 s,  1 s,  0.1 s) as a function of the freestream temperature and droplet size at V∞=75 m/s, ω=1 g/m3, t=232 s, α=4deg
Close modal
Fig. 11
Comparison of the convection heat flux obtained for different time steps (Δt=58 s,  5 s,  1 s,  0.1 s) as a function of the droplet size and freestream temperature at V∞=75 m/s, ω=1 g/m3, t=232 s, α=4deg
Fig. 11
Comparison of the convection heat flux obtained for different time steps (Δt=58 s,  5 s,  1 s,  0.1 s) as a function of the droplet size and freestream temperature at V∞=75 m/s, ω=1 g/m3, t=232 s, α=4deg
Close modal

Figure 11 provides a graphical display of the convective heat flux in terms of the numerical value of the time-step employed in the quasi-steady simulations. It is clear that the 1 s and 0.1 s simulations are able to capture more intricate spatial variations of the convective heat flux that are not fully captured by the 58 s simulations. With the 1 s and 0.1 s simulations of the convective heat flux being virtually similar, our selection of the more conservative 0.1 s simulations for the work presented here is now validated by the convective heat flux data in addition to the earlier validation that has been presented through the ability to reproduce the shape of the ice profile. According to Fig. 11, at T=262K, when the size of the droplet increases from 100μm to 200μm, the heat transfer rate decreases, especially near the ice horn region. Lynch and Khodadoust [9] stated that small droplets tend to hit the airfoil surface near the leading edge while larger droplets tend to have their impact on the surface farther downstream because of their larger inertia. Thus, not all of the impinging droplets might freeze at the leading edge and instead partially splash and partially freeze on the airfoil surface at locations farther away from the ice horn region. Last but not least, Fig. 12 shows the temperature evolution at the ice surface, whereby the surface temperature is increased as time passes. Liu et al. [17] stated that the higher temperature at the stagnation point can be attributed to the fact that as more droplets impinge on the stagnation point, more heat needs to be released. This heat causes the temperature to continue to rise until it reaches the freezing point of water.

Fig. 12
Temperature-position-time history on the iced-up surface of the airfoil obtained using a 0.1 s time-step simulation (Δt=0.1 s) at V∞=75 m/s, T∞=262K, d=50μm, ω=1 g/m3, α=4deg, t=232 s
Fig. 12
Temperature-position-time history on the iced-up surface of the airfoil obtained using a 0.1 s time-step simulation (Δt=0.1 s) at V∞=75 m/s, T∞=262K, d=50μm, ω=1 g/m3, α=4deg, t=232 s
Close modal

4 Conclusions

A quasi-steady simulation of glaze ice accretion on an airfoil in turbulent flow of atmospheric air in the SLD regime has been conducted. The essence of a successful and valid quasi-steady simulation is to assume that steady-state conditions exist for sufficiently small time intervals before advancing to the next time position. The litmus test for the success of the quasi-steady simulation is the ability to reproduce the full time-dependent solution or reproduce experimental data within a very small error band. Results of the ice profile shape presented here were able to capture the glaze ice shape and especially capture the ice horn. We have elected to simulate glaze ice accretion during a 232-second duration by dividing it into 2320 time segments such that steady-state conditions are assumed for a time duration of only 0.1 s at a time. In glaze ice accretion, the shear stress and heat transfer effects control the behavior of the droplet runback and the rate of ice accretion. Results of both the 1 s and 0.1 s simulations appear virtually identical for the slate of relevant variables that have been selected. This validates the notion that the time-step selected for the simulations is indeed” sufficiently small.” Examination of the variation of the angle that the ice horn makes with the airfoil chord line demonstrated a 20% improvement in angle prediction when the time-step is reduced from 116 s to 0.1 s. A mesh morphing method has been employed in the computations to accelerate tracking the moving boundaries and to regenerate the updated mesh before each transition to the next time position. It is important to note here, however, that this technique is suitable only for the lower shot simulations (large time steps). For simulations involving a very large number of shots (very small time steps), remeshing will be required after a certain number of shots have been performed in order to ensure a more accurate application of the quasi-steady approach adopted here. The effect of environmental factors such as freestream velocity and temperature, liquid water content, and droplet size has been investigated and their effect on several quantities such as the local collection efficiency, film thickness, convective heat transfer rate, and ice profile shape has been investigated. Simulations were performed on a two-dimensional NACA 0012 airfoil. Simulation results indicated that increasing the LWC completely converts the so-called “streamwise ice” into a horn-structured glaze ice (sometimes known as” horn ice”). Results also show that the ice horn shrinks in the airfoil leading edge zone and spreads farther downstream of the leading edge, forming a spanwise ridge of ice on the upper airfoil surface as the freestream velocity increases. The analysis also reveals a 12% or 8.5% increase in the maximum collection efficiency, βmax, depending on whether the value of the LWC has been doubled from 0.5 g/m3 to 1 g/m3, or the value of the freestream velocity has been doubled from 75 m/s to 150 m/s, respectively. Thus, the analysis reveals that the LWC has a greater influence on the overall ice buildup than does the freestream velocity. Results also show a close link between the convective heat transfer rate on the one hand and the film and ice thicknesses on the other. Results pertaining to the effect of droplet diameter show a thicker ice with an irregular structure in the SLD regime, with ice horns forming even at lower freestream velocities. At higher freestream temperatures just below the freezing point of water, the convective heat transfer rate decreases, which results in a lower ice accumulation rate.

Funding Data

  • U.S. Department of Energy's Office of Energy Efficiency and Renewable Energy (EERE) under the Advanced Manufacturing Office (Award No. DE-EE0009729; Funder ID: 10.13039/100000015; 10.13039/100006134).

  • Department of Mechanical and Aerospace Engineering at the University of Florida.

Nomenclature

c =

chord length (m)

Cd =

drag coefficient

cp =

specific heat (J kg–1 K–1)

g =

gravitational acceleration (m s–2)

h =

film height (m)

k =

thermal conductivity (W m–1 K–1)

K =

droplet inertial parameter, s

ks =

Equivalent sand-grain roughness, (m)

m˙ =

mass flux rate (kg m–2 s–1)

n =

surface normal vector

p =

order of grid convergence

P¯ =

mean pressure (Pa)

P =

fluctuating pressure (Pa)

q˙ =

heat flux rate (J m–2 s–1)

r =

mesh refinement ratio

R =

surface roughness (mm)

t =

time (s)

T̃ =

film equilibrium temperature (°C)

T¯ =

mean temperature (K)

T =

fluctuating temperature (K)

u¯ =

film average velocity (ms–1)

V =

Velocity vector (ms–1)

V¯ =

mean velocity (ms–1)

V =

fluctuating velocity (ms–1)

Greek Symbols
α =

angle of attack

αd =

volume fraction of droplet

ϵ =

relative error

θ =

angle of ice horn (deg)

μ =

dynamic viscosity (kg m–1 s–1)

ν =

kinematic viscosity (m2 s–1)

ρ =

density (kg m–3)

σ =

droplet surface tension force (N m–1)

τ=

shear stress (Pa)

ω =

liquid water content (g m–3)

Dimensionless Groups
Fr =

Froude number, V/gc

Oh =

Ohnesorge number, μ/ρdσ

Re =

Reynolds number, ρVd/μ

Pr =

Prandtl number, ν/α

Wer =

Weber number, ρ(VaVd)2d/σ

Superscripts and Subscripts
a =

air

d =

droplet

f =

film

I =

impact

P =

postimpact

r =

relative

=

free stream value

References

1.
Messinger
,
B. L.
,
1953
, “
Equilibrium Temperature of an Unheated Icing Surface as a Function of Air Speed
,”
J. Aeronaut. Sci.
,
20
(
1
), pp.
29
42
.10.2514/8.2520
2.
Ruff
,
G. A.
, and
Berkowitz
,
B. M.
,
1990
, “
Users Manual for the NASA Lewis Ice Accretion Prediction Code (LEWICE)
,” NASA, Washington, DC, Report No. 185129.
3.
Hedde
,
T.
, and
Guffond
,
D.
,
1995
, “
ONERA Three-Dimensional Icing Model
,”
AIAA J.
,
33
(
6
), pp.
1038
1045
.10.2514/3.12795
4.
Myers
,
T. G.
,
2001
, “
Extension to the Messinger Model for Aircraft Icing
,”
AIAA J.
,
39
(
2
), pp.
211
218
.10.2514/2.1312
5.
Sherif
,
S. A.
,
Pasumarthi
,
N.
, and
Bartlett
,
C.
,
1997
, “
A Semi-Empirical Model for Heat Transfer and Ice Accretion on Aircraft Wings in Supercooled Clouds
,”
Cold Regions Sci. Technol.
,
26
(
3
), pp.
165
179
.10.1016/S0165-232X(97)00021-9
6.
Naterer
,
G.
,
2003
, “
Temperature Gradient in the Unfrozen Liquid Layer for Multiphase Energy Balance With Incoming Droplets
,”
ASME J. Heat Transfer-Trans. ASME
,
125
(
1
), pp.
186
189
.10.1115/1.1532015
7.
Bourgault
,
Y.
,
Beaugendre
,
H.
, and
Habashi
,
W. G.
,
2000
, “
Development of a Shallow-Water Icing Model in FENSAP-ICE
,”
J. Aircr.
,
37
(
4
), pp.
640
646
.10.2514/2.2646
8.
Zhu
,
C.
,
Fu
,
B. I. N.
,
Sun
,
Z.
, and
Zhu
,
C.
,
2012
, “
3D Ice Accretion Simulation for Complex Configuration Basing on Improved Messinger Model
,”
Int. J. Mod. Phys. Conf. Ser.
,
19
, pp.
341
350
.10.1142/S2010194512008938
9.
Lynch
,
F. T.
, and
Khodadoust
,
A.
,
2001
, “
Effects of Ice Accretions on Aircraft Aerodynamics
,”
Prog. Aerosp. Sci.
,
37
(
8
), pp.
669
767
.10.1016/S0376-0421(01)00018-5
10.
Ashenden
,
R.
,
Lindberg
,
W.
,
Marwitz
,
J. D.
, and
Hoxie
,
B.
,
1996
, “
Airfoil Performance Degradation by Supercooled Cloud, Drizzle, and Rain Drop Icing
,”
J. Aircraft
,
33
(
6
), pp.
1040
1046
.10.2514/3.47055
11.
Honsek
,
R.
,
Habashi
,
W. G.
, and
Aubé
,
M. S.
,
2008
, “
Eulerian Modeling of in-Flight Icing Due to Supercooled Large Droplets
,”
J. Aircr.
,
45
(
4
), pp.
1290
1296
.10.2514/1.34541
12.
Kim
,
I.
,
Bachchan
,
N.
, and
Peroomian
,
O.
,
2016
, “
Supercooled Large Droplet Modeling for Aircraft Icing Using an Eulerian–Eulerian Approach
,”
J. Aircr.
,
53
(
2
), pp.
487
500
.10.2514/1.C033506
13.
Wang
,
C.
,
Chang
,
S.
,
Leng
,
M.
,
Wu
,
H.
, and
Yang
,
B.
,
2016
, “
A Two-Dimensional Splashing Model for Investigating Impingement Characteristics of Supercooled Large Droplets
,”
Int. J. Multiphase Flow
,
80
, pp.
131
149
.10.1016/j.ijmultiphaseflow.2015.12.005
14.
Raj
,
L. P.
,
Lee
,
J.
, and
Myong
,
R.
,
2019
, “
Ice Accretion and Aerodynamic Effects on a Multi-Element Airfoil Under SLD Icing Conditions
,”
Aerosp. Sci. Technol.
,
85
, pp.
320
333
.10.1016/j.ast.2018.12.017
15.
Stebbins
,
S. J.
,
Loth
,
E.
,
Broeren
,
A. P.
, and
Potapczuk
,
M.
, Nov,
2019
, “
Review of Computational Methods for Aerodynamic Analysis of Iced Lifting Surfaces
,”
Prog. Aerosp. Sci.
,
111
, p.
100583
.10.1016/j.paerosci.2019.100583
16.
Cao
,
Y.
,
Huang
,
J.
, and
Yin
,
J.
,
2016
, “
Numerical Simulation of Three-Dimensional Ice Accretion on an Aircraft Wing
,”
Int. J. Heat Mass Transfer
,
92
, pp.
34
54
.10.1016/j.ijheatmasstransfer.2015.08.027
17.
Liu
,
T.
,
Qu
,
K.
,
Cai
,
J.
, and
Pan
,
S.
,
2019
, “
A Three-Dimensional Aircraft Ice Accretion Model Based on the Numerical Solution of the Unsteady Stefan Problem
,”
Aerosp. Sci. Technol.
,
93
, p.
105328
.10.1016/j.ast.2019.105328
18.
Gori
,
G.
,
Parma
,
G.
,
Zocca
,
M.
, and
Guardone
,
A.
,
2018
, “
Local Solution to the Unsteady Stefan Problem for in-Flight Ice Accretion Modeling
,”
J. Aircr.
,
55
(
1
), pp.
251
262
.10.2514/1.C034412
19.
Cao
,
Y.
,
Li
,
G.
, and
Song
,
D.
,
2021
, “
Numerical Simulation of Melting of Ice Accreted on an Airfoil
,”
Aerosp. Sci. Technol.
,
119
, p.
107223
.10.1016/j.ast.2021.107223
20.
Shad
,
A.
, and
Sherif
,
S. A.
,
2022
, “
Numerical Modeling of Atmospheric Rime Ice Accretion on an Airfoil Using an Eulerian Approach
,”
ASME J. Fluids Eng.
,
144
(
6
), p.
061112
.10.1115/1.4054048
21.
Li
,
S.
, and
Paoli
,
R.
,
2019
, “
Modeling of Ice Accretion Over Aircraft Wings Using a Compressible OpenFOAM Solver
,”
Int. J. Aerosp. Eng.
,
2019
, pp.
1
11
.10.1155/2019/4864927
22.
Ozcer
,
I.
,
Switchenko
,
D.
,
Baruzzi
,
G. S.
, and
Chen
,
J.
,
2019
, “
Multi-Shot Icing Simulations With Automatic Re-Meshing
,”
SAE Paper No. 2019-01-1956.
23.
Aliaga
,
C. N.
,
Aubé
,
M. S.
,
Baruzzi
,
G. S.
, and
Habashi
,
W. G.
,
2011
, “
FENSAP-ICE-Unsteady: Unified in-Flight Icing Simulation Methodology for Aircraft, Rotorcraft, and Jet Engines
,”
J. Aircr.
,
48
(
1
), pp.
119
126
.10.2514/1.C000327
24.
Spalart
,
P.
, and
Allmaras
,
S.
,
1992
, “
A One-Equation Turbulence Model for Aerodynamic Flows
,”
AIAA
Paper No.
92
0439
.10.2514/6.1992-439
25.
Pilch
,
M.
, and
Erdman
,
C.
,
1987
, “
Use of Breakup Time Data and Velocity History Data to Predict the Maximum Size of Stable Fragments for Acceleration-Induced Breakup of a Liquid Drop
,”
Int. J. Multiphase Flow
,
13
(
6
), pp.
741
757
.10.1016/0301-9322(87)90063-2
26.
Clift
,
R.
,
Grace
,
J. R.
, and
Weber
,
M. E.
,
1978
,
Bubbles, Drops, and Particles
, 3rd ed.,
Academic Press
,
New York
.
27.
Bai
,
C.
, and
Gosman
,
A.
,
1995
, “
Development of Methodology for Spray Impingement Simulation
,”
SAE
Paper No. 950283.10.4271/950283
28.
Mundo
,
C.
,
Sommerfeld
,
M.
, and
Tropea
,
C.
,
1995
, “
Droplet-Wall Collisions: Experimental Studies of the Deformation and Breakup Process
,”
Int. J. Multiphase Flow
,
21
(
2
), pp.
151
173
.10.1016/0301-9322(94)00069-V
29.
Beaugendre
,
H.
,
2003
, “
A PDE-Based 3D Approach to in-Flight Ice Accretion
,” Ph.D. thesis,
Department of Mechanical Engineering, McGill University
,
Montreal, QC, Canada
.
30.
Pais
,
M.
,
Chow
,
L.
, and
Mahefkey
,
E.
,
1992
, “
Surface Roughness and Its Effects on the Heat Transfer Mechanism in Spray Cooling
,”
ASME J. Heat Transfer-Trans. ASME
,
114
(
1
), pp.
211
219
.10.1115/1.2911248
31.
Shin
,
J.
, and
Bond
,
T.
,
1992
, “
Results of an Icing Test on a NACA 0012 Airfoil in the NASA Lewis Icing Research Tunnel
,” NASA, Washington, DC, Report No. 105374.
32.
Roache
,
P. J.
,
1994
, “
Perspective: A Method for Uniform Reporting of Grid Refinement Studies
,”
ASME J. Fluids Eng.
,
116
(
3
), pp.
405
413
.10.1115/1.2910291
33.
Ozcer
,
I. A.
,
Baruzzi
,
G. S.
,
Reid
,
T.
,
Habashi
,
W. G.
,
Fossati
,
M.
, and
Croce
,
G.
,
2011
, “
Fensap-Ice: Numerical Prediction of Ice Roughness Evolution, and Its Effects on Ice Shapes
,”
SAE
Paper No. 2011-38-0024.10.4271/2011-38-0024
34.
Shin
,
J.
, and
Bond
,
T. H.
,
1992
, “
Experimental and Computational Ice Shapes and Resulting Drag Increase for a NACA 0012 Airfoil
,” NASA, Washington, DC, Report No. 105743.
35.
Bragg
,
M. B.
,
Broeren
,
A. P.
, and
Blumenthal
,
L. A.
,
2005
, “
Iced-Airfoil Aerodynamics
,”
Prog. Aerosp. Sci.
,
41
(
5
), pp.
323
362
.10.1016/j.paerosci.2005.07.001
36.
Gultepe
,
I.
, and
Isaac
,
G.
,
1997
, “
Liquid Water Content and Temperature Relationship From Aircraft Observations and Its Applicability to GCMs
,”
J. Clim.
,
10
(
3
), pp.
446
452
.10.1175/1520-0442(1997)010<0446:LWCATR>2.0.CO;2
37.
Son
,
C.
,
Oh
,
S.
, and
Yee
,
K.
,
2013
, “
Quantitative Investigation Into the Relationship Between Ice Accretion Shape and Ambient Conditions
,”
Trans. Jpn. Soc. Aeronaut. Space Sci.
,
56
(
4
), pp.
175
186
.10.2322/tjsass.56.175