2.1 Air flow evaluation

In this work, Reynold Average Navier-Stokes (RANS) equations are used as the air flow control equations:

(1) \begin{align}\begin{array}{l}{\rm{ }}\nabla \cdot {{\boldsymbol{U}}} = 0{\rm{ }}\\[5pt] \nabla \cdot \left( {{{\boldsymbol{UU}}}} \right) = - \nabla \frac{p}{\rho } + \nabla \cdot \left( {\nu \nabla {{\boldsymbol{U}}} - {{\boldsymbol{R}}}} \right)\end{array}\end{align}

where $\rho$ , U and p are the density, velocity vector and pressure of air, respectively, and R is Reynold stress tensor. The steady-state solver simpleFoam in OpenFOAM, based on SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm, combined with the k-ω-SST turbulence model, is used to solve the RANS equations.

2.2 Droplet impingement computation

The droplet movement and impingement are described through the Lagrangian method with the assumptions for droplets:

1. 1) No collision of droplets exists.

2. 2) Physical properties such as density and temperature are invariable.

3. 3) Droplets distributes evenly and are considered as spheres.

According to Newton’s second law, the movement equations of the droplets affected by the gravity and drag of the air are:

(2) \begin{align}{m_{\rm{d}}}\frac{{d{{{\boldsymbol{U}}}_{\rm{d}}}}}{{dt}} = \frac{1}{2}{\rho _{\rm{d}}}\left| {{{\boldsymbol{U}}} - {{{\boldsymbol{U}}}_{\rm{d}}}} \right|\left( {{{\boldsymbol{U}}} - {{{\boldsymbol{U}}}_{\rm{d}}}} \right){A_{\rm{d}}}{C_D} + {m_{\rm{d}}}{{\boldsymbol{g}}}\end{align}

where ${m_{\rm{d}}}$ , ${{{\boldsymbol{U}}}_{\rm{d}}}$ and ${\rho _{\rm{d}}}$ are the mass, velocity vector and density of the droplets respectively. ${A_{\rm{d}}}$ is the windward area of the droplets. g is the gravity acceleration vector, drag coefficient C _D is acquired via the correlation expression:

(3) \begin{align}\frac{{{C_D}R{e_{\rm{d}}}}}{{24}} = 1.0 + 0.197Re_{\rm{d}}^{0.63} + 2.6 \times {10^{ - 4}}Re_{\rm{d}}^{1.38}\end{align}

Reynold number is:

(4) \begin{align}R{e_{\rm{d}}} = \rho \left| {{{\boldsymbol{U}}} - {{{\boldsymbol{U}}}_{\rm{d}}}} \right|\frac{{{d_{\rm{d}}}}}{\mu }\end{align}

where d _d is the diameter of the droplets.

In order to get a realistic cloud, the injection type of droplets is set to “coneInjection” in the “reactingCloud1Properties” file in the “constant” directory, which can disperse the droplets with a limited number of injection positions. The injected droplets move under the action of gravity and air flow, and eventually impinge onto the wing or escape from it. We can then obtain the local collection efficiency:

(5) \begin{align}\beta = \frac{{d{A_0}}}{{dA}}\end{align}

in which $d{A_0}$ is an infinitesimal area of far-field droplets, and $dA$ is the corresponding impingement area on the wing, as shown in Fig. 1.

Figure 1. The schematic of droplet impingement onto the wing.

2.3 Film flow simulation

The thin water film on the wing surface is formed by the impingement of droplets and flows downstream driven by gravity and air shear force. The thin film (with a thickness less than 0.1mm) is assumed to be laminar. The thin film flow model is based on the thin-film assumption [Reference Meredith, Xin and Vries33]: 1) pressure in the normal direction is constant; 2) normal gradients are considered other than the much smaller tangential gradients; 3) tangential velocity is considered other than the much smaller normal velocity.

2.3.1 Mass continuity equation

The continuity equation of film flow is:

(6) \begin{align}\frac{{\partial h}}{{\partial t}} + \nabla \cdot \left( {h{{{\bar{\boldsymbol{U}}}}_l}} \right) = \frac{{\sum {{m_{\rm{d}}}} }}{{{\rho _l}dtdS}}\end{align}

where h is the thickness of water film, ${\rho _l}$ is the density of the liquid film, average velocity of water film flow is,

(7) \begin{align}{{\bar{\boldsymbol{U}}}_l} = \frac{1}{h}\int_0^h {{{{\boldsymbol{U}}}_l}dz} \end{align}

$\frac{{\sum {{m_{\rm{d}}}} }}{{dtdS}}$ is defined as the mass change of the water film per unit area and per unit time, due to droplet impingement and splash. The splashing model is proposed by Bai [Reference Bai, Rusche and Gosman18].

2.3.2 Momentum transport equation

The momentum equation is:

(8) \begin{align}\frac{\partial }{{\partial t}}\left( {h{{{\bar{\boldsymbol{U}}}}_l}} \right) + {\nabla _{\rm{s}}} \cdot \left( {h\overline {{{{\boldsymbol{U}}}_l}{{{\boldsymbol{U}}}_l}} } \right) = - \frac{1}{{{\rho _l}}}h{\nabla _{\rm{s}}}{p_l} + \frac{1}{{{\rho _l}}}\left( {{{\boldsymbol{\tau }}_{\rm{f}}} - {{\boldsymbol{\tau }}_{\rm{w}}}} \right) + h{{{\boldsymbol{g}}}_{\rm{s}}} + \frac{{\sum {{m_{\rm{d}}}{{{\boldsymbol{U}}}_{{\rm{d}},{\rm{s}}}}} }}{{{\rho _l}dtdS}}\end{align}

where ${\nabla _{\rm{s}}}$ is the vector differential operator tangential to the surface, g _s represents the tangential component of gravity acceleration, $\frac{{\sum {{m_{\rm{d}}}{{{\boldsymbol{U}}}_{{\rm{d}},{\rm{s}}}}} }}{{dtdS}}$ is the tangential momentum change of water film per unit area and per unit time, and ${{{\boldsymbol{U}}}_{{\rm{d}},{\rm{s}}}}$ is the tangential velocity of droplets. The velocity profile in the normal direction is assumed quadratic, then we can get the correlation relationship of $\overline {{{{\boldsymbol{U}}}_l}{{{\boldsymbol{U}}}_l}} $ in convection term:

(9) \begin{align}\overline {{{{\boldsymbol{U}}}_l}{{{\boldsymbol{U}}}_l}} = \frac{1}{h}\int_0^h {{{{\boldsymbol{U}}}_l}{{{\boldsymbol{U}}}_l}dz} = \frac{{{h^2}{\boldsymbol{\tau }}_{\rm{f}}^2}}{{120\mu _l^2}} + \frac{{h{{\boldsymbol{\tau }}_{\rm{f}}}{{{\bar{\boldsymbol{U}}}}_l}}}{{20{\mu _l}}} + \frac{6}{5}{\bar{\boldsymbol{U}}}_l^2\end{align}

The pressure of the liquid film is:

(10) \begin{align}{p_l} = {p_{\rm{a}}} + {p_\sigma } + {p_h} + {p_{\rm{d}}}\end{align}

where ${p_{\rm{a}}}$ is the air pressure, and the pressure caused by surface tension is:

(11) \begin{align}{p_\sigma } = - \sigma {\nabla _{\rm{s}}} \cdot \left( {{\nabla _{\rm{s}}}h} \right)\end{align}

The hydrostatic pressure is

(12) \begin{align}{p_h} = - {\rho _l}{{\boldsymbol{n}}} \cdot {{\boldsymbol{g}}}h\end{align}

where n is the unit normal vector, and the normal dynamic pressure arising from impingement and splash is:

(13) \begin{align}{p_{\rm{d}}} = \frac{1}{2}{\rho _{\rm{d}}}{{\boldsymbol{U}}}_{{\rm{d}},{\rm{n}}}^2\end{align}

where ${{{\boldsymbol{U}}}_{{\rm{d}},{\rm{n}}}}$ is the normal velocity of droplets. The film-gas interface shear stress ${{\boldsymbol{\tau }}_{\rm{f}}}$ and the film-wall interface shear stress ${{\boldsymbol{\tau }}_{\rm{w}}}$ are, respectively:

(14) \begin{align}{{\boldsymbol{\tau }}_{\rm{f}}} = {\mu _l}{\left. {\frac{{\partial {{{\boldsymbol{U}}}_l}}}{{\partial z}}} \right|_{z = h}}\end{align}

(15) \begin{align}{{\boldsymbol{\tau }}_{\rm{w}}} = {\mu _l}{\left. {\frac{{\partial {{{\boldsymbol{U}}}_l}}}{{\partial z}}} \right|_{z = 0}} = - \frac{{{{\boldsymbol{\tau }}_{\rm{f}}}}}{2} + 3{\mu _l}\frac{{{{{\bar{\boldsymbol{U}}}}_l}}}{h}\end{align}

where the film-gas interface shear ${{\boldsymbol{\tau }}_{\rm{f}}}$ can be got from the air flow field solver.

Figure 2. The solution procedure of droplet impingement and film flow using OpenFOAM.

Figure 3. The impingement of droplets onto 2D MS-0317 aerofoil and the contours of droplet dimeters (with unit μm), with free stream velocity of 78.66m/s and MVD of 16μm.

Figure 4. Comparison of the local water collection efficiency of 2D MS-0317 aerofoil with free stream velocity of 78.66m/s and MVD of 16μm.

Figure 5. The impingement of droplets onto 15.04cm diameter sphere and the contours of droplet dimeters (with unit μm), with free stream velocity of 75m/s and MVD of 18.6μm.

Figure 6. Comparison of the local water collection efficiency of 15.04 cm diameter sphere with free stream velocity of 75m/s and MVD of 18.6μm.

Figure 7. Comparison of the film thickness with Nusselt solution for steady continuous water film flow over a plane with inclined angles of 5° and 90°.

Figure 8. Comparison of the film velocity with Nusselt solution for steady continuous water film flow over a plane with inclined angles of 5° and 90°.

Figure 9. Comparison of the pressure coefficient contours of three grids with different cell numbers with free stream velocity of 67m/s and AOA of 0°. (a) The pressure coefficient contours on the symmetry and wall surfaces of the grid with 277,000 cells. (b) The pressure coefficient contours on the symmetry and wall surfaces of the grid with 602,000 cells. (c) The pressure coefficient contours on the symmetry and wall surfaces of the grid with 897,000 cells.

Figure 10. Comparison of the pressure coefficient distributions on the middle section of the half-wing of three grids with different cell numbers with free stream velocity of 67m/s and AOA of 0°.

Figure 11. The surface mesh of MS-0317 straight wing, composed of quadrilateral faces.

Figure 12. The droplets movement around three-dimensional MS-0317 straight wing, covering the region near the spanwise midpoint of the half wing.

Figure 13. The time history of the water film thickness at the leading-edge point of three-dimensional MS-0317 straight wing in 20s, with AOA of 0°, 4°, 8° and 12°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

The film flow control equations are solved using the Finite Area Method (FAM) [Reference Tuković26], which is suitable for solving partial differential equations on curved surfaces. In OpenFOAM, FAM is implemented on a two-dimensional surface mesh that is tangential to the wing surface. The mesh is only one layer normal to the surface and is extruded from the wing boundary mesh. This approach separates the film flow from the air flow, making it more resource efficient than the interface capture method.

2.4 Interaction of the air flow, droplet movement and thin water film flow

The computational process of droplet impingement and film flow is illustrated in Fig. 2. It is considered as a one-way coupling between air flow and droplet movement: the air flow influences the trajectories of the droplets, while the effects of the droplets on the variables of the air flow are negligible. A one-way coupling is also assumed between air flow and water film flow: the air flow shear force drives the water film to flow downstream, without considering the influence of the water film on the air flow. It is considered to be a two-way coupling between droplet movement and film flow: the impingement of droplets onto water film changes the mass and momentum of the water film, while the droplets from the splash return to the air and are tracked again. In OpenFOAM, both droplet movement and water film flow can be simulated using reactingParcelFoam.

3.1 Local collection efficiency

To validate the computation of droplet movement and impingement, the local collection efficiency of a two-dimensional MS-0317 aerofoil and a three-dimensional sphere is calculated and compared with the data from experiments and LEWICE.

3.2 Nusselt solution for film flow

In this section, the Nusselt solution [Reference Nusselt36] is used to verify the thin film flow model. The film thickness and velocity of steady continuous water film flow over a plane with two inclined angles (5° and 90°) are computed and compared with Nusselt solution respectively. The temperature is 298K.

The thickness and velocity of the Nusselt solution have the theoretical correlations with Reynold number respectively as shown below. Comparisons of solution in the present numerical simulation with that of the theoretical Nusselt solution are shown in logarithmic scale in Figs. 7 and 8, which show good agreement.

(16) \begin{align}h = {\left( {\frac{{3{\nu ^2}}}{{g\sin \theta }}} \right)^{1/3}}R{e_f}^{1/3}\end{align}

(17) \begin{align}U = {\left( {\frac{{\nu g\sin \theta }}{3}} \right)^{1/3}}R{e_f}^{2/3}\end{align}

where $\theta $ is the inclined angle.

4.1 Grid independence study

Firstly, for the three-dimensional wing, we use three grids of different cell numbers to investigate the grid independence, with cell numbers of 277,000, 602,000 and 897,000. The three-dimensional wing is obtained by straight extrusion of the MS-0317 aerofoil in the spanwise direction, and the semispan is 1.8288m. The flow velocity is 67m/s, AOA = 0°. Figure 9 shows the pressure coefficient contours on the symmetry and wall surfaces of three grids with 277,000, 602,000 and 897,000, and Fig. 10 is the comparison of the pressure coefficient distributions on the middle section of the half-wing of the three grids. We can see that both the pressure coefficient contours and the distributions of the two grids with 602,000 and 897,000 cells are almost the same. The grid with 602,000, cells is therefore suitable. Figure 11 shows the wing surface mesh of the grid with 602,000 cells.

4.2 Three-dimensional MS-0317 straight wing with different AOAs

The characteristics of the film flow with AOA = 0°, 4°, 8° and 12° are evaluated and compared, respectively. The flow velocity is 67m/s, LWC = 1.03g/m³, MVD = 20μm, and the temperature is 300K. Here the region near the spanwise midpoint of half wing (0.7144m ∼ 1.1144m) is covered by the droplets, as shown in Fig. 12. Without covering the entire wing, firstly the computation cost is reduced greatly, and secondly enough interesting information can be acquired with this region.

Table 1. The maximum value of the film thicknesses with AOA of 0°, 4°, 8° and 12° at the leading-edge point of the spanwise middle section of half wing

Figure 14. The thickness distribution (with unit μm) of the water film on upper surface (left, wing root on top) and lower surface (right, wing tip on top) of the three-dimensional MS-0317 straight wing at 20s, with AOA of 0°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 15. The thickness distributions of the water film on the middle section of three-dimensional MS-0317 straight wing at different time, with AOA of 0°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 16. The chordwise velocity distributions of the water film on the middle section of MS-0317 straight wing at different time, with AOA of 0°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 17. The thickness distribution (with unit μm) of the water film on upper surface (left, wing root on top) and lower surface (right, wing tip on top) of three-dimensional MS-0317 straight wing at 20s, with AOA of 8°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 18. Comparison of the spanwise velocity of the water film on the middle section of three-dimensional MS-0317 straight wing with different AOAs at 20s, with the same free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 19. The spanwise flow velocity of the air near the middle surface of three-dimensional MS-0317 straight wing, with AOA = 8°, free stream velocity of 67m/s.

Figure 20. Comparison of the thickness of the water film on the middle section of three-dimensional MS-0317 straight wing at different time, with different AOAs, with the same free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 21. Comparison of the chordwise velocity of the water film on the middle section of three-dimensional MS-0317 straight wing at different time, with different AOAs, with the same free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 22. The time history of the water film thickness at the leading-edge point of MS-0317 swept wing in 20s, with swept-back angle of 0°, 10°, 20° and 30°, AOA of 0°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 23. The thickness distribution (with unit μm) of the water film on upper surface (left, wing root on top) and lower surface (right, wing tip on top) of three-dimensional MS-0317 swept wing (swept-back angle = 30°) at 20s, with AOA of 0°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 24. Comparison of the spanwise velocity of the water film on the middle section of MS-0317 straight and swept wing (swept-back angle = 30°) at 20s, with AOA of 0°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 25. The spanwise flow velocity of the air near the middle surface of three-dimensional MS-0317 swept wing (swept-back angle = 30°), with AOA = 0°, free stream velocity of 67m/s.

In all simulations, the total computation time is 20s, while the injection duration of droplets is 2s from the start. Figure 13 is the time history of the water film thickness at the leading-edge point of the spanwise middle section of half wing with AOA of 0°, 4°, 8° and 12°. We can see that the thickness gets to the maximum and keeps constant during the injection of the droplets and then the thickness reduces affected by the air shear force once the injection stops but seems to be nonzero at last, and the values eventually converge to the same. This is because during the injection, the transport mass of the thin film flow will equal to the mass added by the droplet impingement in the end, and without injection, the wall shear force eventually gets to equilibrium with the air shear force. As the AOA increases, the water film thickness gradually decreases. The value of the maximum thicknesses with AOA of 0°, 4°, 8° and 12° are shown in Table 1. Although the maximum thicknesses are different, both the history curves of the water film thickness before and after the end of the injection at the leading-edge point can be fitted with logarithmic functions. The following is a detailed study of the water film flow characteristics with AOA = 0° and 8°.

4.3 MS-0317 swept wings with different swept-back angles

For the MS-0317 swept wing, the characteristics of the film flow with swept-back angle = 0°, 10°, 20° and 30° are studied. All the conditions are the same with the first case “AOA = 0°”. Figure 22 is the time history of the water film thickness at the leading-edge point of the spanwise middle section of half wing with swept-back angle of 0°, 4°, 8° and 12°. As the swept-back angle increases, the water film thickness gradually decreases. The history curves of the water film thickness before and after the end of the injection at the leading-edge point can also be fitted with logarithmic functions.

Here is a detailed study of the water film flow characteristics with swept-back angle = 30°. Figure 23 is the thickness distribution of the water film on upper surface (left) and lower surface (right) of the MS-0317 swept wing at 20s. The thickness distribution is also asymmetric particularly on the upper surface, and the film on the upper surface flows outward as a whole.

The comparison of the spanwise velocity of the water film on the middle section between the MS-0317 straight and swept wing with swept-back angle = 30° is shown in Fig. 24. For the swept wing, the film flows outward near the leading edge, the spanwise velocity of the film firstly reduces downstream on the both sides to zero, then the film flows inward, the spanwise velocity increases to the maximum and finally reduces to zero.

The spanwise velocity is caused by the three-dimensional effects of the MS-0317 swept wing. The air velocity around the wing firstly slants to the wing tip, then to the wing root, and finally to the wing tip namely “S-type” streamlines due to the angle of sweepback. Figure 25 is the spanwise velocity around the middle section of the MS-0317 swept wing.

Figure 26 is the comparison of the thickness of the water film on the middle section between the MS-0317 straight and swept wing at different time. Because of the spanwise flow of the water film, we can see that the thickness of the MS-0317 swept wing is smaller than that of the straight one. Figure 27 is the comparison of the chordwise velocity of the water film between the MS-0317 straight and swept wing.

Figure 26. Comparison of the thickness of the water film on the middle section of the MS-0317 straight and swept wing (swept-back angle = 30°) at different time, with AOA of 0°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.

Figure 27. Comparison of the chordwise velocity of the water film on the middle section of MS-0317 straight and swept wing (swept-back angle = 30°) at different time, with the same AOA of 0°, free stream velocity of 67m/s, LWC of 1.03g/m³ and MVD of 18.6μm.