An Unsteady Model for Aircraft Icing Based on Tightly-Coupled Method and Phase-Field Method

Abstract

An unsteady tightly-coupled icing model is established in this paper to solve the numerical simulation problem of unsteady aircraft icing. The multi-media fluid of air and droplets is regarded as a single medium fluid with variable material properties. Taking the droplet concentration as the phase parameter and the droplet resistance coefficient as the interphase force, the mass concentration distribution of the droplet is obtained by solving the Cahn–Hilliard equation. Fick’s law is introduced to improve the Cahn–Hilliard equation to predict the droplet shadow zone. On this basis, the procedure of the unsteady numerical simulation method for aircraft icing is established, including grid generation, the dual-time-step method to realize the unsteady calculation of the air and droplet tightly-coupled mixed flow field, and the improved shallow water icing model. Finally, through the comparative analysis of numerical examples, the effectiveness of the new model in predicting the droplet impact characteristics and the droplet shadow zone are verified. Compared with other icing models, the ice shapes predicted by the unsteady tightly-coupled model were found to be the most consistent with the experiments. In the icing comparison conditions in this manuscript, the prediction accuracy of the ice thickness at the stagnation point of the leading edge was up to 35% higher than that of LEWICE.

1. Introduction

The aircraft icing problem has been one of the most important problems in flight safety [1]. This problem can destroy the aerodynamic shape [2] and increase the resistance of aircraft [3], which is an important hidden danger leading to flight safety accidents. Due to the low cost and high efficiency [4], the ice shape prediction model and its numerical simulation method have become important methods to research the aircraft icing problem.

At present, a number of mature icing numerical simulation codes have been developed and widely used in industry, such as LEWICE developed by NASA [5,6], FENSAP-ICE developed by NTI (Newmerical Technologies International) [7,8], TRAJICED developed by DRA [9], IGLOO2D and IGLOO3D developed by ONERA [10] and so on. Although the numerical simulation methods of these software programs are different, the main numerical process can be divided into four calculation modules, including grid generation, airflow field calculation, droplet trajectory calculation and ice accretion calculation. The numerical simulation methods of droplet trajectory are mainly divided into Lagrangian method and Eulerian method.

The Lagrangian method takes a single droplet as the research object. By tracking its motion trajectory in the flow field, it can judge whether it collides with the surface of the airfoil, and then obtains the droplet collection coefficient. The Lagrangian method is often used to calculate the simple shape of the collection efficiency and is still widely used in LEWICE, ONERA, etc.

The Eulerian method, which has been applied to FENSAP-ICE, takes the continuous distribution of droplets in space as the research object. By solving the governing equation of the water phase, the mass and velocity distribution of droplets in the airflow field are calculated, and then the droplet collection coefficient of the icing surface is obtained. Both the Lagrangian method and Eulerian method ignore the influence of droplets on the airflow field; therefore, they can be considered as loosely-coupled models. The tightly-coupled model proposed in this paper regards the multi-media fluid composed of air and water droplets as a single fluid with variable material properties and focuses on the interaction between air and droplets.

The phase field method, also known as the diffuse interface method [11], has now become a powerful method to simulate multi-media fluid flow problems. It is commonly used in droplet coalescence, fragmentation, rise and change in shear flow, phase separation, contact line dynamics, surfactant adsorption interface dynamics and thermocapillary effects. Compared with other interface capturing methods, the phase field method has the following two advantages: first, it is easy to capture the interface using the numerical method even if the topology changes; secondly, the phase field method can be obtained with the energy-based variational method.

The results show that the phase field method is consistent with the law of energy dissipation, which makes it possible to design an efficient and energy stable numerical scheme. The phase-field method treats multi-phase fluids as one fluid with variable material properties [12]. Specifically, the phase-field method utilizes order parameters, such as the concentration or volume fraction of the fluids, to represent different fluid phases. The order parameter is typically evolved by the Cahn–Hilliard (C-H) equation, the Allen–Cahn equation or other types of dynamic equations to capture the motion of the fluid interface.

In the previous numerical simulation of aircraft icing, the small supercooled droplets with diameters ranging from 10 to 50 μm were regarded as rigid spheres, assuming that there is no deformation in the process of motion and always remains spherical, and their resistance coefficient can be obtained by the classical empirical formula of spherical resistance coefficient [13]. For supercooled large droplets with a diameter of more than 50 μm, Tan et al. [14] observed the process of large droplets impacting the airfoil and found that large droplets will deform or even eventually break when they are close to the airfoil. By adding the SLD deformation effect, they obtained a droplet collection coefficient closer to the experimental results.

Aircraft icing is an unsteady process. At every tiny time step, the growth of the ice layer at the leading edge of the airfoil will change its aerodynamic characteristics, which will affect the air flow field and droplet trajectory and eventually lead to the change of the ice shape [15]. In the current icing calculation methods, this is regarded as a steady-state or quasi steady-state problem and solved by the single-time-step method or multi-time-step method.

The multi-time-step method is introduced into the LEWICE, which can be divided into two categories: the first method [16,17,18] divides the total icing time into several time steps, and calculates the airflow field, droplet trajectory and ice growth in each time step. The specific process is shown in Figure 1a. However, because each time step is a complete and independent freezing process, the calculation consumes a great deal of time. In the second method [19,20], the airflow and droplet trajectory are only calculated at the initial time.

The ice growth process is divided into multi-time steps, and each time step is a complete and independent single-time-step ice growth calculation, which is shown in Figure 1b. While the second method uses an unsteady method to improve the ice accretion model, this method is only suitable for low to moderate ice thicknesses or simple flow in the time step.

Reid et al. [7] developed an unsteady method for numerical simulation of in-flight electrothermal anti-icing or de-icing, using the conjugate heat transfer technique. In this numerical method, both the flow field and the droplet trajectory are solved in the steady state, and only the calculation method of icing and melting is unsteady. Zhang [21] and Ansell [22] successively developed a numerical simulation method of unsteady air flow field to characterize the sources of unsteadiness in iced airfoil aerodynamics.

On the basis of the Myers model [23], Gori [24] established a new model for ice accretion under rime and glaze conditions, which includes an unsteady description of the heat diffusion problem within the ice layer. Based on the numerical solution of the unsteady Stefan problem, Liu [25] established a three-dimensional aircraft icing model, which regards the heat conduction in the ice layer as an unsteady process. Xiao [26] used wall-modeled large-eddy simulation to investigate the unsteady flow characteristics for a 30P30N three-element iced airfoil.

The above unsteady studies of aircraft icing numerical simulation have thus far focused on the unsteady research of airflow field and the unsteady heat transfer of ice accretion model. However, these methods ignore that the aircraft icing is an unsteady process. Zhu [27] proposed an adaptive Cartesian method combined with the ghost cell method for the aircraft unsteady icing problem. However, this numerical method calculates the airflow field by solving Eulerian equations and then uses the Lagrangian method to obtain the droplet trajectories, and the solution of the droplet trajectories is still steady.

In view of this situation, in this paper, we propose a numerical simulation method of the unsteady tightly-coupled icing process. In this method, the multiphase flow of air and droplets are regarded as a single fluid with variable material properties. The mass concentration of droplets is selected as the phase parameter, and Fick’s law is introduced to improve the Cahn–Hilliard equation to predict the droplet shadow zone. At the same time, the dual-time-step method is introduced to calculate the unsteady multiphase flow. The shallow water icing model is also improved to be suitable for unsteady calculation. The flow chart of unsteady icing calculation is shown in Figure 2.

The rest of the paper is organized as follows. In Section 2, the unsteady numerical method of aircraft icing is established, including grid generation, the calculation of air and water droplet tightly-coupled mixed flow field and the improved shallow water icing model. In Section 3, the correctness of the droplet collection coefficient calculated by the tight coupling model is verified. The accuracy of the new model in predicting the droplet shadow zone is presented in Section 4. In Section 5, the ice shapes are provided by the unsteady icing model under a series of icing conditions and compared with the experimental ice shapes and those provided by LEWICE. Finally, the main conclusions of the present work are provided in Section 6.

2. Materials and Methods

2.1. Grid Generation

In order to improve the computational efficiency, the grid partition method is introduced to divide the whole computational domain into several sub regions for grid generation, and the two-dimensional airfoil structure mesh is generated by solving the two-dimensional Poisson equations [28].

For two-dimensional airfoil icing, the icing mainly occurs at the leading edge of the airfoil, and the shape of other parts of the airfoil does not change; therefore, the corresponding mesh does not need to be updated and reconstructed. Thus, as shown in Figure 3, the whole calculation area is divided into the sectorial domain and non sectorial domain. The sectorial domain is the area where icing occurs, and its grid needs to be adjusted continuously with the advance of icing time, while the non sectorial domain only needs to be generated once at the initial time. The computational grid is shown in Figure 4.

In addition, according to our research, when the ice shape change in each physical time step is small enough, the mixed flow variables of air and droplets in the previous time step can be directly transferred to the new mesh without any correction by mesh velocity.

2.2. Tightly-Coupled Mixed Flow Field of Air and Droplets

In order to research the unsteady icing problem, the phase-field method is used to solve the tightly-coupled mixed flow field of air and droplets. In the control volume $\Omega$, the mass concentration of supercooled droplet is defined as the phase variable:

$$\varphi ={\varphi }_w=\frac{{\overline{\rho }}_w}{\rho },{\varphi }_a=1-{\varphi }_w.$$

(1)

Then, the average density of the gas–liquid multiphase flow is defined by harmonic interpolation:

$$\frac{1}{\rho }=\frac{\varphi }{{\rho }_w}+\frac{1-\varphi }{{\rho }_a}.$$

(2)

The viscosity coefficient $\eta$ and thermal conductivity k of multiphase flow are also calculated by harmonic interpolation:

$$\frac{1}{\eta }=\frac{\varphi }{{\eta }_w}+\frac{1-\varphi }{{\eta }_a},\frac{1}{k}=\frac{\varphi }{k_w}+\frac{1-\varphi }{k_a}.$$

(3)

where ${\eta }_w$ and ${\eta }_a$ are, respectively, the viscosity coefficients of water and air, $k_w$ and $k_a$ are, respectively, the thermal conductivities of water and air.

For any calculation unit, assuming that the two fluids move at different velocities $u_i(x,t)$, the mass conservation equation of each fluid is expressed as follows:

$$\frac{\partial {\overline{\rho }}_i}{\partial t}+\nabla ·m_i=0.$$

(4)

where $m_i$ is the mass flow rate in the calculation unit, including the mass flow caused by convection and the mass flow caused by diffusion, that is $m_i={\overline{\rho }}_i{\mathbf{U}}_i+{\rho }_i{d}_i$, and $d_i$ is the diffusion flow rate of fluid i. Two kinds of fluids are characterized by i = 0 and 1, respectively. When i = 0, it means fluid 0, and when i = 1, it means fluid 1.

According to Equation (4), the mass conservation equation of multiphase flow can be obtained by:

$$\frac{\partial ({\overline{\rho }}_a+{\overline{\rho }}_w)}{\partial t}+\nabla ·({\overline{\rho }}_a{\mathbf{U}}_a+{\overline{\rho }}_w{\mathbf{U}}_w+{\rho }_a{d}_a+{\rho }_w{d}_w)=0.$$

(5)

where ${\mathbf{U}}_a$ is the velocity of air , and ${\mathbf{U}}_w$ is the velocity of droplet. The average density $\rho$ and the average velocity $\mathbf{U}$ in the calculation unit are defined as follows:

$$\rho ={\overline{\rho }}_a+{\overline{\rho }}_w,\rho \mathbf{U}={\overline{\rho }}_a{\mathbf{U}}_a+{\overline{\rho }}_w{\mathbf{U}}_w.$$

(6)

Taking Equation (6) into Equation (5), the mass conservation equation of multiphase flow can be obtained as follows:

$$\frac{\partial \rho }{\partial t}+\nabla ·(\rho \mathbf{U}+{\rho }_a{d}_a+{\rho }_w{d}_w)=0.$$

(7)

The diffusion mass flow rates of the two fluids should be equal in absolute value and opposite in direction, ${\rho }_a{d}_a+{\rho }_w{d}_w=0$.

To sum up, the calculation formula of the mass equation of multiphase flow can be described as:

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathbf{U}\right)=0.$$

(8)

In the following derivation, the multiphase flow is regarded as a single fluid moving at velocity $\mathbf{U}$.

The Cahn–Hilliard (C-H) equation and free energy density formula are as follows:

$$\frac{\mathrm{D}\varphi }{\mathrm{D}t}=M_c\Delta \mu .$$

(9)

where $M_c$ is the nonnegative permeability of the phase field. The free energy density $\mu$ can be obtained by:

$$\Delta \mu =\frac{\varphi (u-u_w)}{M_c}.$$

(10)

For the simple airfoil icing problem, this effect can be ignored; however, for the complex aircraft icing problem, the droplet impact characteristics and droplet collection coefficient of the object surface near the trailing edge of the airfoil may be affected, resulting in numerical inaccuracy. Considering the diffusion effect caused by concentration differences, the range of the droplet shadow zone can be predicted more accurately. In order to take the diffusion effect of droplets into account, a macroscopic motion law is introduced to describe the phenomenon of matter diffusion, namely Fick’s law. According to the definition of Fick’s law, the phase equation of a water droplet can be written as follows:

$$\frac{\mathrm{D}\varphi }{Dt}=D\Delta \varphi .$$

(11)

where D is the diffusion coefficient.

In conclusion, the phase equation can be written as follows:

$$\frac{\mathrm{D}\varphi }{Dt}=M_c\Delta \mu +D\Delta \varphi .$$

(12)

The force exerted by air on droplets is mainly composed of pressure difference resistance and friction resistance, and the portion of pressure versus resistance forces is different for different shaped objects. According to the literature, whether it is a spherical body or a disk-shaped object, the proportion of pressure difference resistance to the total resistance is more than 90%. Therefore, when we study the drag coefficient of water droplets, we typically only consider the pressure drop resistance and ignore the friction resistance. According to Newton’s third law, the force of water droplets on air is equal to that of air on droplets and directed to contrary parts as shown in Figure 5.

According to Newton’s third law of motion (force and acceleration), the force formula of a droplet on air can be written in the form [16]:

$${\mathbf{F}}_{wa}=-{\mathbf{F}}_{aw}={\rho }_w{V}_{w}K({\mathbf{U}}_w-{\mathbf{U}}_a).$$

(13)

where, $V_w$ is the volume of the droplet, and K is the air–droplet interaction factor, which is specified in the form:

$$K=\frac{18{\eta }_{a}f}{{\rho }_w{d}^2}.$$

(14)

where the factor f is the drag function and defined differently for different exchange coefficient models. In this paper, according to the Schiller–Naumann model, the drag function is given as:

$$f=\frac{C_{D}Re}{24}.$$

(15)

When the diameter of the droplet is small, the aerodynamic force of the droplet is smaller than the surface tension of the droplet. Therefore, it is generally considered that the droplet can remain spherical in the process of movement without deformation, and its drag coefficient can be obtained by the formula of the spherical resistance coefficient. However, when the diameter of the water droplet is large, the aerodynamic force of the droplet will exceed its surface tension and become the dominant factor in the process of movement, which makes it difficult to keep the droplet spherical and deform in the process of movement, , eventually leading to the change of aerodynamic force and motion state. Therefore, the deformation of large droplets can no longer be ignored.

Honsek [29] regarded the deformation of droplets as a spring system, and considered that the shape of droplets in the deformation process is between the rigid sphere and the disk, and the resistance coefficient between them can be obtained as:

$$C_D=f_e{C}_{D,Sphere}+(1-f_e)C_{D,Disk}.$$

(16)

where, $C_{D,Sphere}$ is the drag coefficient of a spherical droplet and $C_{D,Disk}$ is the droplet drag coefficient of the equivalent disk. The eccentric action coefficient $f_e$ is written in the form:

$$f_e={(1+0.007\sqrt{We})}^{-6}.$$

(17)

$We$ is the dimensionless relative Weber number, which is used to characterize the ratio of the surface tension to inertial force of water droplets and is given as

$$We=\frac{{\rho }_a\mathrm{MVD}{({\mathbf{U}}_a-{\mathbf{U}}_w)}^2}{\sigma }.$$

(18)

where ${\rho }_a$ is the density of air, $\mathrm{MVD}$ is the droplet diameter, and $\sigma$ is the surface tension coefficient.

The calculation formulas of the drag coefficient of a spherical water drop and that of a water drop equivalent to a disc are given as a wide range function of the Reynolds number [30,31]:

$$C_{D,Sphere}=\frac{24}{Re}(1+0.15R{e}^{0.687}).$$

(19)

$Re$ is the relative Reynolds number of water droplets, and the calculation formula can be given as follows:

$$Re=\frac{{\rho }_a|{\mathbf{U}}_a-{\mathbf{U}}_w|d}{{\mu }_a}.$$

(20)

By inserting the force into the original calculation model of air flow and droplet [32], the improved continuity equation and energy equation can be obtained as follows:

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }{\mathbf{W}}^{n+1}d\Omega +\frac{\partial }{\partial \tau }{\int }_{\Omega }{\mathbf{W}}^{n+1}d\Omega +\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\partial \Omega }({\mathbf{F}}_c-{\mathbf{F}}_v)dS=\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }\mathbf{Q}d\Omega ,$$

(21)

$$\mathbf{W}=\left[\begin{array}{c}\rho \\ \varphi \\ {\rho }_a{u}_a\\ {\rho }_a{v}_a\\ {\rho }_w\alpha u_w\\ {\rho }_w\alpha v_w\\ \rho E\end{array}\right],{\mathbf{F}}_c=\left[\begin{array}{c}\rho V\\ \varphi \mathbf{U}\\ {\rho }_a{u}_a{\mathbf{U}}_a+n_{x}p\\ {\rho }_a{v}_a{\mathbf{U}}_a+n_{y}p\\ {\rho }_w\alpha u_w{\mathbf{U}}_w\\ {\rho }_w\alpha v_w{\mathbf{U}}_w\\ \rho H\mathbf{U}\end{array}\right],{\mathbf{F}}_v=\left[\begin{array}{c}0\\ 0\\ n_x{\tau }_{xx}+n_y{\tau }_{xy}\\ n_x{\tau }_{yx}+n_y{\tau }_{yy}\\ 0\\ 0\\ n_x{\Theta }_x+n_y{\Theta }_y+M_c\mu \Delta \mu \end{array}\right],$$

$$\mathbf{Q}=\left[\begin{array}{c}0\\ 0\\ {\rho }_w\alpha K(u_w-u_a)+{\rho }_a{g}_x\\ {\rho }_w\alpha K(v_w-v_a)+{\rho }_a{g}_y\\ {\rho }_w\alpha K(u_a-u_w)+{\rho }_w\alpha g_x\\ {\rho }_w\alpha K(v_a-v_w)+{\rho }_w\alpha g_y\\ 0\end{array}\right].$$

In which t is the physical time step, and $\tau$ is the virtual time step. The intermediate variables of viscous stress ${\Theta }_x$ and ${\Theta }_y$ are expressed as follows:

$$\begin{array}{cc}\hfill {\Theta }_x& =u{\tau }_{xx}+v{\tau }_{xy}+w{\tau }_{xz}+k\frac{\partial T}{\partial x},\hfill \\ \hfill {\Theta }_y& =u{\tau }_{yx}+v{\tau }_{yy}+w{\tau }_{yz}+k\frac{\partial T}{\partial y}.\hfill \end{array}$$

(22)

At present, the most widely used unsteady calculation methods in computational fluid dynamics can be divided into two categories: the physical time iteration method and dual-time-step method [33]. When the physical time step iterative method is used to deal with unsteady flow problems, due to the limitation of stability conditions, the time step needs to be smaller, which results in a longer calculation time.

Faced with this issue, the dual-time-step method proposed by Jameson for unsteady calculation has incomparable advantages. In the dual-time-step method, the whole calculation process is divided into several physical time steps, and the virtual time step is adopted in each physical time step to convert the unsteady problem into the steady problem. At the same time, the local time step and residual smoothing method can be used to accelerate the convergence and improve the calculation efficiency.

The whole calculation process is divided into several physical time steps. Then, on the $n+1$ physical time step, the fully implicit form of the governing equation is as follows:

$$\frac{\mathrm{d}}{\mathrm{d}t}\left({\Omega }^{n+1}{\mathbf{W}}^{n+1}\right)+\mathbf{R}\left({\mathbf{W}}^{n+1}\right)=0.$$

(23)

The second order backward difference scheme is used for the time term of Equation (23), and its expression is written as follows:

$$\frac{\mathrm{d}}{\mathrm{d}t}\left({\Omega }^{n+1}{\mathbf{W}}^{n+1}\right)=\frac{3{\Omega }^{n+1}{\mathbf{W}}^{n+1}-4{\Omega }^n{\mathbf{W}}^n+{\Omega }^{n-1}{\mathbf{W}}^{n-1}}{2\Delta t}.$$

(24)

Introducing Equation (24) into (23), we can obtain:

$$\frac{3{\Omega }^{n+1}{\mathbf{W}}^{n+1}-4{\Omega }^n{\mathbf{W}}^n+{\Omega }^{n-1}{\mathbf{W}}^{n-1}}{2\Delta t}+\mathbf{R}\left({\mathbf{W}}^{n+1}\right)=0.$$

(25)

Here, we define a new residual value ${\mathbf{R}}^{*}\left({\mathbf{W}}^{*}\right)$ as follows:

$${\mathbf{R}}^{*}\left({\mathbf{W}}^{*}\right)=\mathbf{R}\left({\mathbf{W}}^{*}\right)+\frac{3{\Omega }^{n+1}{\mathbf{W}}^{*}-4{\Omega }^n{\mathbf{W}}^n+{\Omega }^{n-1}{\mathbf{W}}^{n-1}}{2\Delta t}.$$

(26)

where ${\mathbf{W}}^{*}$ is the approximation of ${\mathbf{W}}^{n+1}$. By introducing the virtual time step $\tau$, the implicit unsteady equation is transformed into the steady form of the virtual time:

$$\frac{\partial }{\partial \tau }\left({\Omega }^{n+1}{\mathbf{W}}^{*}\right)+\mathbf{R}\left({\mathbf{W}}^{*}\right)=0.$$

(27)

In the iterative process of virtual time step, the fourth-order Runge–Kutta method is used to solve the problem explicitly. The virtual time step is limited by the stability condition, and its expression is as follows

$$\Delta \tau =min[\frac{CFL·\Omega }{({\Lambda }_c^I+{\Lambda }_c^J+{\Lambda }_c^K)+C({\Lambda }_v^I+{\Lambda }_v^J+{\Lambda }_v^K)},\frac{2}{3}\Delta t].$$

(28)

The finite volume method based on the cell centered form is applied for the discretization in space. An improved Jameson center scheme is used to calculate the inviscid fluxes. In order to suppress the numerical oscillation, the artificial dissipation term is also introduced [34]. The S-A model is chosen for the turbulent model.

2.3. Ice Accretion Model

The shallow water icing model (SWIM) [35] is an improvement over the Messinger model. Its greatest advantage is the unsteady nature of the SWIM and the way it models the water film for more complex configurations compared with a simple 2D airfoil. Therefore, the SWIM is introduced into the icing calculation and needs to be improved to be used in unsteady model.

As shown in Figure 6, in the shallow water icing model, the water film flow is simplified as a Coutte flow under the action of air shear stress, and the velocity of water film can be approximately described in the form of a linear distribution:

$$\mathbf{U}(s,y)=\frac{y}{{\mu }_w}{\tau }_a\left(s\right).$$

(29)

where s is the arc direction of the surface, and y is the normal direction of the surface. Then, the average flow rate of the water film can be defined as:

$$\overline{\mathbf{U}}\left(s\right)=\frac{1}{h_w}{\int }_{h_i}^{h_w+h_i}udy=\frac{h_w}{2{\mu }_w}{\tau }_a\left(s\right).$$

(30)

In which $h_w$ is the thickness of the water film, $h_i$ is the thickness of the ice layer, and ${\tau }_a\left(s\right)$ is the shear stress of the air flow on the water film. The calculation formula of the mass conservation equation in the control volume is as follows:

$$\frac{\mathrm{d}}{\mathrm{d}t}\left({\rho }_w{h}_w\right)+\frac{\mathrm{d}}{\mathrm{d}s}{\int }_{h_i}^{h_w+h_i}{\rho }_w\mathbf{U}dy={\dot{m}}_{imp}-{\dot{m}}_{evp}-{\dot{m}}_{ice}.$$

(31)

In which $\frac{\mathrm{d}}{\mathrm{d}s}{\int }_{h_i}^{h_w+h_i}{\rho }_w\mathbf{U}dy$ is the mass balance of inflow and outflow in the control volume, and $\frac{d}{dt}\left({\rho }_w{h}_w\right)$ is the represents the increase of water film height with physical time. For the unsteady model used in this paper, it can be described as:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left({\rho }_w{h}_w\right)& ={\rho }_w\frac{h_w^{n+1}-h_w^n}{\Delta t},\hfill \\ \hfill \frac{\mathrm{d}}{\mathrm{d}s}{\int }_{h_i}^{h_w+h_i}{\rho }_w\mathbf{U}dy& ={\dot{m}}_{out}-{\dot{m}}_{in}.\hfill \end{array}$$

(32)

where $h_w^{n+1}$ is the water film height of the current physical time step. $h_w^n$ is the height of residual water calculated in the previous physical time step, ${\dot{m}}_{in}$ and ${\dot{m}}_{out}$ are respectively the inflow and outflow of liquid water in the control volume.

In Equation (20), ${\dot{m}}_{imp}$ is the water droplet collection rate. ${\dot{m}}_{evp}$ is the evaporation rate, which is obtained by analogy with the convection heat transfer term. ${\dot{m}}_{ice}$ is the freezing rate. The unknown terms are specified as follows [36]:

$$\begin{array}{cc}\hfill & {\dot{m}}_{imp}=\mathrm{LWC}·U_{\infty }·\beta ,\hfill \\ \hfill & {\dot{m}}_{evp}=\frac{htc}{C{p}_a{\rho }_a{R}_{v}Le{w}^{2/3}}[\frac{e_s\left(T_w\right)}{T_w}-\frac{e_s\left(T_{\infty }\right)}{T_{\infty }}],\hfill \\ \hfill & {\dot{m}}_{ice}=f({\dot{m}}_{imp}+{\dot{m}}_{in}).\hfill \end{array}$$

(33)

where $\beta$ is the droplet collection coefficient. $htc$ is the convective heat transfer coefficient. $C{p}_a$ is the specific heat capacity of air. $Lew$ is the Lewis criterion number, which is related to the $\mathrm{LWC}$, with a general value of 1. $R_v$ is the gas constant of the water vapor. $e_s\left(T_w\right)$ and $e_s\left(T_{\infty }\right)$ are, respectively, the saturated vapor pressures corresponding to the internal surface temperature $T_w$, and $T_{\infty }$. f is a dimensionless freezing parameter, which is used to characterize the ratio of frozen water to entering water in the control volume.

$$\begin{array}{cc}\hfill & f=\frac{C{p}_w}{Lf}·(A+\frac{B·htc}{C{p}_w·U_{\infty }·\mathrm{LWC}·\beta }),\hfill \\ \hfill & A=-(T_{\infty }-T_w)-\frac{\left|\right|{\mathbf{U}}_d{\left|\right|}^2}{2C{p}_w},\hfill \\ \hfill & B=T_w-T_{\infty }-0.89·\frac{U_d}{2C{p}_w}+\frac{\lambda ·C·Le}{C{p}_a·htc},\hfill \\ \hfill & C=(\frac{e_s\left(T_w\right)}{T_{\infty }}-\frac{P_t·e_s\left(T_{\infty }\right)}{P_{\infty }}·T_t)/(\frac{P_t}{0.622T_t}-\frac{T_w}{T_{\infty }}).\hfill \end{array}$$

(34)

In which $Le$ and $Lf$ are the latent heat of water and ice. $C{p}_w$ is the specific heat capacity of water. ${\mathbf{U}}_d$ is the velocity of the impinging droplets. $\lambda$ is the the thermal conductivity of air. $T_t$ is the total temperature of the incoming stream, and $P_t$ is the total pressure of the incoming stream, which can be specified as follows:

$$\begin{array}{cc}\hfill & T_t=T_{\infty }·(1+\frac{\gamma -1}{2}·M{a}^2),\hfill \\ \hfill & P_t=P_{\infty }·{(1+\frac{\gamma -1}{2}·M{a}^2)}^{\frac{\gamma }{\gamma -1}}\hfill \end{array}$$

(35)

where $\gamma$ is the heat capacity ratio of air and $Ma$ is the inlet Mach number.

The energy conservation equation in the control volume can be written as follows:

$$q_{cond}=(q_{conv}+q_{evp}+q_{rad}+q_{war})-(q_{air}+q_{kin}+q_{fre}).$$

(36)

where $q_{cond}$ is the heat flux of heat conduction between the water film and the ice layer. $q_{evp}$ is the energy density taken away by evaporation or sublimation on the surface of the water film. $q_{rad}$ is the radiation heat flux of the water film. $q_{war}$ is the heat flux of heating. $q_{conv}$ is the convective heat transfer density between the water film and the air flow. $q_{air}$ is the heat flux of pneumatic heating. $q_{kin}$ is the heating heat flux converted from kinetic energy when water droplets impact. $q_{fre}$ is the latent heat flux of liquid water freezing. The unknown terms are specified as follows:

$$\begin{array}{cc}\hfill & q_{cond}={(-{\lambda }_w\partial T_w/\partial y)}_{y=h_w},\hfill \\ \hfill & q_{rad}={\sigma }_r(T_w^4-T_{\infty }^4),\hfill \\ \hfill & q_{conv}-q_{air}=htc(T_w-T_{rec}),\hfill \\ \hfill & q_{war}-q_{kin}=U_{\infty }·\mathrm{LWC}·\beta ·[C{p}_w(T_w-T_d)-\frac{\left|\right|{\mathbf{U}}_d{\left|\right|}^2}{2}],\hfill \\ \hfill & q_{fre}={\dot{m}}_{ice}C{p}_i{T}_w.\hfill \end{array}$$

(37)

In which $T_d$ is the temperature of the impinging droplets. $T_{rec}$ is the boundary layer recovery temperature.

In icing simulations, the convective heat transfer process $htc$ has a great influence on the ice shape and is the key to solving the convective heat transfer process. At present, most icing simulations use the boundary layer integral functions proposed by LEWICE, which consider the surface roughness and velocity variation to solve the convection heat transfer coefficient [37].

$$\begin{array}{c}\hfill \begin{array}{c}\hfill htc=\left\{\begin{array}{cc}0.296·\frac{{\lambda }_s}{{\nu }_a^{0.5}}{\left(V_e^{-2.87}{\int }_0^s{V}_e^{1.87}ds\right)}^{-0.5},\hfill & R{e}_k<R{e}_c\hfill \\ \frac{C_f/2}{Prt+\sqrt{C_f}/S{t}_k}·{\rho }_a·V_e·C{p}_a,\hfill & R{e}_k\ge R{e}_c\hfill \end{array}\right.\end{array}\end{array}$$

(38)

where ${\lambda }_s$ is the thermal conductivity of air. ${\nu }_a$ is the kinematic viscosity of air. $V_e$ is the velocity of air at the outer boundary of the boundary layer. s is the surface distance. $Prt$ is a Prandtl number, which is approximated as a constant 0.9. $R{e}_k$ and $R{e}_c$ are, respectively, the Reynolds number at the rough surface and the critical Reynolds number. $C_f$ is the friction coefficient between the air flow and the wall surface. $S{t}_k$ is a rough Stanton number, which can be obtained by

$$S{t}_k=0.8·{[\frac{Ve·k_s{C}_f}{2{\nu }_a}]}^{-0.2}Pr{t}^{-0.44},$$

(39)

where $k_s$ is the roughness height, which can be obtained by Shin’s improved equivalent sand roughness height model [38].

The liquid water flow on the surface is considered in SWIM, and the water film flow velocity is related to the water film thickness. In order to ensure the operation of the calculation, it is necessary to assume that there is a very thin water film on the surface at the initial time. The traditional SWIM does not consider the real situation of the initial water film, and manually sets the initial water film thickness as a constant (for example, assuming that the initial water film thickness is ${10}^{-6}$ m), and solves SWIM on this basis.

As the physical time step of the unsteady icing numerical simulation method used in this paper is small, it can be considered that at the initial time, the icing surface will form a water film and then freeze into ice, and the initial liquid water mass is the difference between the mass of the impact droplets and the amount of evaporation:

$${\left({\rho }_w{h}_w\right)}^0=\Delta t·({\dot{m}}_{imp}+{\dot{m}}_{in}-{\dot{m}}_{evp}-{\dot{m}}_{out}).$$

(40)

For quasi-steady icing calculations, due to the long physical time step, if the initial water film is simply used, it will lead to the problem that the initial water film is too thick and the water film flow velocity is too large, which will affect the calculation of icing shape. Therefore, for the quasi-steady icing calculation, the initial water film of each physical time step needs to be obtained before solving the SWIM, which is assumed in this paper that there is a short period of time at the initial time during which the surface liquid water does not freeze, and then the surface liquid water content distribution is obtained. However, for unsteady icing calculations, except that the initial water film of the first physical time step is obtained using the above method, the initial water film of each subsequent physical time step can be obtained by inheriting the calculation results of the previous physical step.

3. Verification of the Impact Characteristics

Before comparing the ice shapes, the accuracy of the tightly-coupled method to calculate the droplet collection coefficient should be verified. The summary of the icing aerodynamic database is shown in

Table 1. Icing conditions for the SSD (Supercooled Small Droplet)-case and SLD (Supercooled Large Droplet)-case.

Verification Case	Pressure Pa	MVD μm	Chord m	AOA deg	Velocity m/s	LWC g/m3	Temperature °C
SSD-case	101325	16	1	5	138.88	1.0	17
SLD-case	101325	154	0.9144	2.5	78.23	1.44	25.9

.

Figure 7 shows the contour of droplet volume fraction around the airfoil. Due to the angle of attack, the impact range of droplets on the lower surface of the airfoil is larger than that on the upper surface. The reference results from experiment [39], FENSAP and the Lagrangian method are used for comparison with the result from the tightly-coupled method in the SSD case, which is shown in Figure 8.

As can be seen from the figure, the results from the tightly-coupled method are essentially coincident with the experimental results, and are similar to those provided by the Eulerian method. The maximum droplet collection coefficient calculated by the tightly-coupled method and other methods are shown in

Table 2. The maximum collection coefficient.

	Maximum Collection Coefficient
Experiment	0.5695
Fensap	0.5339
Lagrangian	0.5715
Tightly-coupled	0.5648

. It can be seen from the figure that the maximum droplet collection coefficient predicted by the tightly-coupled model is the closest to the experimental value.

Figure 9 shows the comparison of the results provided by the icing codes and the experiment in the SLD case. In this case, the normal model indicates that the classical spherical resistance model is adopted without considering the dynamic mechanical effect of droplets, while the SLD model indicates that the droplet resistance coefficient model considering deformation [29] is adopted, and the crushing [40] and splashing [41] models are introduced.

According to the droplet diameter and droplet velocity, the broken droplet diameter and final droplet collection efficiency are modified by the crushing model and splashing model, respectively. It can be seen that the droplet collection coefficient from the SLD-tightly-coupled method agrees with the experimental results. In conclusion, the droplet collection efficiency predicted by the tightly-coupled model is in good agreement with the experimental results, which verifies the correctness of the model.

4. Droplet Shadow Zone

In order to study the influence of the droplet shadow zone, a cylinder with a diameter of 1 cm was selected as the research object: the air velocity $\mathbf{U}$ = 40 m/s, the angle of attack AOA = 0°, the flight height h = 0 km, the ambient temperature T = 273.15 K, the liquid water content LWC = 1 g/m³, and the droplet diameter MVD = 40 μm. The calculation result and the experimental result of reference [42] are shown in Figure 10. It can be seen from the figure that the boundary of the sheltered area obtained by the numerical simulation is the same as the experimental results, which proves the accuracy of the tightly-coupled model in this paper in predicting the shadow zone of droplets, and also shows the correctness of the calculation results of water droplet trajectory.

Then, the NACA0012 airfoil is selected as the research object, with chord length c = 1 m, incoming Mach number Ma = 0.382, angle of attack AOA = 6°, flight height h = 0 km, ambient temperature T = 300 K, liquid water content LWC = 1 g/m³ and water droplet diameter MVD = 20 μm. Figure 11 shows the droplet shadow zone calculated with the Eulerian method and the droplet trajectory obtained by the Lagrangian method. It can be seen from the figure that the two areas are essentially the same and extend downstream for too far, which is inconsistent with the actual situation in reference [42]. Therefore, both methods cannot accurately describe the location and scope of the droplet shadow zone. Figure 12 and Figure 13 show, respectively, the cloud image of droplet concentration and the droplet shadow zone calculated by the tightly-coupled model in this paper. The blue area is the droplet shadow zone, and the liquid water content in this area is 0, which is far less than the calculation results of the Eulerian method and Lagrangian method. This is due to the diffusion effect of the droplet flow field in the tail region of the airfoil, and the droplet shadow zone is greatly reduced compared with the former two methods. The results show that the droplet shadow zone can reach 10% behind the trailing edge of the airfoil.

5. Verification of Ice Shape

In order to compare the experimental shapes and the shapes provided by different icing codes, a total of six ice shapes and the icing conditions were selected from the literature, which are listed in

Table 3. Icing conditions for the SSD-case and SLD-case.

Verification Case	Airfoil	Time s	MVD μm	Chord m	AOA deg	Velocity m/s	LWC g/m3	Temperature °C
Case1 [43]	NACA0015	600	19	0.3531	0	85.2	0.75	−9.9
Case2 [43]	NACA0012	360	30	0.5334	4	102.9	1.8	−10.8
Case3 [44]	NACA0012	600	20	0.5334	4	67.1	1	−13.3
Case4 [44]	NACA0012	300	20	0.5334	0	58.1	2.1	−9.7
Case5 [45]	NACA0012	804	70	0.5334	0	51	0.91	−19.6
Case6 [45]	NACA0012	336	160	0.5334	0	77	1.04	−19.2

. For the purposes of comparison, the stagnation thickness and horn angle were chosen to be the ice shape characteristics, which were determined for the reference and prediction ice shapes.

Through these characteristic values illustrated below in Figure 14, it is easy to calculate the deviations between the experimental shapes and the prediction shapes provided by different icing codes. At the same time, in order to verify the accuracy of the unsteady tightly-coupled model, the ice shapes calculated by the new model are compared with those calculated by the one-step LEWICE icing code based on the Messinger model.

5.1. Ice Shape Predication for Unsteady Model and Quasi-Unsteady Model

Before verifying the unsteady tightly-coupled icing model, it is necessary to clarify the difference between the unsteady model and the quasi-steady model. The comparison of ice shapes in Case 6 shows that the smaller the physical time step is, the closer the ice shape is to the experimental result, as shown in Figure 15. The results in

Table 4. Deviations of stagnation thickness.

	Stagnation Thickness	Deviation
Experiment	0.038
Unsteady-Tightly-Coupled	0.0364	−3.6%
Onestep-Tightly-Coupled	0.027	−29.5%
Sixstep-Tightly-Coupled	0.031	−19.8%
LEWICE	0.30	−20.4%

show that the stagnation thickness is closer to the experimental results with the decrease of the physical time step.

For the quasi-steady model, because it is a complete and independent icing calculation process in each physical time step, the amount of calculation will increase with the decrease of physical time step. Therefore, for the sake of computational efficiency, the physical time step of the quasi-steady model is usually not defined as too small. The unsteady model can maintain the computational efficiency even though the time step is very small.

In this paper, the physical time step in the unsteady model was selected as 1 s. Except that the first time step is a complete steady icing calculation process, the remaining time steps only need to calculate 20 steps. When the physical time step is 1 s, the unsteady calculation process of Case 6 only takes 1 h, while the quasi-steady calculation process takes 12 h.

5.2. Verification Analysis of Unsteady Tightly-Coupled Model

The icing conditions calculated in this section are shown in Table 3, in which Case 1 to Case 4 focus on the icing of SSD, and Case 5 and Case 6 focus on the icing of SLD. Figure 16 show the comparison of ice shapes from the experiment, unsteady tightly-coupled model and LEWICE 3.0 under SSD icing conditions.

Table 5. Stagnation thickness of the SSD cases.

	Case 1		Case 2		Case 3		Case 4
Experiment	0.036		0.026		0.023		0.0167
Unsteady-Tightly-Coupled	0.035	−0.7%	0.026	1.0%	0.024	7.8%	0.0144	−13.4%
LEWICE	0.024	−35.3%	0.027	6.4%	0.014	−37.9%	0.0125	−25%

and

Table 6. Horn angle of the SSD cases.

	Case 1				Case 2		Case 3
	Top Horn Angle		Bottom Horn Angle		Top Horn Angle		Top Horn Angle
Experiment	63.2°		−54.1°		34.5°		51.4°
Unsteady-Tightly-Coupled	59.6°	−5.7%	−59.9°	10.8%	57.04°	65.5%	50.8°	−2.7%
LEWICE	64.2°	1.5%	−67.1°	24.1%	76.3°	121.5%	46.3°	−10.1%

show, respectively, the deviation of stagnation thickness and horn angle provided by the unsteady tightly-coupled model LEWICE and the experiment.

This set of data demonstrate that the stagnation thickness and horn angle from tightly-coupled model are closer to the experiment than those from LEWICE. Combined with the ice shape comparison results, the ice shape predicted by the tightly-coupled model is more similar to the experimental ice shape than LEWICE.

Figure 17 shows the comparison of ice shapes from the experiment, unsteady Eulerian method and LEWICE under SLD icing conditions. Figure 17a shows the comparison of numerical icing simulation under Case 5 meteorological conditions. In this icing condition, the ice shapes predicted by LEWICE form a concave angle at the leading edge of the airfoil.

The ice shape at the leading edge of the airfoil obtained by the numerical simulation of the unsteady tightly-coupled icing model can be considered to have no concave angle formation, and the result is closer to the experimental result. Figure 17a shows the comparison of numerical simulation of icing under Case 6. Under this icing condition, a convex angle is formed at the leading edge of the airfoil by the three methods.

Compared with the calculation results of the LEWICE 3.0 ice model, the convex angle formed by the unsteady tightly-coupled model is larger, and the ice thickness at the stagnation point is more consistent with the experimental result. As shown in

Table 7. Stagnation thickness of SLD cases.

	Case 5		Case 6
Experiment	0.0614		0.0378
Unsteady-Tightly-Coupled	0.0443	−27.8%	0.0364	−3.6%
LEWICE	0.0408	−33.6%	0.0301	−20.8%

, the deviation of stagnation thickness provided by the tightly-coupled model and the experiment is smaller than that by LEWICE. The ice shape predicted by the unsteady tightly-coupled model is closest to the experimental ice shape. In particular, in Case 5, only the unsteady tightly-coupled model can capture the ice shape characteristics.

In conclusion, from the comparison of ice shape characteristics, the calculation results of the tightly-coupled model demonstrated higher accuracy than those of LEWICE, and the deviations with the experimental results were smaller.

5.3. Ice Shape Predication for Tightly-Coupled Model and Loosely-Coupled Model

Aircraft icing is an unsteady process of two-way interaction between water and air. The influence of droplets on the air will increase with time, and thus the droplet collection coefficients calculated by the Eulerian model and tightly-coupled model are different, and the droplet collection coefficient has an important influence on the icing ice shape, which leads to the difference of the final calculated ice shape and the change of aerodynamic force. In this paper, the Eulerian model was chosen as a loosely-coupled model.

The unsteady Eulerian model and unsteady tightly-coupled model were used to predict the ice shapes in Case 1 and Case 5. The comparisons between the ice shapes calculated by the above two icing models and the ice shapes from experiment and predicted by LEWICE [45] are shown in Figure 18 and Figure 19. Figure 20 is an enlarged view of the leading edge ice shape. Under Case 1, the ambient temperature is −9.9 °C and the ice shape predicted by the unsteady loosely-coupled model is quite different from that predicted by LEWICE, which is closer to the experimental shape.

Table 8. Stagnation thickness of Case1 and Case 5.

	Case 1		Case 5
	Stagnation Thickness	Deviation	Stagnation Thickness	Deviation
Experiment	0.036		0.0614
LEWICE	0.024	−35.3%	0.0408	−33.6%
Unsteady-Tightly-Coupled	0.035	−0.7%	0.0443	−27.8%
Unsteady-Eulerian	0.031	−13.2%	0.0366	−40.3%

shows that the stagnation thickness provided by the unsteady tightly-coupled model is closest to the experiment.

Under Case 5, the ambient temperature is −19.6 °C, the leading edge of the experimental ice shape is convex, but the ice shape predicted by the unsteady-Eulerian model forms a concave angle at the leading edge of the airfoil, just like the ice shape predicted by LEWICE. The ice shape calculated by the unsteady tightly-coupled icing model does not form a concave angle at the leading edge, which is closer to the experimental ice shape.

In addition, the ice shape predicted by the unsteady-Eulerian method is essentially consistent with that predicted by LEWICE. Table 8 shows that the stagnation thickness provided by unsteady tightly-coupled model is closer to the experiment than that predicted by the unsteady-Eulerian model and LEWICE.

In summary, when the temperature is low, the icing model has the smaller impact, and the tightly-coupled model has the larger impact. At high temperatures, SWIM has the larger impact, and the tightly-coupled model can further improve the accuracy on this basis.

The tightly-coupled model and the Eulerian model were used to calculate the droplet collection coefficient of Case 5. Figure 21 and Figure 22 show the calculation results at the initial time and t = 804 s, respectively. At the initial time, the results of the two models essentially coincide. With the increase of time, when t = 804 s, the difference between the two models is significant. Especially in the leading edge of the airfoil, the droplet collection coefficient calculated by the tightly-coupled model is significantly greater than that calculated by the Eulerian model.

6. Conclusions

In this paper, we proposed a numerical simulation method for unsteady tightly-coupled icing model. The air and droplets were regarded as a single fluid, and the unsteady calculation was realized by using the dual-time-step method. The shallow water icing model was improved to adapt to the unsteady ice settlement.

The droplet collection coefficients of NACA0012 airfoil under different icing conditions were calculated and analyzed using the tightly-coupled model. The results show that the tightly-coupled model could accurately predict the droplet collection coefficient, which was the most consistent with the experimental results compared with the Lagrangian method and FENSAP-ICE. At the same time, the tightly-coupled model could correctly predict the droplet shadow zone.

The ice shapes under different icing conditions, including SSD and SLD icing conditions, were predicted and analyzed. The results show that, with the decrease of the physical time step, the ice shape predicted by the unsteady model was closer to the experimental result, and the unsteady model maintained the computational efficiency compared with the quasi-steady model. As time increased, the cumulative effect of droplets on the airflow field gradually appeared, the difference between the results of the tightly-coupled model and the loosely-coupled model increased, and the results of the tightly-coupled model were more consistent with the experimental ice shapes. When the temperature was low, there was little difference between the loosely-coupled model and LEWICE in predicting ice shape, indicating that the influence of icing model was small and the tightly-coupled model had the larger impact. At high temperatures, SWIM had the larger impact, and the tightly-coupled model could further improve the accuracy on this basis.

Through comparative analysis of the numerical simulations and the experimental results, the ice shapes predicted by the unsteady tightly-coupled model were more consistent with the experimental ice shapes under different icing conditions. Compared with LEWICE, the unsteady tightly-coupled model had higher accuracy. For glaze ice, especially horn ice, where the existing icing methods, such as LEWICE, cannot accurately capture the ice shape characteristics, the unsteady tightly-coupled model could more accurately predict the ice shapes and successfully capture the ice shape characteristics.