A Computational Modeling for Semi-Coupled Multiphase Flow in Atmospheric Icing Conditions

AbstrAct: A two-dimensional second-order positivitypreserving finite volume upwind scheme is developed for a semi-coupled algorithm involving the air and droplet flow fields in the Eulerian frame, which shares the grid for each phase. Special emphasis is placed on the computational modeling, which is induced from a strongly-coupled algorithm that satisfies the strict hyperbolicity and its numerical scheme based on the Harten-Lax-van Leer-Contact solver preserving the positivity to handle multiphase flow in the Eulerian frame. The proposed modeling associated with the semi-coupled algorithm including the Navier-Stokes and droplet equations takes into account different boundary conditions on the solid surface for each phase. The verification and validation studies show that the new scheme can solve the air and droplet flow fields in fairly good agreement with the exact analytical solutions and experimental data. In particular, it accurately predicted the maximum value of the droplet impingement intensity near the stagnation region and the droplet impingement area.


IntroductIon
In the atmosphere, an in-flight aircraft icing occurs when supercooled water droplets freeze on impact with any part of the external structure of an aircraft during flight in icing conditions, such as a cloud (Lynch and Khodadoust 2001).Icing also occurs on wind turbine blades exposed to low temperatures and water droplets (Parent and Ilinca 2011).Such flows can be found in the air-mixed supercooled water droplet fields.Exposure to this condition for a considerable period can be extremely dangerous to aircraft, as the built-up ice can distort the smooth flow of air over the wing.This serious aircraft safety concern can be involved in diminishing the wing's maximum lift coefficient, reducing the angle of attack for maximum lift coefficient, adversely affecting aircraft handling qualities and significantly increasing the drag (Lynch and Khodadoust 2001;Bragg et al. 2005).In addition, aircraft accidents due to icing have been occasionally reported in articles.Representatively, the ATR 72 during an En-route is crashed near Roselawn, Indiana, USA, by the loss of controllability due to icing in 1994 (Gent et al. 2000).A prediction how the water droplets impinge on external structure in such cases can be essential to design proper ice protection systems that prevent or remove ices accumulated on these structures (Kind et al. 1998;Ahn et al. 2015).The design of ice protection systems requires knowledge of the local and total impingement intensities in order to determine the heat inputs for removing built-up ice and the extent of the surface area to be protected.Useful data on the local and total impingement intensities can be obtained via icing wind tunnel tests, or state-of-the-art computational fluid dynamics (CFD) simulations.Although the icing wind tunnel tests can provide detailed information, they are very expensive and there are scaling limitations of test models due to wind tunnel sizes in most cases.On the other hand, the CFD-based simulations are relatively cheap and can handle icing problems without limitations (Anderson 2000).Therefore, CFD-based simulations have been increasingly spread and have found their way into the main stream of computational methods for obtaining information about droplet impingement.The numerical simulations to predict droplet impingements have been traditionally based on the two approaches of modeling particle transport in CFD simulations, the Eulerian and the Lagrangian approaches.The Eulerian approach treats the particle phase as a continuum and develops its conservation equations on a control volume basis in a similar form as that for the fluid phase.The Lagrangian approach considers particles as a discrete phase and tracks the pathway of each individual particle, such as the NASA LEWICE and ONERA codes (Hamed et al. 2005;Wright et al. 1997;Morency et al. 1999;Silveira et al. 2003;Vu et al. 2002).In the Lagrangian formulation, the trajectory of each particle is computed using a force balance equation.Meanwhile, the Eulerian approach is more flexible for three-dimensional complex geometries.
In general, multiphase flows are simulated using the strong coupled algorithm in which each phase is closely interacted.Nevertheless, the multiphase flows around the aircraft flying into the cloud can be simulated using a weakly coupled algorithm, because the mass loading ratio of the bulk density of the droplets over the bulk density of air is on the order of 10 -3 under icing conditions.It means that the effect of droplets in the air flow can be ignored.A number of phenomena and forces may be considered while modeling air-droplets flows, but the following assumptions are sensible in an aircraft icing: the spherical droplets without any deformation or breaking; no droplets collision, coalescence or splashing; no heat and mass exchange between the droplets and the surrounding air; turbulence effects on the droplets can be neglected; the only air drag, gravity and buoyancy forces acting on the droplets (Bourgault et al. 2000).This observation can justify a weakly coupled algorithm in which the calculations for each phase can be independently conducted for the atmosphere icing conditions mixed with the air and a supercooled water droplet.This assumption is however valid for steady state simulations.For unsteady state simulations, such as a moving body, the weakly coupled algorithm may not be valid since, at every time step, air solver should provide the time-dependent primitive variables, such as density and velocity of air flow, to the droplet solver taking into account the drag and buoyancy of droplet.In such cases, multiphase flows around the aircraft may be treated in strongly-coupled or semi-coupled manners so that the equations of each phase are interactively solved or the air flow affects the droplet flow strongly.For the strongly-coupled manner, the full system of equations is based on the strong interactions between the phases divided by volume fractions.Interestingly, the droplet flow does not affect the air flow due to the aforementioned assumption, mass loading ratio of the bulk density between each flow.It results that the multiphase flows on the aircraft icing do not require the strongly-coupled manner for unsteady simulations.In addition, the full system of equations of strongly-coupled algorithm is very complex and requires a lot of computational resource compared with semi-coupled algorithm.Consequently, the semi-coupled algorithm can be an alternative way to simulate the unsteady simulations.
At present, however, critical computational issues remain regarding the CFD-based air-mixed droplet simulations as a unified solver in which air and droplet solvers are integrated using the single numerical scheme.The issues concern the accuracy of numerical solutions to the droplet equations and its hyperbolic nature.In particular, the liquid water content (LWC) is observed around a solid surface after the droplets impinge on the surface (Jung and Myong 2013;Jung et al. 2011), as illustrated in Fig. 1.Moreover, it consists of various non-linear wave regions: shock fronts, rarefactions, contact discontinuities, and transitional layers at the interface of the shadow and non-shadow areas.Without proper positivity-preserving schemes, the LWC near the surface may become negative, resulting in the numerical breakdown (Kim et al. 2013).The reason is that the convective system of the droplet equations is not strictly hyperbolic.Even though the system has real eigenvalues, λ i = 1, 2, 3 = u, they are not distinct and it results that the system of equations does  not satisfy the hyperbolic conservation law.This means that the well-known numerical methods based on the well-posed strictly hyperbolic system such as the approximate Riemann solver may not be applicable, and the special schemes based on the kinetic approximation (Bouchut et al. 2003) are required.To deal with this issue, a two-dimensional second-order positivity-preserving finite volume upwind scheme only for the droplet equations was developed in the previous study (Jung and Myong 2013).The scheme is based on the Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver (Toro 2001;Toro et al. 1994) for the so-called shallow water droplet model (SWDM) (Jung and Myong 2013).The key component of the scheme is a simple technique based on splitting of the original mathematical system into a well-posed hyperbolic part and the source term in which a vector term involving the divergence of the artificial term is added and subtracted to the momentum equation for a purely numerical purpose (Jung and Myong 2013).
In this study, a new two-dimensional second-order positivitypreserving finite volume upwind solver for the semi-coupled algorithm in Eulerian frame has been developed using the well-posed hyperbolic part.It will be shown that this numerical scheme satisfies the density positivity condition.In addition, the second-order unified solver of the semi-coupled algorithm in aircraft icing conditions can serve as a basic concept, which is extendable to the moving body simulations.In "Governing Equations" section, the mathematical model of the semi-coupled algorithm for an aircraft icing from a compressible multiphase flow (strongly-coupled multiphase algorithm) introduced by Saurel and Abgrall (1999) is reviewed, and the well-posed hyperbolic part of the eigen structure of the system is presented.In "Numerical Methods" section, the second-order finite volume method based on the HLLC approximate Riemann solver is presented.In "Verification and Validation" section, numerical results of one-and two-dimensional test problems, as the verification and validation of the new scheme, are presented.

GovernInG equatIons
Atmospheric multiphase flow for an aircraft icing can be treated as a compressible and dilute multiphase flow ignoring the droplet effects in air flows since the mass loading ratio of the bulk density of the droplets over the bulk density of air is on the order of 10 -3 in flows.Although the flows can be simulated using a weakly-coupled algorithm for a steady state simulation, it exposes a limitation for an unsteady state simulation since all primitive variables at every time-marching step cannot be provided into the droplet solver, through the source terms of the droplet equations, due to enormous data size.Consequently, strongly-coupled or semi-coupled algorithms should be employed to simulate unsteady-depended aircraft icing.In this study, the semi-coupled algorithm is proposed by the aforementioned reason in introduction section.The convection-type equations of semi-coupled algorithm for an aircraft icing can be induced by the compressible multiphase flow (strongly-coupled algorithm), which is introduced by Saurel and Abgrall (1999).The full system of equations of compressible multiphase flow, except for the viscous effects of gas phase and the buoyancy and gravity forces acting on the droplet of liquid phase, is described as: where: superscripts g and l denote the gas and liquid phases, respectively; E is the total energy; P means static pressure; I is the identity matrix; V and α are an interfacial velocity and a volume fraction, α g + α l = 1; ρ and u are the density and velocity of each phase, respectively; µ is the viscosity; λ refers to the drag coefficient.The subscript i means an interface between two phases.Using the aforementioned assumptions (α g ~ 1, α g ρ g = ρ a , α l ρ l = ρ d = LWC, ρ l P l = 0 and V i = P i = 0), the system of equations are summarized as: (1) (2) where: -λ(u l -u g ) denotes a drag force of droplet due to the air flow; the subscripts a and d represent the air and the droplet, respectively.
Finally, the full system of equations involving the viscous effects of air flow and the buoyancy and gravity forces acting on the droplet is described as: where: a a is a sound of speed of air flow.The matrix A(W) has real eigenvalues, λ 1 = u a -a a , λ 2 = u a , λ 3 = u a + a a , λ 4 = λ 5 = u d , with corresponding right eigenvectors: (3) where: Here, the non-conserved variables τ and Q denote the viscous shear stress tensor and the heat flux, respectively.In Eq. ( 3), the symbol [Vu a ] (2) in the viscous shear stress tensor stands for the traceless symmetric part of the tensor; k is the thermal conductivity and depend on the air temperature.For the air flow, the ideal equation of state p = ρRT is used in order to close the equations.S b = ρg[0, 0, 1 -ρ a /ρ w ] T , being ρ a and ρ w the density of air and water, is the resultant force of the gravity and buoyancy of droplets.The A u (u a -u d ) denotes the drag acting on droplets caused by the airflow, and g is acceleration due to gravity.Also, the MVD is the mean volume diameter of the droplet.Re u and C Du are the Reynolds number of droplets and the drag coefficients of the spherical droplets, respectively.The drag coefficient can be obtained from Lapple (1951) as follows: which is valid for Re u < 1,000.In Eq. ( 3), the first three equations are the well-known Navier-Stokes equations for air flow and the last two equations are the droplet equations for droplet flow, respectively.

semi-coupled Algorithm with shAllow wAter droplet model for strict hyperbolicity
A shallow water droplet model based on the Eulerian framework developed in the previous work (Jung and Myong 2013) was employed for the present semi-coupled algorithm.In the present semi-coupled algorithm of Eq. ( 3), the strict hyperbolic conservation law is not satisfied due to a lack of distinctive eigen systems.It means that the well-posed numerical schemes based on the approximated Riemann solvers may not be applied to the equations.In this study, the eigen systems of convective fluxes of Eq. ( 3) are introduced using the primitive formulation of a one-dimensional semi-coupled algorithm in order to proof the lack of distinctive eigen systems: Even though the matrix A(W) has real eigenvalues, fourth and fifth eigenvalues and eigen vectors exactly corresponding to droplet equations are not distinctive.Thus, the formulation based on splitting of the droplet equations into the well-posed hyperbolic part and the source term, which circumvents the nonstrictly hyperbolic nature of the droplet equations, is applied.The modified full system of equations involving the previous study can be written as: and its primitive form with corresponding eigen systems is summarized as: where: A Computational Modeling for Semi-Coupled Multiphase Flow in Atmospheric Icing Conditions a d is (gD) 1/2 .D and g are a diameter named by the MVD of droplet and gravity, respectively.The matrix A(W) has real eigenvalues, λ 1 = u a -a, λ 2 = u a , λ 3 = u a + a, λ 4 = u d -a d , λ 5 = u d + a d , with corresponding right eigenvectors: There are six distinct waves for the well-posed hyperbolic part of Eqs. 6 and 7.The solutions in the Star region for each phase consist of two constant states separated from each other by a middle wave of speed.In the exact Riemann solver, the middle wave speed S * can be estimated for the particle velocity in the middle state (S * = u * ; S is a wave speed and u * is the x-directional particle velocity in the middle state).The well-posed hyperbolic part of governing equations in two-dimensional and x-spirit case can be rewritten as follows with a change of notation ψ = v (v is the y-directional particle velocity): The tangential velocity component ψ(x, t) represents the concentration of a pollutant or some other passive scalar.The convective flux of Eq. 8 has real and distinct eigenvalues, λ 1 = u a -a a , λ 2 = λ 3 = u a , λ 4 = u a + a, λ 5 = u d -a d , λ 6 = u d , λ 7 = u d + a d .The quantity ψ for each phase gives rise to the middle eigenvalue, λ 3 = u a and λ 6 = u d , respectively.For this hyperbolic system, the following HLLC flux for the approximated Godunov method can be defined at the interface of left and right cells as: Note that Eqs. 6 and 7 satisfy the hyperbolic conservation law.In the following section, a numerical method is described using the semi-coupled algorithm based on Eqs. 6 and 7.

numerIcal methods hllc ApproximAted riemAnn solver for the semi-coupled Algorithm
For the semi-coupled algorithm, a vacuum state caused by the rarefaction waves traveling in opposite directions is analogous to the shadow area, which refers to the very low density of the droplets around a solid surface after the droplets impinge on the solid surface.Although various numerical methods for preserving the positivity and avoiding the negative density have been published in the past, only a few methods satisfy the positivity in a vacuum state.In particular, the HLLC approximate Riemann solver developed by Toro (2001) and Toro et al. (1994) has shown good behavior and satisfies the positivity.In fact, the HLLC approximate Riemann solver for droplet equations can be archived in two-or three-dimensional cases because the HLLC solver basically requires at least three equations to complete wave structures reflecting left, right and shear waves.For this reason, a two-dimensional architecture for the semi-coupled algorithm is derived.
Figure 2 illustrates the assumed wave structure in the HLLC approximate Riemann solver for the semi-coupled algorithm.
(8) (9) (10) where: superscript k denotes each phase, such as air and droplet phases, and F * is an inter-cell flux in the star region for each phase, as defined by F *w = F w + S w (U *w -U w ) and a subscript w indicates left or right directions.The state, U *w , for each phase is given by: where: wave speeds, S L and S R (L and R subscripts denote the left and right states, respectively), for air flows, employ the direct wave speed estimation, as suggested by Davis (1988), which is the most well-known approach for estimating bounds for the minimum and maximum signal velocities in the solution of the Riemann problem: and, for droplet, Even though the minimum and maximum principals for wave speed estimations of left and right sides in Eq. ( 11) can be imposed to droplet equations, the principals may not ensure the positivity of density in a vacuum state.For this issue, the depth positivity condition proposed by Toro is applied to circumvent numerical difficulties via an assumption of density in a star region: In Eq. 12, q w is defined by: where: * is an estimate of the exact solution for ρ d in the star region.From the depth positivity condition, ρ * in the star region is derived as follows: The middle wave speeds (S * ) in the HLLC Riemann solver are based on the Rankine-Hugoniot conditions across each of the waves, which are basically assumed that all wave speeds are known.The speeds S * , purely in terms of the speeds S L and S R , defined in Eqs.11 and 12 for each phase, can be expressed as: for air,

two-dimensionAl finite volume formulAtion
The present finite volume formulation for the semicoupled algorithm is based on a cell-centered scheme and structured grid.The complete set for the two-dimensional semi-coupled algorithm can be written as: In the two-dimensional finite volume formulation, the line integral dl on the right-hand side can be approximated by a sum of the fluxes crossing the faces of the control volume.It is usually assumed that the flux is constant along the individual interface and is evaluated at the mid-point of the interface.Finally, the following discretized equation is: The value of droplet density at the interface in source terms is determined as follows: The air flow indeed should take into account viscous effects.In this study, the Spalart-Allmaras turbulence model is employed to simulate the turbulent effect in air flow fields and an ideal state of equation is used to close the system of equations.

boundAry conditions
A collection efficiency, a main feature in the in-flight and atmospheric icing, is determined by the droplet density and velocity near the solid surface.It is the fraction of the liquid water, which is deposited as ice on that component while flying in icing conditions.Therefore, the setting of the numerical boundary condition on the solid surface is an important factor in solving droplet flow fields.When the projection of a normal vector on a solid surface and the droplet velocity in an adjacent cell on the solid surface are positive, the droplets should not collide with the solid surface.On the contrary, when the projection is negative, the droplets should collide with the solid surface.On the basis of this observation, the following simple boundary conditions used in the previous study (Jung and Myong 2013) can be imposed for the droplet size less than 40 μm (see Fig. 3): In Eq. 19, the numerical flux H k at k-th interface is determined by the HLLC approximate Riemann solver.Although the implementation of explicit schemes is much easier than that of implicit schemes, explicit schemes require careful time step selection in order to fulfill the stability requirement.The semicoupled algorithm in this study is calculated under the maximum allowable time steps used in the work of Erduran et al. (2002).The local time stepping for the temporal discretization is achieved using the fifth stage Runge-Kutta scheme.In case of unsteady simulation, the minimum time step, after convergences for each phase, should be chosen to guarantee the stability for each phase.
where: n = (n x , n y ) denotes the normal vector of the solid surface.For the boundary conditions on the solid surface of air flows, a non-slip condition is applied.For the far-fields, Riemann invariant conditions for each phase are applied.

verIFIcatIon and valIdatIon
Two types of problems are selected to verify and validate the second-order positivity-preserving finite volume upwind ∫ scheme of the semi-coupled algorithm for air-mixed droplet flow.
The first problem is intended to compare the exact analytical and numerical solutions as a verification study.The second problem is the air-mixed droplet flow around an airfoil taken from Morency et al. (1999) as a validation study.The conditions, which represent the shadow (very low density of droplet) and high LWC (very high density of droplet) regions in Fig. 1, as the Riemann problems, for the first problem are summarized in Table 1.The exact and numerical solutions to the onedimensional well-posed hyperbolic part of the semi-coupled algorithm are compared in Figs. 4 and 5.The computational domain size 50 and 100 grid points with the CFL number 0.2 are used in these computations.Figure 4 (case 1) shows the positivity, which is most vulnerable to the negative density, between two identical rarefaction waves traveling in opposite directions.The HLLC scheme satisfies the density positivity condition of droplet; the density of droplet remains very low but always positive in the vacuum region.The new scheme is found to be very accurate for resolving rarefaction waves.The gap in the velocity of droplet is, nonetheless, found in the vicinity of the very low density of droplet.The numerical reason behind this shortcoming for the HLLC scheme is well understood (Toro 2001) and, under the more pressing need of the positivity property, the issue is left for future study.Another Riemann problem (case 2) and its numerical solutions are illustrated in Fig. 5.The case is considered in order to investigate the evolution of non-linear waves after the collision of two streams, which represents the local high LWC region around a leading edge of an airfoil where droplets of the free stream collide with the solid surface.It is clearly shown in Fig. 5 that the density of droplet increases drastically in the star region, indicating a locally high LWC region.The new scheme for this case is found to be very accurate in predicting the shock location and resolving the shock discontinuities.Prior to validate the new scheme for air-mixed droplet flow in the in-flight and atmospheric icing conditions, the grid sensitivity tests for the new scheme have been conducted by using the fine and coarse grids.Figures 6  and 7 show the grid topologies and convergence histories, respectively.The coarse grid requires less iteration than the fine grid; however, the accuracy due to grid topologies shows that the pressure coefficient of fine grid is more accurate than the coarse grid in Fig. 8. Nevertheless, the new scheme shows the robust characteristics with regard to the grid topologies.In order to validate the new scheme, experimental test cases are selected from the literature (Morency et al. 1999;Silveira et al. 2003;Vu et al. 2002).From the literature survey, the local collection efficiency is defined as the normalized influx of droplets at a given location, (Bourgault et al. 1999), where the velocity    are used.Figure 9 shows the structure mesh distribution around the NACA0012 airfoil and the mesh near the solid surface are clustered to capture viscous effects.Figure 10 shows the distributions of pressure coefficient and the collection efficiency on the solid surface for each phase, respectively.The present computational results of each phase are in close agreement with experimental data, demonstrating that the present semicoupled algorithm code is capable of producing the air and droplet information.AoA: angle of attack.
Also, the validation tests are extended and simulated via the existing experimental results, which are conducted by the NASA in USA, in order to compare the present results, including other existing codes, which are provided from the literature (Vu et al. 2002).In particular, the effects of higher angles of attack, different droplet sizes and airfoil shapes are presented and the test conditions are summarized in Table 2. Figure 13 shows the pressure coefficients on the surfaces of MS-0317 and NLF-0414 airfoils.In case of the LEWICE, the angles of attack in computations were slightly adjusted by approximately -2.0° to 2.5° depending on the test model, in order to match the experimental pressure.The proposed method, however, does not match the experimental pressure.Nevertheless, the present results are in relatively good agreement with experimental results.Figure 14 shows the effects of higher angles of attack of the MS-0317 airfoil with other fixed test conditions.The droplet impingement limits are moved from the upper surface to the lower surface around a leading edge, and maximum impingement intensities are slightly decreased in case of an angle of attack of 8°.Regarding the MVD effects, the maximum impingement intensity and droplet impingement limit are proportional to the size of droplet diameters.Figure 15 clearly shows that small droplet size is less effective than the large droplet size.Figures 16 and 17 present the airfoil shape effects with other fixed test conditions.In fact, the shape effects may be considered as an influence of the leading edge radius shown in Fig. 16, where the MS-0317 airfoil has a larger leading-edge

conclusIons
A semi-coupled algorithm was proposed and its twodimensional second-order positivity-preserving finite volume upwind solver was developed for the semi-coupled algorithm of air and droplet flows in atmospheric icing conditions.Much attention was paid to the mathematical descriptions of semi-coupled algorithm and preservation of density positivity in the numerical scheme, as well as the verification and validation study of the scheme.From the previous study (Jung and Myong 2013), the positivity-preserving property was achieved by the splitting of the original mathematical system into the well-posed hyperbolic part and the source term and then applying the positivity-preserving HLLC approximate Riemann solver to the semi-coupled algorithm.
In the verification and validation studies, it was shown that the new scheme can solve the air and droplet flow fields in fairly good agreement with the exact analytical solutions and experimental data, respectively.For the verification, the vacuum state and the shock wave due to two rarefaction waves and the collision of two waves were tested.In particular, a shock wave due to the collision of two waves could be considered as a droplet path line intersection, which represents the locally high LWC region around a leading edge of an airfoil where droplets of the free stream collide with the solid surface.It is clearly shown in Fig. 5 that the density of droplet increases drastically in the star region, indicating a locally high LWC region.The new scheme for this case is found to be very accurate in comparisons of exact solutions with the shock location and the shock discontinuities.In addition, Fig. 11 shows the locally high LWC distributions around the leading edge of the airfoil.From Figs. 5 and 11, the present scheme for the droplet path line intersection may circumvent the local accuracy limitations due to the numerical instability in the Eulerian framework.For the validation, various test models were chosen and the new scheme showed the applicability for air-mixed in-flight icing conditions in comparisons of existing wind tunnel droplet impingement test results.In particular, it accurately predicted the maximum value of the droplet impingement intensity near the stagnation region and the droplet impingement area.Such capability may be due to the second-order accuracy of numerical scheme in spatial discretization, while satisfying the positivity-preserving requirement.It should be noted that, without proper positivitypreserving schemes, to slow the appearance of breakdown due to negative density, the numerical code should remain firstorder instead of second-order, typical in CFD, which in turn may cause the underprediction of impingement intensity near the stagnation region.
In this study, a construction of the semi-coupled algorithm and its numerical approaches were mainly considered; nonetheless, there is still a room of improvement in the computational modeling of the moving body techniques for unsteady simulations.For example, periodically up and down motions of airfoils in order to simulate a section of helicopter rotor blades.The study of the problem will be the subject of a future paper.

Figure 1 .
Figure 1.Liquid water content distribution around an airfoil and identification of the Riemann problems (Jung and Myong 2013).
for Semi-Coupled Multiphase Flow in Atmospheric Icing Conditions

Figure 2 .
Figure 2. The HLLC approximate Riemann solver for the well-posed hyperbolic part of the semi-coupled algorithm.
In order to enhance the accuracy, Monotone Upstream Scheme for Conservation Law (MUSCL) proposed by van Leer (1979) together with van Albada limiter(van Albada et al. 1982) is employed:where: W = [ρ a , u a , ψ a , P a , ρ d , u d , ψ d ] T denotes primitive variables.The van Albada limiter is defined as Ψ L/R = 0.5[(1 + κ)r L/R + (1 -κ)]Ф L/R and Ф(r) = 2r/(1 + r 2 ).κ represents an extrapolation parameter and r R = (W i+1 -W i )/(W i+2 -W i+1 ) and r L = (W i+1 -W i )/(W i -W i-1).A general procedure to precisely solve the Riemann problem for two-dimensional and x-split air and droplet equations is referenced in Toro's work for air flows and previous work for the shallow water droplet model(Jung and Myong 2013), respectively.
represents the bounding surface of the control volume Ω.The two-dimensional finite volume formulation in general non-Cartesian domains can be derived by exploiting the rotational invariance property of the droplet equations, cosθF(U) + sinθG(U) = T -1 F(TU), where T = T(θ) represents the rotation matrix.Eqs.18 can then be expressed as:
The verification cases where x 0 and t represent the position of the initial discontinuity and the output time in seconds, respectively.Rarefaction (case 1) and shock (case 2) waves traveling in opposite directions.

Figure 4 .
Figure 4.The HLLC and the exact solutions of the density and velocity of each phase.Air (a and c); Droplet for semi-coupled algorithm -test case 1 (b and d).

Figure 5 .
Figure 5.The HLLC and the exact solutions of the density and velocity of each phase.Air (a and c); Droplet for semi-coupled algorithm -test case 2 (b and d).

Figure 7 .
Figure 7. Convergence graphs for fine (a) and coarse grid (b), and the convergence limit for each variable is equal to 10 -3 .rho: density.
Figure11shows the pressure and LWC distributions around the test model with an angle of attack 0°.In particular, a very low density of LWC is shown in the shadow area near the upper and lower side of the airfoil.The area of shadow zone at both the upper and lower sides of the airfoil is almost the same, mainly due to the effect of the free stream angle of attack and symmetrical shape.Furthermore, in order to check the low LWC results compared with free stream flow LWC, the low LWC values in the shadow area around the middle point of airfoil are presented in Fig.12.In the shadow areas, the low LWC values show the 1.4465e-05 g/m 3 and 7.3195e-06 g/m 3 around the upper and lower surfaces, respectively.Those values are less than 0.0002% compared with the free stream flow LWC.

Figure 8 .
Figure 8. Comparisons of fine and coarse grids with the experiment data.

Figure 10 .
Figure 10.Present and experimental results of pressure coefficient (a) and collection efficiency (b) for the NACA 0012.

Figure 12 .Figure 14 .Figure 15 .
Figure 12.The low LWC values around middle points of upper and lower sides of the NACA0012 airfoil in the shadow area.

table 2 .
Test models and conditions.