Numerical Analysis of Droplet Impacting on an Immiscible Liquid via Three-Phase Field Method

In this work, we establish a two-dimensional axisymmetric simulation model to numerically study the impacting behaviors between oil droplets and an immiscible aqueous solution based on the three-phase field method. The numerical model is established by using the commercial software of COMSOL Multiphysics first and then validated by comparing the numerical results with the previous experimental study. The simulation results show that under the impact of oil droplets, a crater will form on the surface of the aqueous solution, which firstly expands and then collapses with the transfer and dissipation of kinetic energy of this three-phase system. As for the droplet, it flattens, spreads, stretches, or immerses on the crater surface and finally achieves an equilibrium state at the gas–liquid interface after experiencing several sinking-bouncing circles. The impacting velocity, fluid density, viscosity, interfacial tension, droplet size, and the property of non-Newtonian fluids all play important roles in the impact between oil droplets and aqueous solution. The conclusions can help to cognize the mechanism of droplet impact on an immiscible fluid and provide useful guidelines for those applications concerning droplet impact.

As for the impacting hydrodynamics of droplets, it has accumulated the interest of researchers for many years. Now, the impacting behaviors of droplets on various receiving surfaces ranging from solid surfaces to fluidic interfaces have been studied [24][25][26]. Recently, with the advancements of high-speed imaging techniques (HSIs) and particle image velocimetry (PIV), the impacting details of droplet can be finer captured and observed [27]. Although these studies are interesting and meaningful for understanding the droplet-impacting mechanism, most of the present analyses consider a particular class of drop impact where the droplet is made of the same fluid as one of the other two phases, which essentially sets up a biphasic drop impact problem [28,29]. As for the impact between

Problem Description
Oil-in-air droplet impacting with an aqueous solution refers to three-phase fluids. To numerically analyze the droplet-impacting behaviors, a two-dimensional simulation model based on the three-phase field method is established by using the commercial software of COMSOL Multiphysics. The axisymmetric geometrical model with triangular meshes shown in Figure 1 is adopted. The upper computational domain with a depth of H 1 is the air phase (i = a). The lower computational domain with a height of H 3 is the aqueous solution phase (i = w). As for the oil droplet phase (i = o), its initial release height is H 2 . The width of the computational domain is 0.5 W. The specific values are H 1 = 20 mm, H 2 = 10 mm, H 3 = 20 mm, and W = 40 mm, respectively. The droplet diameter D ranges from 4 mm to 12 mm in this study. To mimic the freefall of a droplet from a long distance without increasing the unnecessary computational time, the droplet has an initial downward velocity of u and then migrates under the gravity effect.

Numerical Method
The impact of oil-in-air droplets on an aqueous solution is a problem concerning three immiscible fluids, and accurate phase interface tracking is a key problem during simulation. Up to now, there are many methods for tracking fluidic interfaces, such as Lattice Boltzmann, level set, phase field, and volume of fluids [40][41][42]. Among these approaches, the phase field method is attractive because of its better mass conservation, rapid computation speed, and satisfactory algorithm stability [43]. Therefore, in this study, we adopt the three-phase field method to establish a two-dimensional simulation model. To track the fluidic interfaces of three immiscible fluids, the following Cahn-Hilliard equations are solved [44,45]: (1) (2) (3) (4) (5) where is the phase field variable, u is the velocity vector, is the molecular mobility parameter, is the generalized chemical potential, is the control parameter of interface thickness, σi,j is the interfacial tension, and is the free energy of the three-phase system.
The phase field variables vary between 0 and 1 and are a measure of the concentration of each phase. At each point, the phase field variables satisfy the following equation: (6) To ascertain continuous variation across the fluidic interface, the density and dynamic viscosity of the fluid mixture are defined as:

Numerical Method
The impact of oil-in-air droplets on an aqueous solution is a problem concerning three immiscible fluids, and accurate phase interface tracking is a key problem during simulation. Up to now, there are many methods for tracking fluidic interfaces, such as Lattice Boltzmann, level set, phase field, and volume of fluids [40][41][42]. Among these approaches, the phase field method is attractive because of its better mass conservation, rapid computation speed, and satisfactory algorithm stability [43]. Therefore, in this study, we adopt the three-phase field method to establish a two-dimensional simulation model. To track the fluidic interfaces of three immiscible fluids, the following Cahn-Hilliard equations are solved [44,45]: where ϕ i is the phase field variable, u is the velocity vector, M ϕ is the molecular mobility parameter, η i is the generalized chemical potential, is the control parameter of interface thickness, Σ i,j is the interfacial tension, and F(φ) is the free energy of the threephase system. The phase field variables vary between 0 and 1 and are a measure of the concentration of each phase. At each point, the phase field variables satisfy the following equation: To ascertain continuous variation across the fluidic interface, the density and dynamic viscosity of the fluid mixture are defined as: where ρ w , ρ o , and ρ a are the densities of the aqueous solution, oil, and air, respectively, µ w , µ o and µ a are the dynamic viscosities of these three phases, respectively. In our study, the material of the droplet is silicone oil, and the aqueous solution is a mixture of 60 wt% glycerol, 38 wt% deionized water, and 2 wt% polyvinyl alcohol. Their basic physical properties of fluids are listed in Table 1. Table 1. Physical properties of fluids used in this study [46][47][48]. The fluids also meet the Navier-Stokes and the continuity equations:

Fluid
where ρ is the volume averaged density, p is the pressure, µ is the dynamic viscosity of the fluid, g is the gravitational acceleration, and F st is the capillary force per unit volume. The capillary force is described as follows:

Initial and Boundary Conditions
The fluidic phases are initially considered to be quiescent and at uniform pressure, u = 0, and p = constant at t = 0 The droplet with an initial velocity of u migrates under gravity. These conditions are the initial conditions employed for the study, whereas the following boundary conditions are employed in order to solve the governing equations.
At the sidewall, the free-slip boundary condition is set, where n is the normal vector. At the bottom wall, we have employed the following no-slip and impermeability conditions: A pressure point constraint is employed to enforce a zero reference gauge pressure at the central point of the air-aqueous solution interface.
The boundary conditions for Equation (1) are, in general, homogeneous Neumann type. For the chemical potential η i , this condition ensures that there is no mass diffusion through the boundary, For the parameter (ϕ i ), wetted wall boundary conditions are incorporated on the side and bottom walls [32], The contact angle of the interface in between phases 'i' and 'j' is denoted by θ ij .

Solution Methodology and Validation
A commercial software of COMSOL Multiphysics 6.0 is used to solve the aforementioned time-dependent partial differential equations. The variables for flow and phase fields are calculated by using PARDISO solvers. In the simulation, the environmental temperature is fixed at 293.15 K. The droplet, with an initial downward velocity, migrates under the effect of gravity and finally impacts the aqueous solution. Before formal numerical analyses, we conduct a grid dependence test and compare our results with the experiments conducted by Jain et al. [30] to demonstrate the accuracy of the simulation model, as shown in Section S1 of the Supplementary Materials. The suitable agreement between simulations and experiments indicates the accuracy and feasibility of this numerical model. What is more, the simulation model can keep the conservation of mass at different time points, as shown in Figure 2. Specifically, when the impacting velocity of the droplet is 0.64 m/s, the impacting states between the droplet and liquid film under different time points are given in Figure 2a (see Video S1). It can be found that the droplet migrates downward first and is bounced back afterward; finally, it reaches the balanced state at the air-water interface because of the transfer and dissipation of kinetic energy (see Section S2 of the Supplementary Materials). During the impacting process, the volume fractions of the three phases are listed in Figure 2b. Clearly, the mass of the three phases is conservative at different time points, which is also evidence of the feasibility of our numerical model. Thus, the simulation results can help us to understand the droplet-impacting hydrodynamics.
The boundary conditions for Equation (1) are, in general, homogeneous Neumann type. For the chemical potential , this condition ensures that there is no mass diffusion through the boundary, (16) For the parameter ( ), wetted wall boundary conditions are incorporated on the side and bottom walls [32], (17) The contact angle of the interface in between phases 'i' and 'j' is denoted by .

Solution Methodology and Validation
A commercial software of COMSOL Multiphysics 6.0 is used to solve the aforementioned time-dependent partial differential equations. The variables for flow and phase fields are calculated by using PARDISO solvers. In the simulation, the environmental temperature is fixed at 293.15 K. The droplet, with an initial downward velocity, migrates under the effect of gravity and finally impacts the aqueous solution. Before formal numerical analyses, we conduct a grid dependence test and compare our results with the experiments conducted by Jain et al. [30] to demonstrate the accuracy of the simulation model, as shown in Section S1 of the Supplementary Materials. The suitable agreement between simulations and experiments indicates the accuracy and feasibility of this numerical model. What is more, the simulation model can keep the conservation of mass at different time points, as shown in Figure 2. Specifically, when the impacting velocity of the droplet is 0.64 m/s, the impacting states between the droplet and liquid film under different time points are given in Figure 2a (see Video S1). It can be found that the droplet migrates downward first and is bounced back afterward; finally, it reaches the balanced state at the air-water interface because of the transfer and dissipation of kinetic energy (see Section S2 of the Supplementary Materials). During the impacting process, the volume fractions of the three phases are listed in Figure 2b. Clearly, the mass of the three phases is conservative at different time points, which is also evidence of the feasibility of our numerical model. Thus, the simulation results can help us to understand the droplet-impacting hydrodynamics.

Results and Discussions
The impacting behaviors between an oil-in-air droplet and an aqueous solutions are influenced by droplet kinetic energy, inertial (gravity) force, interfacial tensions, and viscous force. The magnitudes of these forces are affected by various factors, such as fluid viscosity, density, interfacial tensions, impacting velocity, droplet diameter, and Newtonian and non-Newtonian properties of fluids. To understand their effects on

Results and Discussions
The impacting behaviors between an oil-in-air droplet and an aqueous solutions are influenced by droplet kinetic energy, inertial (gravity) force, interfacial tensions, and viscous force. The magnitudes of these forces are affected by various factors, such as fluid viscosity, density, interfacial tensions, impacting velocity, droplet diameter, and Newtonian and non-Newtonian properties of fluids. To understand their effects on droplet impacting, we will conduct systematic numerical analyses in the following sections.

Influence of Fluid Viscosity on Droplet Impacting
In the droplet (phase o)-aqueous solution (phase w)-air (phase a) system, to analyze the influences of the dynamic viscosities of oil (µ o ) and aqueous solution (µ w ) on droplet impacting, the initial downward velocity and diameter of the droplet are fixed at 1 m/s and 4 mm, respectively. Moreover, the impacting velocity will reach V = 1.08 m/s at the moment of impacting under the effect of gravity. Figure 2a indicates that the impacted droplet would migrate downward first and bounce back afterward; finally, it achieves the balanced state at the air-water interface. So, we use the maximum sinking states of the droplet to judge the influences of fluid viscosity on the impacting behaviors between the oil droplet and aqueous solution. The maximum sinking states of oil droplets after impacting with an aqueous solution under different viscosities of aqueous solution are exhibited in Figure 3a (also see Video S2). Clearly, with the increase in aqueous solution viscosity, the sinking depth of the droplet and the size of the water crater both decrease. The configurations of maximum craters under different viscosities of aqueous solution are shown in Figure 3b. Then, we can obtain the quantitative curves of crater sizes versus aqueous solution viscosity, as shown in Figure 3c. It can be found that the crater depth d 1M decreases from 18.61 mm to 6.79 mm as the viscosity of the aqueous solution rises from 0.01 Pa·s to 10 Pa·s. As for the width of crater w, it decreases by 93.8% when the viscosity of the aqueous solution is improved from 0.01 Pa·s to 10 Pa·s. This phenomenon is due to the fact that the viscosity of the aqueous solution has a negative effect on fluid flow, resulting in the impact kinetic energy between the droplet and liquid film having a weak effect on droplet sinking and liquid crater formation as the viscosity of the aqueous solution rises.

Influence of Fluid Density on Droplet Impacting
In the above analyses, the densities of oil and aqueous solution are fixed at 970 kg/m 3 and 1150 kg/m 3 , respectively. Obviously, the density ratio of oil to water δ = ρo/ρw will affect the gravity and buoyancy force acting on the droplet. To understand the effect of fluid density on droplet impacting, we fix the density of the aqueous solution at 1150 kg/m 3 and then adjust the density ratio δ to conduct simulations of droplet impacting. When the impacting velocity and density ratio are 0.64 m/s and 1.4, respectively, the impacting states between the oil-in-air droplet and aqueous solution under different time points are shown in Figure 4a (Video S3). After impacting, the droplet first migrates downward until the water crater reaches the maximum size (the impacting state at t = 0.026 s). Then, the aqueous solution begins to bounce back, since the upward force of the aqueous solution applied on the droplet is smaller than the gravity of the droplet; the The maximum sinking states of the oil droplet after impacting with an aqueous solution under different oil viscosities are given in Figure 3d (see Video S2). We can find that the maximum sinking depth of the droplet rises with the rising of oil viscosity. The configurations of the formed liquid crater at this state are shown in Figure 3e. Accordingly, the statistical depth d 1M and width w of the liquid crater under different oil viscosities are listed in Figure 3f. It is clear that the crater depth d 1M rises from 4.95 mm to 6.64 mm as the oil viscosity rises from 0.01 Pa·s to 10 Pa·s. In contrast, the crater width w is negatively affected by the increased oil viscosity. Specifically, when the oil viscosity increases from 0.01 Pa·s to 10 Pa·s, the w decreases from 19.57 mm to 17.00 mm. This can be explained by the following statements. With the rising of oil viscosity, the deformability of the oil droplet in the horizontal direction decreases during the impacting process, resulting in the impact kinetic energy becoming more concentrated at the impact position. Obviously, the more concentrated impact kinetic energy will lead to an increase in liquid crater depth and a decrease in crater width.

Influence of Fluid Density on Droplet Impacting
In the above analyses, the densities of oil and aqueous solution are fixed at 970 kg/m 3 and 1150 kg/m 3 , respectively. Obviously, the density ratio of oil to water δ = ρ o /ρ w will affect the gravity and buoyancy force acting on the droplet. To understand the effect of fluid density on droplet impacting, we fix the density of the aqueous solution at 1150 kg/m 3 and then adjust the density ratio δ to conduct simulations of droplet impacting. When the impacting velocity and density ratio are 0.64 m/s and 1.4, respectively, the impacting states between the oil-in-air droplet and aqueous solution under different time points are shown in Figure 4a (Video S3). After impacting, the droplet first migrates downward until the water crater reaches the maximum size (the impacting state at t = 0.026 s). Then, the aqueous solution begins to bounce back, since the upward force of the aqueous solution applied on the droplet is smaller than the gravity of the droplet; the partial oil droplet will pass through the gas-liquid interface and enter the aqueous solution, and the residual oil droplet will be located at the gas-liquid interface. However, at δ = 0.6, the total oil droplet will be bounced back because the upward force of the aqueous solution applied on the droplet is higher than the gravity of the droplet, as shown in Figure 4b (Video S3).

Influence of Interfacial Tension on Droplet Impacting
The configurations of fluidic interfaces are deeply affected by the interfacial tensions of fluids. In order to understand the influence of interfacial tension on droplet-impacting behaviors, we further simulate the impacting hydrodynamics of three phases. The above studies have demonstrated that the droplet would be balanced at the air-liquid interface at the final state. At this stage, the interfacial tension forces acting on the droplet are shown in Figure 5a, and these three interfacial tensions meet Young's relation. According to Young's relation, we can conclude that the relative magnitudes of these three interfacial tensions are related to the droplet state at the air-liquid interface.
When conducting simulations, the impacting velocity and droplet diameter are fixed at 0.64 m/s and 4 mm, respectively. The influences of σo-a on droplet impact and the size of the formed liquid crater are first analyzed, as shown in Figure 5b (see Video S4). It can be seen that the maximum and final depths of the liquid crater are both positively affected by the increased σo-a. Concretely, the maximum crater depth d1M increases from 3.90 mm to 5.65 mm when the σo-a rises from 0.02 N/m to 0.06 N/m. The final crater depth d1F at the stable state increases by 86.2% as the σo-a is improved from 0.02 N/m to 0.06 N/m. This is due to the fact that the increased σo-a makes the improvement of the interfacial free energy, so to maintain the stability of the three-phase fluidic system, more surfaces of the oil droplet will come into contact with the aqueous solution. That is to say, The curve of the maximum liquid crater depth d 1M versus density ratio is given in Figure 4c. The crater depth d 1M is positively affected by the density ratio δ. Specifically, d 1M rises from 2.8 mm to 7.0 mm when δ is increased from 0.6 to 1.8. During the change of density ratio, the oil droplets in some conditions will be completely bounced back, whereas the partial oil droplet enters an aqueous solution, and the residual one bounces back when the oil density is much higher than that of water. To quantify this influence, we further analyze the ratios of the volume of the droplet entering the aqueous solution to the total volume of the oil droplet under different fluid density ratios, as shown in Figure 4d. It can be found that the bouncing state of the droplet varies when the density ratio δ is higher than 1.3. The volume ratio increases from 0 to 0.83 as δ rises from 0.6 to 1.8.

Influence of Interfacial Tension on Droplet Impacting
The configurations of fluidic interfaces are deeply affected by the interfacial tensions of fluids. In order to understand the influence of interfacial tension on droplet-impacting behaviors, we further simulate the impacting hydrodynamics of three phases. The above studies have demonstrated that the droplet would be balanced at the air-liquid interface at the final state. At this stage, the interfacial tension forces acting on the droplet are shown in Figure 5a, and these three interfacial tensions meet Young's relation. According to Young's relation, we can conclude that the relative magnitudes of these three interfacial tensions are related to the droplet state at the air-liquid interface.

Influence of Impacting Velocity on Droplet Impacting
The impacting kinetic energy is another factor affecting the impacting behaviors between the oil droplet and the aqueous solution. So, we have set different initial downward velocities for the droplet to conduct simulations. When the impacting velocity and droplet diameter are fixed at 2.08 m/s and 4 mm, respectively, the impacting states of the droplet under different time points are shown in Figure 6a (see Video S5). It is clear that the droplet will migrate downward first and be bounced back afterward; finally, it reaches the balanced state at the air-water interface. During this process, the droplet will have apparent sinking and bouncing distances (d1M and d2M). If the impacting velocity is decreased to 1.08 m/s, the sinking and bouncing distances of the droplet and the liquid crater become much smaller, as shown in Figure 6b (Video S5).
To quantify the impacting behaviors of the droplet under different impacting velocities, we have plotted the droplet displacement under different time points in Figure  6c. It is clear that the droplet will experience a sinking-bouncing circle before reaching the equilibrium state at V = 1.08 mm/s. However, when the impacting velocity is increased to 2.08 m/s, the droplet will undergo two sinking-bouncing circles before reaching the equilibrium state, and the sinking and bouncing distances in the second circle are smaller than that in the first circle due to the dissipation of energy. The configurations of the airliquid interfaces when the crater is in the deepest and highest positions are shown in Figure 6d. Obviously, the impacting velocity has a significant influence on droplet When conducting simulations, the impacting velocity and droplet diameter are fixed at 0.64 m/s and 4 mm, respectively. The influences of σ o-a on droplet impact and the size of the formed liquid crater are first analyzed, as shown in Figure 5b (see Video S4). It can be seen that the maximum and final depths of the liquid crater are both positively affected by the increased σ o-a . Concretely, the maximum crater depth d 1M increases from 3.90 mm to 5.65 mm when the σ o-a rises from 0.02 N/m to 0.06 N/m. The final crater depth d 1F at the stable state increases by 86.2% as the σ o-a is improved from 0.02 N/m to 0.06 N/m. This is due to the fact that the increased σ o-a makes the improvement of the interfacial free energy, so to maintain the stability of the three-phase fluidic system, more surfaces of the oil droplet will come into contact with the aqueous solution. That is to say, the droplet will have a deeper sinking state, and the crater depth will be higher.
When the impacting velocity, σ w-a , and σ w-o are fixed at 0.64 m/s, 0.050 N/m, 0.020 N/m, respectively, we further obtain a series of droplet-impacting states under different σ w-o , as exhibited in Figure 5c (also see Video S4). It can be seen that the maximum crater depth d 1M rises from 3.90 mm to 4.67 mm as the σ w-o increases from 0.035 N/m to 0.060 N/m. As for the final crater depth d 1F , it increases by 11.3% when the σ w-o changes from 0.035 N/m to 0.060 N/m. This variation trend can be explained by the following statements. With the rising of σ w-o , in order to decrease the interfacial free energy of the fluidic system, the droplet becomes more repellent to contact with an aqueous solution. That is, the contact area between the droplet and the aqueous solution decreases with the increased σ w-o . Accordingly, the decreased contact area leads to increased contact stress during impact, resulting in an increase in droplet sinking and liquid crater depth.
Furthermore, the effect of σ w-a on the impacting behaviors between the oil droplet and the aqueous solution is investigated, as illustrated in Figure 5d (see Video S4). With the increase in the liquid-air interfacial tension σ w-a , the energy barrier for droplets crossing the gas-liquid interface also increases. Thus, the maximum and final depths of the liquid crater are both negatively affected by the increased σ w-a . Specifically, the maximum crater depth d 1M reduces from 5.68 mm to 3.90 mm as σ w-a rises from 0.02 N/m to 0.05 N/m. The final crater depth at the stable state decreases by 21% when the σ w-a is increased from 0.02 N/m to 0.05 N/m.

Influence of Impacting Velocity on Droplet Impacting
The impacting kinetic energy is another factor affecting the impacting behaviors between the oil droplet and the aqueous solution. So, we have set different initial downward velocities for the droplet to conduct simulations. When the impacting velocity and droplet diameter are fixed at 2.08 m/s and 4 mm, respectively, the impacting states of the droplet under different time points are shown in Figure 6a (see Video S5). It is clear that the droplet will migrate downward first and be bounced back afterward; finally, it reaches the balanced state at the air-water interface. During this process, the droplet will have apparent sinking and bouncing distances (d 1M and d 2M ). If the impacting velocity is decreased to 1.08 m/s, the sinking and bouncing distances of the droplet and the liquid crater become much smaller, as shown in Figure 6b (Video S5). impact. That is, the crater depth and bounced distance of the droplet were both positively affected by the increased impacting velocity. Specifically, the maximum crater depth d1M rises from 1.80 mm to 6.80 mm as the impacting velocity increases from 0.1 m/s to 2 m/s (Figure 6e). Moreover, the maximum bouncing height of the droplet d2M increases from 0.02 mm to 3.37 mm as the impacting velocity is improved from 0.1 m/s to 2 m/s ( Figure  6e).

Influence of Droplet Diameter on Impacting States
When two droplets with different sizes have the same velocity, it is obvious that the bigger droplet has the larger kinetic energy compared with that of the smaller one. So the droplet diameter will also influence the impacting behaviors between the oil droplet and the aqueous solution. To understand the specific effect of droplet diameter on droplet To quantify the impacting behaviors of the droplet under different impacting velocities, we have plotted the droplet displacement under different time points in Figure 6c. It is clear that the droplet will experience a sinking-bouncing circle before reaching the equilibrium state at V = 1.08 mm/s. However, when the impacting velocity is increased to 2.08 m/s, the droplet will undergo two sinking-bouncing circles before reaching the equilibrium state, and the sinking and bouncing distances in the second circle are smaller than that in the first circle due to the dissipation of energy. The configurations of the air-liquid interfaces when the crater is in the deepest and highest positions are shown in Figure 6d. Obviously, the impacting velocity has a significant influence on droplet impact. That is, the crater depth and bounced distance of the droplet were both positively affected by the increased impacting velocity. Specifically, the maximum crater depth d 1M rises from 1.80 mm to 6.80 mm as the impacting velocity increases from 0.1 m/s to 2 m/s (Figure 6e). Moreover, the maximum bouncing height of the droplet d 2M increases from 0.02 mm to 3.37 mm as the impacting velocity is improved from 0.1 m/s to 2 m/s (Figure 6e).

Influence of Droplet Diameter on Impacting States
When two droplets with different sizes have the same velocity, it is obvious that the bigger droplet has the larger kinetic energy compared with that of the smaller one. So the droplet diameter will also influence the impacting behaviors between the oil droplet and the aqueous solution. To understand the specific effect of droplet diameter on droplet impact, the impacting states of droplets with different sizes are simulated, as shown in Figure 7a,b (see Video S6). Apparently, when the impacting velocity is fixed at 0.64 m/s, a bigger liquid crater will form when a larger oil droplet impacts with the aqueous solution. What is more, the droplets with different sizes will experience absolutely different sinkingbouncing circles. According to the curves of droplet displacement versus time plotted in Figure 7c, we can find that the droplet will undergo a sinking-bouncing circle before reaching a stable state when its diameter is 3 mm. In contrast, the drop with a diameter of 8 mm experiences several sinking-bouncing circles before reaching the equilibrium state. This is due to the fact that the droplet with the bigger diameter has a higher kinetic energy. The gradually damped sinking and bouncing of the droplet are related to the transfer and dissipation of kinetic energy.

Influence of Non-Newtonian on Droplet Impacting
In the above analyses, the air, oil, and aqueous solution are all Newtonian fluids; that is, their viscosities keep constant under different shear rates of fluids. Non-Newtonian fluids are common in our daily life, such as water and silicone oil. In To quantify the influence of droplet diameter on the impacting behaviors, we have plotted the air-liquid configurations under different droplet diameters when the droplet sinks to the deepest position, as shown in Figure 7d. It is apparent that the maximum crater size rises with the droplet diameter. Specifically, the maximum crater depth d 1M increases from 2.70 mm to 5.57 mm as the droplet diameter increases from 3 mm to 8 mm ( Figure 7e). As for the crater width w, it improves from 11.17 mm to 35.31 mm when the droplet diameter is increased from 3 mm to 8 mm (Figure 7e).

Influence of Non-Newtonian on Droplet Impacting
In the above analyses, the air, oil, and aqueous solution are all Newtonian fluids; that is, their viscosities keep constant under different shear rates of fluids. Non-Newtonian fluids are common in our daily life, such as water and silicone oil. In contrast, there are also lots of fluids whose viscosities are shear-rate-dependent in the fields of biology, chemistry, and material synthesis, such as human blood and starch solution. After analyzing the Newtonian fluids, we try to further study the droplet behaviors in non-Newtonian fluids. According to the relationship between viscosity and shear rate, non-Newtonian fluids can be classified as many kinds, such as the power-law fluid, Carreau fluid, Bingham fluid, and Casson fluid [49][50][51]. In these non-Newtonian fluids, the power-law fluid is a simple yet typical kind. Many studies concerning non-Newtonian fluids have taken it as an example [52,53]. So, in our study, the viscosities of non-Newtonian oil and aqueous solution also meet the power-law model.
The dynamic viscosity of the power-law fluid is defined by: where K is the consistency coefficient, and its value is fixed at 0.01 Pa·s n in our simulation. n is the power-law index.
. γ is the rate of the strain tensor, and it can be expressed by the following equation: The value of power-law index n has a significant influence on fluid property:      0 < n < 1 shear thinning fluid n = 1 Newtonian fluid n > 1 shear thickening fluid (20) Therefore, we have systematically analyzed the effect of the power-law index on the droplet impact in this section. The influence of the power-law index of aqueous solution n w on droplet impact is first studied. When conducting simulations, the impacting velocity of the droplet is fixed at 0.5 m/s, and the viscosity of the oil is set to be 0.01 Pa·s under different shear rates. The droplet-impacting states under two different n w are portrayed in Figure 8a,b, respectively (also see Video S7). When n w is 0.7, the aqueous solution is a shear-thinning fluid; that is, the shear rate will decrease the viscosity of water. So a large water crater will form under the impact of the droplet (Figure 8a). In contrast, the aqueous solution becomes shear thickening fluid at n w = 2.0; that is, the shear rate will increase fluid viscosity. So the crater size is much smaller compared with that at n w = 0.7 (Figure 8b). The plot of maximum crater depth d 1M versus the power-law index n w of water is exhibited in Figure 8c. The d 1M decreases from 3.55 mm to 1.46 mm as n w rises from 0.7 to 2. As for the power-law index n o of the oil phase, its effect on crater depth is shown in Figure 8d. It can be seen that the crater depths rise from 3.22 mm to 3.39 mm as n o is improved from 0.7 to 2.0. Obviously, this increasing trend is very small. This is due to the fact that the increased oil viscosity will improve water crater depth (Figure 3f), but n o ranging from 0.7 to 2.0 can only slightly increase the viscosity of the oil phase during impacting. nw of water is exhibited in Figure 8c. The d1M decreases from 3.55 mm to 1.46 mm as nw rises from 0.7 to 2. As for the power-law index no of the oil phase, its effect on crater depth is shown in Figure 8d. It can be seen that the crater depths rise from 3.22 mm to 3.39 mm as no is improved from 0.7 to 2.0. Obviously, this increasing trend is very small. This is due to the fact that the increased oil viscosity will improve water crater depth (Figure 3f), but no ranging from 0.7 to 2.0 can only slightly increase the viscosity of the oil phase during impacting.

Conclusions
In this article, to understand the physics of droplet impacting on another immiscible fluid, we establish a two-dimensional axisymmetric simulation model based on the three-phase field method to conduct numerical analyses. The accuracy of the established numerical model is first validated by the comparison between simulations and experiments. When carrying out simulations, the oil droplet with an initial downward velocity is released into the air phase first. Then, it migrates down under the gravity force, and finally, it impacts the aqueous solution. Under the impact of the oil droplet, a crater will form on the surface of the aqueous solution, which first expands and then collapses with the transfer and dissipation of kinetic energy of this three-phase system. As for the droplet, it flattens, spreads, stretches, or immerses on the crater surface and finally achieves an equilibrium state at the gas-liquid interface after experiencing several

Conclusions
In this article, to understand the physics of droplet impacting on another immiscible fluid, we establish a two-dimensional axisymmetric simulation model based on the three-phase field method to conduct numerical analyses. The accuracy of the established numerical model is first validated by the comparison between simulations and experiments. When carrying out simulations, the oil droplet with an initial downward velocity is released into the air phase first. Then, it migrates down under the gravity force, and finally, it impacts the aqueous solution. Under the impact of the oil droplet, a crater will form on the surface of the aqueous solution, which first expands and then collapses with the transfer and dissipation of kinetic energy of this three-phase system. As for the droplet, it flattens, spreads, stretches, or immerses on the crater surface and finally achieves an equilibrium state at the gas-liquid interface after experiencing several sinking-bouncing circles. The simulations indicate that the increased oil viscosity, density, interfacial tension σ o-a , σ w-o , impacting velocity, and droplet diameter will enhance the droplet penetration depth, water bouncing height, maximum water crater size, and the number of sinking-bouncing circles. In contrast, the impacting effects between droplet and liquid are negatively affected by the improved water viscosity, density, interfacial tension σ w-a , and power-law index n w . For instance, when the impacting velocity and dynamic viscosity of the aqueous solution are fixed at 1.08 m/s and 0.01 Pa·s, respectively, the maximum crater depth during impacting increases from 4.95 mm to 6.64 mm as the oil viscosity rises from 0.01 Pa·s to 10 Pa·s. These conclusions can help to know about the mechanism of droplet impacting on an immiscible fluid and provide useful guidelines for those applications concerning droplet impacting, such as deep-water oil spills and inkjet printing. Moreover, in the future, we will conduct experiments to analyze the impacting behaviors between oil droplets and aqueous solutions.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/mi14050951/s1, Figure S1: Grid dependence test at V 0 = 0.64 m/s, D = 4 mm and t = 0.1 s; Figure S2:Comparison between simulation and experiments; Video S1: The process of a liquid-in-air droplet impacting on an-other immiscible solution; Video S2: The influence of fluid viscosity on droplet impacting; Video S3: The influence of fluid density on droplet impacting; Video S4: The influence of fluidic interfacial tensions on droplet impacting; Video S5: The influence of impacting velocity on droplet impacting; Video S6: The influence of droplet diameter on droplet impacting; Video S7: The influence of non-Newtonian properties of fluids on droplet impacting.