## 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.5\u2009\u2009g/m3$ to $1\u2009\u2009g/m3$, or the value of the freestream velocity has been doubled from $75\u2009\u2009m/s$ to $150\u2009\u2009m/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\u2009\mu 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\u2009\mu 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.34\u2009\u2009cm$ 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.

*a*) can be written as follows:

where *ρ*, *ν*, *c _{p}*, and

*k*are the density, kinematic viscosity, specific heat, and thermal conductivity of the air phase. The quantities $P\xaf,\u2009V\xaf,\u2009T\xaf$ represent the mean pressure, velocity, and temperature, respectively, and $V\u2032$ and $T\u2032$ are the fluctuating velocity and temperature. The quantities $V\u2032V\u2032\xaf$ and $V\u2032T\u2032\xaf$ 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.

*d*) in the Eulerian framework can be expressed as follows:

*C*is the droplet drag coefficient and Re

_{d}*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}*ω*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:

where **n** is the surface normal vector.

_{cr}) for the onset of droplet breakup is as given below according to Pilch and Erdman [25]

**f**, applied to the vicinity of the surface

*δ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]:

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.

*y*is normal to the surface, $\tau wall$ is the shear stress from the air movement, and

*h*is the film height [7]. The resulting PDE formulation for mass conservation can be written as follows:

_{f}where $q\u02d9ice,\u2009q\u02d9evap,\u2009q\u02d9conv,\u2009q\u02d9rad$, and $q\u02d9cond$ 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 $\Vert Vd\Vert 2$ is the square norm of the droplet impact velocity and $T\u0303=T\u2212Tref$ is the equilibrium temperature of the film, where *T*_{ref} is a reference temperature whose value was set as the triple point of water in Kelvin. The three unknowns ($hf,T\u0303f,m\u02d9ice$) 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 ($\alpha =4\u2009deg$) 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.

### 2.5 Ice Roughness.

*k*, 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

_{s}*c*denotes the chord length

### 2.6 Grid Convergence Study.

_{53}and from an intermediate to a coarse mesh size GCI

_{31}. For comparisons over three or more grids, the GCI can be defined as follows:

*ϵ*is the relative error for two subsequent values of the studied quantity. The quantity $r\u22481.3$ is the ratio of grid refinement and

*p*=

*2.26 denotes the order of convergence as given below:*

*f*here is the investigated quantity (in this case

*C*), 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:

_{d}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.

Scheme | Time interval, s | Computational time, hour |
---|---|---|

2 shots | 116 | $\u2248$ 0.1 |

3 shots | 77 | $\u2248$ 0.25 |

4 shots | 58 | $\u2248$ 0.3 |

6 shots | 38 | $\u2248$ 0.5 |

46 shots | 5 | $\u2248$ 3.8 |

232 shots | 1 | $\u2248$ 19 |

2320 shots | 0.1 | $\u2248$ 98 |

Scheme | Time interval, s | Computational time, hour |
---|---|---|

2 shots | 116 | $\u2248$ 0.1 |

3 shots | 77 | $\u2248$ 0.25 |

4 shots | 58 | $\u2248$ 0.3 |

6 shots | 38 | $\u2248$ 0.5 |

46 shots | 5 | $\u2248$ 3.8 |

232 shots | 1 | $\u2248$ 19 |

2320 shots | 0.1 | $\u2248$ 98 |

Scheme | Ice horn angle (θ) | Percentage deviation |
---|---|---|

2 shots | 37.1 deg | $\u2248$37% |

3 shots | 37.5 deg | $\u2248$36% |

4 shots | 41.3 deg | $\u2248$30% |

6 shots | 42.2 deg | $\u2248$29% |

46 shots | 45.7 deg | $\u2248$22% |

232 shots | 49.1 deg | $\u2248$17% |

2320 shots | 49.5 deg | $\u2248$16% |

Scheme | Ice horn angle (θ) | Percentage deviation |
---|---|---|

2 shots | 37.1 deg | $\u2248$37% |

3 shots | 37.5 deg | $\u2248$36% |

4 shots | 41.3 deg | $\u2248$30% |

6 shots | 42.2 deg | $\u2248$29% |

46 shots | 45.7 deg | $\u2248$22% |

232 shots | 49.1 deg | $\u2248$17% |

2320 shots | 49.5 deg | $\u2248$16% |

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\u2009\mu m$), which may result in some deviations. The glaze icing conditions, which have been considered are listed in Table 3.

### 3.2 Ice Accretion Versus Free Stream Velocity ($V\u221e$) 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 ($y\u2245\u22120.01\u2009\u2009m$), 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 $2\u2009g/m3$, the film thickness starts to develop another peak located just after the horn (around $y\u22450.02\u2009\u2009m$), 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.

While temperature is a critical factor in influencing the shape of the final ice profile, the shape also depends on the freestream velocity ($V\u221e$) 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/m^{3} to 1 g/m^{3}, 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.

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 $\omega =0.5\u2009g/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\u221e=75\u2009m/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\u2009\mu 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.

### 3.4 Ice Accretion Versus Free Stream Temperature ($T\u221e$).

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\u221e$) 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\u221e=262K$ transforms to a thinner ice layer that is more spread out at the higher temperature $T\u221e=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.

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\u221e=262K$, when the size of the droplet increases from $100\u2009\mu m$ to $200\u2009\mu 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.

## 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/m^{3} to 1 g/m^{3}, 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)

*C*=_{d}drag coefficient

*c*=_{p}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

*k*=_{s}Equivalent sand-grain roughness, (m)

- $m\u02d9$ =
mass flux rate (kg m

^{–2}s^{–1})*n*=surface normal vector

*p*=order of grid convergence

- $P\xaf$ =
mean pressure (Pa)

- $P\u2032$ =
fluctuating pressure (Pa)

- $q\u02d9$ =
heat flux rate (J m

^{–2}s^{–1})*r*=mesh refinement ratio

*R*=surface roughness (mm)

*t*=time (s)

- $T\u0303$ =
film equilibrium temperature (°C)

- $T\xaf$ =
mean temperature (K)

- $T\u2032$ =
fluctuating temperature (K)

- $u\xaf$ =
film average velocity (ms

^{–1})=*V*Velocity vector (ms

^{–1})- $V\xaf$ =
mean velocity (ms

^{–1})- $V\u2032$ =
fluctuating velocity (ms

^{–1})

*α*=angle of attack

*α*=_{d}volume fraction of droplet

*ϵ*=relative error

*θ*=angle of ice horn (deg)

*μ*=dynamic viscosity (kg m

^{–1}s^{–1})*ν*=kinematic viscosity (m

^{2}s^{–1})*ρ*=density (kg m

^{–3})*σ*=droplet surface tension force (N m

^{–1})- $\tau $=
shear stress (Pa)

*ω*=liquid water content (g m

^{–3})