A mathematical model based on unstructured mesh for ice accretion

This paper proposes a three-dimensional model to simulate ice accretion on the surface of an aircraft. The model is developed using a prism grid cell, which enables iterations with an unstructured mesh. The model assumes that the impinged water droplets on the surface can form a thin water film, and ice accretion is caused by the flow and solidification of the film. The model is developed by analyzing the conservations of mass, momentum, and energy for the thin water film on the surface. The mass and energy conservation equations of the water film are discretized on the prism grid cell, while the momentum conservation equation, which is simplified using lubrication theory, is integrated along the surface normal direction to obtain an analytical solution for the film velocity. The thicknesses of the water film and ice layer are calculated by iterating the discretized mass and energy conservation equations with the analytical solution of the film velocity. The model is verified by comparing its results with published experimental and numerical data. The results show that the model provides accurate predictions of ice accretion under different conditions, even when using an unstructured grid.

NOMENCLATURE A = area of the cell surface Cw = heat capacity of water F = body force Hice = thickness of the ice layer Hw = thickness of the water film L f = latent heaṫ m = mass flow rate ⌢ n = unit vector for the surface normal direction p = pressure q = heat flow rate T = temperature u = velocity of the water film α = volume fraction β = water catch coefficient μ = dynamic viscosity ρ = density λ = heat conduction coefficient σ = surface tension τ = shear stress

I. INTRODUCTION
Ice accretion is a severe threat to flight safety 1,2 and has drawn significant attention from researchers around the world since the early 1950s. 1 Ever since, large efforts have contributed to numerical simulations of ice accretion, because of its advantages of lower monetary and time resources compared with experimental works. Messinger 3 first proposed a mathematical model for ice accretion in 1953. The model is based on the analysis of flow and solidification of unfrozen water in a two-dimensional (2D) rectangular control volume. In this model, all the unfrozen water in the control volume is assumed to flow downstream, and the ice layer thickness is calculated by iterating the mass and energy conservation equations. Based on this model, NASA developed the famous ice accretion simulation ARTICLE scitation.org/journal/adv code LEWICE. 4,5 Though Messinger's model was primarily developed for 2D simulations of ice accretion, many efforts have extended the model to three-dimensional (3D) simulations. Some have suggested using the original Messinger model to simulate ice profiles on different cutting-slices of a 3D surface and reconstruct the 3D ice profile from these slices. 6 Others proposed an extended method that uses the shear stress to distribute the unfrozen water on a 3D surface. 7,8 Thus, Messinger's model still has wide applicability when studying ice accretion. More sophisticated mathematical models were proposed in the late 1990s to simulate ice accretion. Bourgault et al. [9][10][11] developed a Shallow Water Icing Model (SWIM) to simulate ice accretion on a 3D surface. The SWIM approach assumes that not all the unfrozen water flows outside of the control volume and back downstream. Instead, the flow direction and mass flow rate of the water film are decided by the film velocity, which can be calculated using Newton's shear law. Therefore, they assume that the film velocity is decided solely by the shear stress on its surface. As there is a simple expression for the film velocity, the icing control volume can be either a prism or hexahedron, and the SWIM algorithm can be applied for both structured and unstructured grids. The famous commercial simulation code FENSAP is developed based on the SWIM architecture.
However, the flow of the thin film on a surface is naturally very complex; 12-14 different kinds of forces rather than just shear stress may contribute to its flow pattern. Myers et al. 15,16 believed that the velocity of unfrozen water is affected by more forces other than just the shear stress. They used a set of simplified Navier-Stokes expressions as the momentum equations of the water film to calculate its velocity. In the momentum equations, the effects of the shear stress, pressure gradient, gravity, and surface tension are considered. The temperature gradient along the direction of the film thickness direction is also considered in the model. The model was initially developed based on a Cartesian coordinate system and was only suitable for very simple Cartesian grids. 15 The model was later extended to an arbitrary surface using orthogonal curvilinear coordinates. 17 However, it is still difficult to apply the model to unstructured grids, because the surface grid requires orthogonalization to calculate the analytical solution of the water film flux. Chen et al. 18 improved Myers's model using nonorthogonal curvilinear coordinates, making it more suitable to simulate ice accretion on complex geometries. However, this model also needs to calculate the analytical solutions in the two main directions tangent to the surface, which adds difficulties in the application for unstructured grids.
Both Myers's 17 and Chen's 18 models require calculating the Lamé parameters introduced with curvilinear coordinates, which increases the calculation complexity and may introduce extra errors in the calculations. Villalpando et al. 19 extended Myers's model to use unstructured grids in the commercial code. However, they assumed that all the unfrozen water film in a grid cell runs downstream, which is the same as Messinger's model. In their paper, they suggested the gravitational force and air shear stress act on the water surface should be considered in future work.
In real flight conditions, ice can accrete on a complexgeometrical surface, which may be difficult to discretize using a structured grid. Therefore, it is necessary to develop an icing model that can be applied to unstructured grids and can also consider the effects of the shear stress, pressure gradient, gravity, and surface tension on the movement of the film as done in Myers's model. In this paper, a numerical model based on an unstructured mesh is proposed to simulate ice accretion. The model is developed based on analyses of the mass, momentum, and energy conservation equations of an unfrozen water film on a surface. The effects of the shear stress, pressure gradient, gravity, and surface tension are considered in the momentum conservation equation, and the temperature gradient along the film thickness is considered in the energy conservation equation. The accuracy of the current model is validated by comparing its numerical results with experimental data found in the literature.

II. MATHEMATICAL MODEL AND SOLUTION SCHEME
A. Simulation of the air/droplet flow field To simulate ice accretion, the mass of the impinged water droplets on the surface is needed, which can be obtained by simulating the air/droplet flow field. The Eulerian method 21 is used in this work to simulate the flow of the air and water droplets, which is not discussed in detail for conciseness. After the air/droplet flow field is simulated, the water catch coefficient is calculated as where ρw is the water density, α and v normal are the volume fraction and the normal velocity of the water droplets on an iced surface, respectively, and LWC and U∞ are the liquid water content and velocity of the water droplets for the incoming flow, respectively.

B. Mathematical model for the flow and solidification of the water thin film
The icing model is developed by analyzing the mass, momentum, and energy conservation equations for the water film using a prism grid cell, as shown in Fig. 1. The cell is divided into three layers: the bottom layer is the ice that contacts the solid surface; the middle layer is the water film that flows out of the control volume through the surrounding surface; and the top layer of the cell is the air/droplet flow, from which the water droplets impinge into the water film. The prism cell is chosen to ensure the model is applicable to unstructured grids.

Mass conservation
In the control volume shown in Fig. 1, the water droplets impact the surface of the water film with a mass ofṁimp ⋅ A sub , where A sub is the area of the bottom surface of the control volume (direct contact with the solid surface). The mass of the water that freezes into ice isṁice ⋅ A sub , and the evaporated mass of water isṁevp ⋅ A sub . The total mass of water that flows out of the control volume is n is the unit vector of the surface normal direction and ⃗ u is the velocity of the water film. The mass change of the water inside the control volume is expressed as ρw ∂H w ∂t A sub . Therefore, the mass conservation equation of the water film is

Momentum conservation
The flow of the shallow water film on the iced surface is considered to be an incompressible laminar flow. As the film is very thin, its momentum equation can be approximated using lubrication theory as where z is the direction normal to the solid surface, μ is the dynamic viscosity of water, p is the pressure of the water film, and ⃗ F is the mass force (gravity in stationary coordinates) that contains gravity, centrifugal forces, and Coriolis forces for rotating coordinates. Because the size of the impinged water droplet is very small (usually smaller than 50 μm), the effect of the impact of the water droplet on the momentum equation of the film can be neglected. 15 The momentum equation is given in differential form because the film velocity can be calculated analytically and does not need to be solved numerically which will be shown in Sec. II C.
The bottom of the thin water film (z = 0) directly contacts the solid surface, indicating the nonslip velocity boundary condition can be applied. The thin water film is driven by the shear stress at its top (z = Hw, where Hw is the thickness of water film), suggesting the Neuman boundary condition can be applied. Therefore, the velocity boundary condition for Eq. (3) can be written as where ⃗ τ is the shear stress exerted on the top of the water film by the air flow. As the water film is relatively thin, the component of the water film velocity at the normal direction of the surface can be neglected, Both the air pressure and the surface tension contribute to the pressure at the interface between the water film and the air. Therefore, the boundary condition for the pressure at the top of the water film is where pσ is the capillary pressure induced by the surface tension, which can be calculated using the Yong-Laplace equation. The thinness of the water film suggests that the pressure does not change along the z-direction.

Energy conservation
By analyzing the energy conservation of the water film in the control volume, we obtain where Cw is the heat capacity of water,Tw = 1 Twdz is the mean temperature in the water film,qevap =ṁevp ⋅ Le is the energy removed via evaporated water where Le is the latent heat of evaporation,q cond is the heat transferred from the water film to the substrate (zero in the current icing simulation),qconv = h( Tw| z=H w − Ta) is the heat transferred from the water film to the outside airflow via convection,qimp =ṁimp ⋅ is the heat that enters the film with the impinged water, andqice = ρiceL f ∂H ice ∂t is the heat generated by ice accretion.
To calculate the energy conservation equation from Eq. (7), we need to know the temperature inside the water film. As the water film is very thin, we assume that the temperature changes linearly along the z-direction, which gives The bottom surface of the water film is the interface between the water and the ice; therefore, its temperature should be at freezing which is The top surface of the water film contacts the air/droplet flow field, and the Robin boundary condition can be applied where λw is the heat conduction coefficient. Integrating Eq. (8) with the boundary conditions of Eqs. (9) and (10), we can express the temperature inside the water film as where n is the iteration number and Ai is the net mass of water that flows out of the control volume, which depends on the mean water film velocity on the surface around each cell. The mean velocity of the water film can be calculated by integrating the momentum conservation equation of Eq. (3) along z with the boundary condition of Eq. (4), which gives As can be seen, the film velocity does not change linearly along the normal direction of the surface as assumed by Bourgault,9 because the effects of the pressure and body force are considered. According to Eq. (5), the component of the water film along the z direction can be neglected, To calculate the net mass of water that flows out of the control volume, we need to know the mean velocity of the water film on the surface. The mean velocity along the film thickness direction at each cell center is calculated as The mean water film velocity on the surface between two cells is calculated using a second-order upwind scheme as In Eq. (15), ∇p is the gradient of the pressure in the center of a grid cell, which can be calculated using Green-Gauss theorem as where p f is the pressure on the surface around the cell and ⃗ A f is the surface area in the normal direction that points away from the cell. The pressure at the surface center can be calculated based on the pressure at the nodes of the surface as pn, node is the value of pressure on the node of the cell which contains both the air pressure and capillary pressure caused by the surface tension.
Calculating the film velocity in the form of vectors Eq. (15) with the boundary conditions Eq. (14) can ensure that the current model does not need to find the two principal axis that are tangent to the surface as done by Myers et al. 17 and Chen et al. 18 This allows the model to be calculated in an unstructured grid. The calculation schemes of the film velocity on the interface Eq. (16) and pressure gradient Eqs. (17) and (18) also make the model be able to applied for unstructured grid.
The energy conservation equation in Eq. (7) is discretized as According to Eqs. (15) and (11), both the velocity of the water film ⃗ u and its temperature depend on its thickness, so (19) can be written as f (Hw) and g(Hw), respectively. Therefore, the discretized mass conservation equation in Eq. (12) and the energy conservation equation in Eq. (19) can be rewritten using the Crank-Nickson scheme as ρw whereTw is calculated using Eq. (11). By solving Eqs. (20) and (21), we can obtain the thicknesses of the water film and ice layer. As shown in Eqs. (20) and (21), the convective terms are calculated on both time step n and n + 1 to achieve a stable and second-order accuracy solution. When it is rime ice, all the impinged water will freeze into ice and the ice thickness is H n+1 ice = H n ice +ṁ imp ρ ice Δt. The icing model proposed in this manuscript can be applied for meshes formatted in unstructured scheme. It can be applied for tetrahedral grids with prism boundary layer or hexahedral mesh. However, it cannot be applied for fully tetrahedral mesh.

D. Implementation into commercial code
The model proposed can be implemented into commercial Computational Fluid Dynamic (CFD) code, such as FLUENT. 20 One can use the User Defined Function (UDF) function in FLUENT to find all the boundary cells on the wall, and then calculate Eqs. (20) and (21) in these boundary cells. The detailed steps to simulate ice accretion in FLUENT with the proposed icing model are shown in Fig. 2. The first step is to use the Fluent-Solver to simulate the dry air flow. The second step is to simulate the water droplets using the Eulerian method 21  time step Δt in Eqs. (20) and (21). The choice of Δt needs to satisfy the CFL (Courant-Friedrichs-Lewy) condition. The fourth step is to update the geometry and mesh according to the increased ice mass, which can be accomplished using the dynamic-mesh function. After these four steps are completed, if the time has reached the set-time, the calculation is ended, otherwise the new geometry and mesh are used, and the process is repeated from the first step.

A. Validation of ice accretion on an NACA0012 airfoil
The proposed model requires verification, which is achieved using ice accretion on an NACA0012 airfoil. The chord length of the airfoil is C = 0.5334 m, and three different icing conditions from the literature are used. The detailed conditions are shown in Table I The main purpose of this paper is to provide an icing model that is suitable for unstructured grids. Therefore, an unstructured grid with a prism boundary layer mesh is adopted for validation (see Sec. A in supplementary material). Mesh independence test was done by using four meshes with different sizes. The results can be found in Sec. A in the supplementary material, and the mesh of 781 thousands nodes is used for analysis in this section.
An initial time step of Δt = 0.05 s which can vary with the film velocity to meet the CFL condition is used in the current study. Δtime = 105 s is used for Case 1 and 3, while Δtime = 90 s is used for case 2. Therefore, the dynamic mesh function is called for four times for all these three cases. A minimum length scale of 1 mm is used for the deformation zone in the dynamic mesh function. Figure 3 depicts the calculated water collection efficiency β using the Eulerian method for Case 1 and its comparison with the data from Ref. 21, which are calculated using the Lagrangian method. In Fig. 3, S is the arc length of the airfoil profile, which begins form the stagnation point (S = 0) with S < 0 on the pressure side and S > 0 on the suction side. As observed, β reaches its maximum of near 1 when S/C = 0, because the moving direction of the water droplets does not readily change at this position. It is also seen in Fig. 3 that the impact limit on the side for S/C < 0 is larger than that on the side for S/C > 0, which is due to the attack angle being nonzero. Figure 3 also shows that the calculated β matches well with the data in Ref. 21.
The heat transfer coefficient (HTC) is an important parameter in the simulation of ice accretion. Thus, the predicted ice profile can be quite different from the actual scenario if the predicted HTC is inaccurate. This paper adopts the method suggested in Ref. 23 to calculate the HTC, which uses a realizable k-ε turbulence model together with fixed temperature boundary conditions. The roughness of the surface is calculated by using the model provided by NASA. 4 Figure 4 shows the calculated HTC for Case 1 and its comparison with the data from Ref. 21, which is calculated using the LEWICE algorithm. These results match well for most of the airfoil, except that the maximum value of the HTC is slightly smaller than the LEWICE code results around S/C = 0.1. Figure 5 shows the calculated thickness and velocity vectors of the unfrozen water film, and Fig. 6 shows the shear stress and its direction exerted by the air on the surface, both for Case 1. As seen, the magnitude and direction of the film velocity are coincident with those of the shear stress, which implies that the shear stress is the dominant factor impacting the film velocity in the condition of Case 1. However, when the film thickness is large, the pressure gradient can also affect the flow of the film (see Sec. B in supplementary material). And when the surface is rotating, the centrifugal force can draw the film flow in the span direction. Therefore, it is necessary  to include the pressure gradient and body force term in Eq. (15). Figure 5 also shows that the water film is thickest at the stagnation point of the airfoil because the amount of impinged water is large (see Fig. 3) while the heat transfer coefficient is small (see Fig. 4). This causes a significant amount of water to remain unfrozen. The film velocity at the stagnation point is also relatively small (see Fig. 5), which means little unfrozen water will flow downstream with most remaining at the stagnation point. Figure 7 shows the change of thickness of water film and ice layer with time on the airfoil for Case 1. As can be seen, the film thickness at the stagnation point keeps at the maximum value for different time steps, while it increases with time goes by. The film covered area also increases with time. This is becauseṁimp,ṁice anḋ mevp can reach a balance with the net mass flux flow out of the cell at the stationary point. However, water will be driven into the cells in the downstream from the stationary point. Therefore, the thickness and the covered area of water film in the downstream increase with time. Figure 7 also shows that the ice-covered area at the up surface (S/C > 0) matches that covered by the water film, which implies it is glaze ice at the up surface. However, at the down surface (S/C < 0) the ice-covered area is same as the impact-area of droplets (see in Fig. 3) and is larger than that covered by the water film. Therefore, the results show that rime ice and glaze ice occur at the same time on the airfoil.
The details of the ice profile at t = 420 s are seen from the blue line with solid circles in Fig. 8. Because the refined boundary mesh is used near the wall, the differences in ice shape on cutting slices of different locations are unrecognizable. The ice profile shows an obvious ice horn feature, and the existence of the attack angle causes the upper horn to be larger than the lower horn.

B. Validation of ice accretion on a cylinder
The experimental data for ice accretion on a cylinder provided by Koss 24 is adopted here to further validate the proposed model and its solution scheme. Three cases are chosen with the icing conditions  , and the mesh of 572 thousands nodes is used for analysis in this section. Figure 11 shows the calculated thickness of the unfrozen water film on the cylinder surface of Case 1. The largest water film thickness occurs at the stationary point of the cylinder because there is a larger amount of impinged water there (see Fig. 12). As the water droplet inertia is large, it is difficult to change its moving direction, and most of the water droplets impinge onto the leading edge of the cylinder, which causes β to be large there. Figure 13 compares the calculated ice profile with the experimental data from Case 1. The calculated ice profile matches well with the experimental data, except that the calculated ice thickness is slightly larger than the experimental case near the stationary point of the cylinder, which may be due to deviations between the calculated HTC and the actual value. Figure 14 compares the calculated ice profiles with the experimental data for Cases 2 and 3, respectively. The calculated ice profiles match well for both cases, which means the current model and its solution scheme have good precision. As the incoming velocity for Case 3 is smaller than that for Case 2, the ice thickness on the cylinder for Case 3 is smaller than in Case 2.
Though the profiles of the current two geometries selected for validation are simple, simulations with 3D unstructured grid have been done and the results can show the precision and capacity of   the current model. Simulation of ice accretion on a more complex geometry can be find in Sec. C in the supplementary material, which provides a future validation of the current model.

IV. CONCLUSIONS
This paper provides a 3D numerical model to simulate ice accretion, which can be applied to unstructured grids, and can also consider the effects of the shear stress, pressure gradient, gravity, and surface tension on the movement of the film. The model is developed based on the analyses of the mass, momentum, and energy conservation equations for an unfrozen water film with a prism grid cell. In the momentum conservation equations, the effects of the shear stress, pressure gradient, gravity, and surface tension are considered. The analytical solution for the water film velocity is obtained by integrating the momentum conservation equation. The mass and energy conservation equations for the thin water film are discretized using the prism grid cell. The thickness of the water film and ice layer are calculated by iterating the discretized mass and energy conservation equations with a combined Crank-Nickson scheme. To verify the current model, different conditions for ice accretion are tested on a NACA0012 airfoil and a cylinder, and the results are compared with data from the literature. The comparisons show that the model can be used with an unstructured grid to provide accurate predictions of ice accretion under different conditions.

SUPPLEMENTARY MATERIAL
See Sec. A in supplementary material for the complete mesh independence analysis.
See Sec. B in supplementary material for the complete analysis of the pressure gradient on the film flow.
See Sec. C in supplementary material for the complete simulation on a complex geometry.