Numerical Study on the Formation and Solidification of LMPA Microdroplet in a Microfluidic Device

LMPA droplets or particles have contributed to many fields such as the application of sensors and valves, and droplet-based microfluidics has been applied to the preparation of LMPA droplets. Understanding how flow rate, interfacial tension, and temperature affect the formation and solidification of droplets is helpful to design a microfluidic platform. In this study, a coupled VOF and enthalpy-porosity method will be used to numerically simulate how these factors affect the LMPA droplet formation and solidification process. We find that increasing the velocity of the continuous phase or decreasing the interfacial tension will reduce the LMPA droplet size and simultaneously increase the frequency of droplet formation. In addition, increasing the interfacial tension will decrease the required solidification time of LMPA droplets, and the solidification time of droplets will first increase and then decrease with the growth of continuous-phase velocity. On the other hand, increasing the continuous-phase temperature or cooling wall temperature will reduce the solidification time of LMPA droplets, but has no obvious influence on the size and frequency of droplet generation.

LMPA droplets or particles have contributed to many fields such as the application of sensors and valves, and droplet-based microfluidics has been applied to the preparation of LMPA droplets. Understanding how flow rate, interfacial tension, and temperature affect the formation and solidification of droplets is helpful to design a microfluidic platform. In this study, a coupled VOF and enthalpy-porosity method will be used to numerically simulate how these factors affect the LMPA droplet formation and solidification process. We find that increasing the velocity of the continuous phase or decreasing the interfacial tension will reduce the LMPA droplet size and simultaneously increase the frequency of droplet formation. In addition, increasing the interfacial tension will decrease the required solidification time of LMPA droplets, and the solidification time of droplets will first increase and then decrease with the growth of continuous-phase velocity. On the other hand, increasing the continuous-phase temperature or cooling wall temperature will reduce the solidification time of LMPA droplets, but has no obvious influence on the size and frequency of droplet generation.

INTRODUCTION
In recent years, low-melting-point alloy (LMPA) particles have been widely used in many fields because of their unique advantages. For example, micro LMPA particles are useful for conductive inks (Boley et al., 2015;Cano-Raya et al., 2019), sensors (Won et al., 2015;Tang et al., 2016), and composite materials (Mohammed et al., 2014;Bartlett et al., 2016). Applications such as mixers (Shamberger et al., 2017), pumps , and valves (Sander and Weigand, 2009) usually only need a small amount of LMPA particles, but the particle diameter needs to be accurately controlled. One method of controlling the particle size is to directly produce monomer or liquid droplets approximately equal to the desired particle size, and then to polymerize or solidify the droplets. Typical methods of forming micro-and nano-scale droplets include on-chip (Tang et al., 2016) or off-chip sonication treatment (Zhang et al., 2014) or molding (Mohammed et al., 2014). However, it is difficult to continuously produce micron-sized droplets in molding, and sonication cannot accurately control the droplet size. T-junction or flow-focusing microfluidic systems are able to continuously generate micro-scale droplets with a uniform diameter (Stone et al., 1986;Fu et al., 2012;Jamalabadi et al., 2017;Li et al., 2019a;Outokesh et al., 2022). Thus, in recent years, several experimental and numerical methods have been used to study and evaluate the droplet formation process and its influencing process    based on these devices (Li et al., 2019a;Cui et al., 2019;Lei et al., 2020;Outokesh et al., 2022). Computational fluid dynamics (CFD) simulation provides a method for an indepth understanding of this complicated process. Typically, the level-set method (LS) (Yap and Chai, 2012;Liu et al., 2015;Shin et al., 2018;Helland et al., 2019;Valle et al., 2020), volume of fluid (VOF) method (Soh et al., 2016;Chaves et al., 2020;Sattari et al., 2020;Garoosi and Hooman, 2022), and lattice Boltzmann method (Kim and Pitsch, 2015;Reijers et al., 2016;Wu et al., 2021;Li et al., 2022) are used to model the multiphase flow. Among them, the VOF technique is the most commonly employed technique to simulate the free surface in academic and commercial CFD software packages (Hirt and Nichols, 1981). One of the advantages of the VOF method is volume conservation because for a specific phase in each unit, it is the volume fraction which is calculated and tracked in this model instead of the interface. Based on the VOF technique, Hoang et al. (2013) carefully investigated the splitting processes of droplets in microchannels and the results showed that the rapid pinching of droplets was driven by the surface tension, which was consistent with previous experiments (Stone et al., 1986). Nekouei and Vanapalli (2017) extended the study of droplet formation with a high viscosity ratio using the VOF method and found that the squeezing time for droplet generation increased with the increase in the high viscosity ratio. Li et al. (2019b) used the VOF method to numerically address the droplet-generation process in ordinary and improved T-junction structures, and they successfully decreased the droplet size while improving monodispersity. However, all of these simulations and experiments studies in the aforementioned articles mainly focused on the generation process of water droplets or oil droplets. Compared with these conventional fluids, liquid LMPA has obviously higher interfacial tension and density, so it is necessary to study how different factors influence the formation of LMPA droplets. Some scholars have experimentally demonstrated the feasibility of generating LMPA droplets using the microfluidic device (Hutter et al., 2012;Tian et al., 2017;Tian et al., 2021). For example, Hutter et al. (2012) reported the generation of spherical and nonspherical eutectic gallium indium (EGaIn) liquid metal droplets in a microfluidic device around room temperature. The results show that the volume of droplets was determined by the channel size and the applied flow rates, and nonspherical LMPA droplets could be produced in oxidized silicone oil. Tian et al. (2017) proposed a liquid metal-based microfluidic device to produce liquid droplets and actively sort the liquid metal droplets using dielectrophoresis. The results showed that the microsystem could produce uniform spherical metal droplets and dominate the flow orientation of the droplets. However, these literature works only focused on the formation of liquid droplets, without paying attention to the solidification behavior of droplets in microfluidic devices. In order to account for the phase change, some numerical models based on the phase-field (PF) method (Tourret and Karma, 2015;Takaki et al., 2016;Dong et al., 2020), cellular automaton (CA) approach (Yuan and Lee, 2012;Hu et al., 2020), or enthalpy-porosity method (Brent et al., 1988;Le Bars and Worster, 2006;Chakraborty and Dutta, 2011;Seredyński and Banaszek, 2021) have been explored and applied to simulate the solidification of LMPA fluid. For example, Hu et al. (2020) employed a 2D CA-finite different method model to study the microstructure transformation of columnar dendrites during the directional solidification of LMPA. The results showed that several dendrite competitive growth patterns were reproduced between two divergent and convergent dendrite series; in addition, with the rising of cooling rate and temperature gradient, the tertiary dendrites located in the dispersed grain boundary preferred to develop into new primary dendrite arms. Takaki et al. (2016) and Tourret and Karma (2015) presented 2D PF models for simulating the growth competition of columnar dendrites and examined the effects of the solute interaction, grain bicrystallography, and temperature gradient on the microstructural variation of bicrystal dendrites during the directional solidification period. Chakraborty and Dutta (2011) studied the solidification process of LMPA using the enthalpy-porosity method, and the results showed that the solid movement has a The enthalpy-porosity method converges rapidly and can produce precise results for the morphology and position of the phase change front at different times with a lower computational cost (Brent et al., 1988). Therefore, in this study, a coupled VOF and enthalpy-porosity method will be used to predict the formation and solidification properties of LMPA droplets in a flow-focusing device. The effects of flow velocity, interfacial tension, continuous-phase temperature, and wall temperature on the LMPA droplet formation and solidification process will be examined in detail.

GOVERNING EQUATIONS AND BOUNDARY CONDITIONS
For the numerical analysis of the droplet generation and solidification process, the following assumptions were made: (1) the two-phase fluid was Newtonian and incompressible; (2) the flow in the microchannel was unsteady and laminar, ignoring viscous dissipation; (3) the thermal properties of the liquid and solid phases were assumed to be constants; (4) the 3D convection was not taken into account, as a 2D model was used in this study.
The numerical method in this study is based on the VOF and enthalpy-porosity method of the commercial code FLUENT. The governing equations for the mass conservation and momentum conservation are given as follows: where ρ is the density; V is fluid flow velocity; and t, P, and μ are, respectively, the time, pressure, and dynamic viscosity. F is a source term and equal to the interfacial tension force in this study since acceleration can be ignored in this micrometer scale and F F σ . The interface tracking between the phases was achieved by solving the following volume fraction continuity equation: The density and dynamic viscosity in equation Eq. 2 can be calculated as follows: where α is the proportion of the fluid phase in the unit grid and subscript C and D represent the continuous phase and discrete phase, respectively. The continuum surface force (CSF) method was applied to achieve the interfacial surface force term in Eq. 2, where: where σ is the surface tension coefficient and κ is the local curvature which is given by: wheren is the surface unit norm defined as follows: Defining the contact angle at the channel wall θ to consider wall adhesion and calculate the surface normal at the reference cell close to the channel wall was performed using the following formula:n n w cosθ +t w sinθ , wheret w andn w are the unit vectors tangential and normal to the wall, respectively. Energy equation can be written as Eq. 9: where T is temperature, k is the thermal conductivity, and H is the specific enthalpy. The specific enthalpy can be received by the sum of the enthalpy change caused by the phase change γL and the sensible enthalpy, H H ref + where H ref is the reference enthalpy, T ref , C P is the specific heat, L is the specific melting enthalpy, and γ is the liquid fraction during the phase transition which happens in theT S < T < T L region, defined by Eq. 10: , where T S and T L are the liquidus temperature and solidus temperature of the dispersed phase, respectively, which are 285.85 and 288.85, respectively, as shown in Table 1.
The source term S in the momentum equation is defined as follows: where A is a constant, which is used to describe how sharply the velocity is suppressed during phase material solidification, and it can be normally changed from 10 4 to 10 7 (10 6 is used in the present study). ε is a small value (0.001) to prevent divergence. The boundary conditions are as follows: Inlet: Outlet: Wall: where V C and V D are, respectively, the velocity of the continuous phase and dispersed phase; T C and T D are the temperature of the continuous phase and dispersed phase, respectively; P atm is the atmospheric pressure; and T W is the cooling wall temperature. It is complex and time-consuming to conduct a 3D simulation for the process of droplet generation. Although the 2D model is vulnerable to some errors, its computational cost is lower than that of the 3D model; the 2D plot of the dispersed phase versus continuous versus contours in the literature (Serra et al., 2007) indicated that the depth curvature can be negligible and it is sufficient for the 2D assumption to eliminate the general features of the system. The 2D numerical method used in this study is validated in Section 4 by comparing it with the experimental results of the LMPA droplet formation in a 3D microchannel junction proposed by Hutter et al. (2012).
The equations are solved using Ansys Fluent 2021R1 software. A two-dimensional model is established. The PRESTO! scheme is employed to compute the pressure equation. The pressure-velocity coupling term is discretized using the PISO scheme and the second-order upwind scheme is adopted for the momentum equation and the interpolation of the interface geometry is solved using the Geo Reconstruct method.

SIMULATION SETUPS
A schematic view for a cross-flow-focusing device case model used for simulation in the present study is shown in Figure 1A. Material properties and the initial values used in the simulation are given in Table 1. Liquid LMPA which has been heated to 290 K was used as the dispersed phase and silicone oil was used as the continuous phase. The LMPA liquid phase entered from inlet 1, and the silicone oil phase entered from inlets 2 and 3. In this study, the rectangular mesh was used as shown in Figure 1B.
The dependence on the cell size and time steps was examined by using four mesh sizes (4,802, 12,046, 23,498, and 42,418) and three time steps (0.0001, 0.00001, and 0.000001 s), and the droplet-generation frequency was almost the same with 23498 and 42418 mesh elements or with 0.00001 and 0.000001 s time steps as shown in Tables 2, 3. Therefore, a step size of 0.00001 s and the model with 23,498 cells were used for the balance of precision and computational cost.

MODEL VALIDATION
As there is no suitable reported experimental result for simultaneously realizing the formation and solidification process of the LMPA droplet in the microfluidic device, in order to verify our simulation method, we only simulated the formation process of droplets and compared it with the relevant experimental results.
The coupled VOF and enthalpy-porosity method were used to numerically simulate the formation of LMPA liquid droplets in a flow-focusing microchannel and verified by the experimental results of the microfluidic device (Hutter et al., 2012). Figure 2A shows the design of the microfluidic device with 80-µm nozzle geometry, the depth of the channel was 75 μm, and Figure 2B shows the geometric structure for simulation. The droplets were made of EGaIn (75.5% Ga and 24.5% In by weight) liquid metal solution and the continuous phase was silicone oil. As the working fluids and geometry of this experimental study were very similar to our designed microchannel, these experimental results were selected for verification in the present study. In the simulation, the initial temperature and the inlet temperature of both phases were kept at 300 K, and the wall was insulated. Therefore, the droplets had no heat transfer in the microfluidic channel. A value of 0.09 N/m for the interfacial tension between the LMPA phase and silicone oil phase tested by Hutter et al. (2012) was used in this simulation, and other material properties used are given in Table 1. In order to achieve the numerical results, the simulations were carried out by considering several flow ratios (the ratio of the volumetric flow rate of continuous phase Q C to the volumetric flow rate of dispersed phase Q d ). The frequency was obtained by the quantity of the produced droplets between two similar frames and the time between these two frames. The optimum contact angle should be ascertained before the simulation. Therefore, in this study, different contact angle values of θ, including 0°, 45°, 90°, 135°, 150°, 160°, and 180°, were carefully simulated at a given flow rate before performing the verification. For an assigned QC QD 1 and the total flow rate Q Q C + Q D 2000 μL h −1 , when the contact angle is 160 • , the droplet-generation frequency is 312.5, which is almost the same as the value of 315.97 obtained in the experiment (Hutter et al., 2012) and the simulation results were the most consistent with the experimental results in terms of shape and droplet length, as shown in Figure 3. Therefore, the contact angle was assumed to be 160°in all simulations.
After ascertaining the optimal contact angle value, the simulations were carried out. In Figure 4A, the dropletgeneration frequency obtained by the numerical simulation was in good agreement with the experimental results. Figure 4B showed that the droplet lengths obtained by simulations were slightly smaller than those achieved by the experiment. It is possible that some variations of experimental physical properties (such as small fluctuations of temperature, flow rate, or incomplete straight channel wall) might cause the deviation between observed droplets and simulated sizes. The roughness of  the downstream channel might also influence the length of larger droplets. In this part, the effect of the mesh size was examined using the number of mesh elements of 4,802, 12,046, 22,823, and 40,441, and the droplet-generation frequency was almost the same with 22,823 and 40,441 mesh elements as shown in Table 4. Therefore, an optimal grid resolution containing 22,823 elements was used for the balance of precision and computational cost. The effect of time steps was also checked, and it is shown that the droplet-generation frequencies were 0.00001 and 0.000001. Thus, simulations were carried out using a time step of 0.00001 s.

RESULTS AND DISCUSSION
To study the effect of the continuous-phase velocity, interfacial tension, continuous-phase temperature, and cooling wall temperature on the droplet formation and solidification process, we use various combinations of parameters for simulation, and the different cases are summarized in Table 5. Figures 5A,B, respectively, show the simulation results of the LMPA droplet formation and solidification process in the flow-focusing microfluidic device for case 2, where the dispersed-phase temperature at the inlet is set as 290 K in all considered cases. The mixed phase in the phase-transition state is described by a liquid fraction. As shown in Figure 5B, the continuous phase is always in a liquid state; for the dispersed phase, liquid fraction F 1 represents liquid LMPA, 0 < F < 1 represents solidifying LMPA, and F 0 represents the solidified LMPA droplet. In the early stage, there is no solidified LMPA, as the dispersed-phase temperature inside the droplet is still higher than solidifying temperature and the LMPA droplet is still liquid. As the formed droplet goes along the microchannel, the LMPA solidification ratio increases gradually. As shown in Figure 5C, the temperature of the liquid droplets in the microchannel is obviously higher than that of the continuous phase. In addition, the temperature of the liquid droplets gradually decreases along the microchannel, and there exists a critical position in the microchannel where the liquid droplets are completely solidified. Figure 6 shows the temperature curves of the three different monitor points (monitor points 1, 2, and 3 as shown in Figure 1A) in the microchannel with time. It can be seen that after a certain period of time, the temperatures of the three monitor points change periodically with basically fixed frequency and amplitude. This is because there exists a relatively large temperature difference between the continuous phase and the dispersed phase. When there is no liquid droplet observed at the monitor point, the temperatures at the monitor points are greatly influenced by the cooling wall temperature and the continuous phase temperature, which will result in a lower temperature. However, when the LMPA droplets pass by the monitor points, the liquid droplets have not received enough time to reach the thermal equilibrium with the surrounding continuous phase. Therefore, the temperatures are still high in the monitor positions. However, as can be seen from Figure 6, the maximum and minimum temperature values of monitoring points 1, 2, and 3 decrease in turn, which is because the droplets exchange heat with the surrounding continuous phase during passing through these monitor points. As only when the droplets reach the monitor points, the temperature at the monitor points will rise, while the temperature at the monitor points will drop after the droplets leave. Therefore, the required time for producing a single droplet can be calculated by Δt t 1 − t 2 , as shown in Figure 6, where t 1 and t 2 are the times when two consecutive droplets pass through monitor point 2, respectively. The generation frequency of LMPA droplet is the reciprocal of Δt.

Effects of Continuous-Phase Velocity
In this section, we changed the velocity of the continuous phase from 0.02 m/s to 0.08 m/s and kept all the other parameters unchanged as Table 5 (cases 1, 2, 3, and 4). The visual comparison for the generation and solidification process of LMPA droplets at different continuous-phase velocities is given in Figure 7A. It is can be seen that the continuous-phase velocity greatly affects the droplet-generation process, and the droplet length gradually decreases with increasing continuous-phase velocity, which is due to the shear force of the continuous phase as a pressure acting on the dispersed phase to detach the droplet. Table 6 shows that the frequency of droplet generation increases as the continuous-phase velocity increases. This is because the shear force leads to faster splitting into the continuous phase, thus increasing the droplet formation frequency. These results are consistent with the previous experimental data (Kovalchuk et al., 2018). In addition, Figure 7A also shows that at a lower flow rate, the droplets are completely solidified closer to the inlet, which means that increasing the continuous-phase velocity can increase the distance between the first solidified droplet and the inlet of the microchannel under given conditions. Figure 7B indicates that the temperature distribution of the continuous phase in the channel is almost unaffected by the continuous-phase velocity. However, it is clearly given in Figure 8 that the time required for a single droplet from initial formation to complete solidification first increases and then decreases with the increase in the continuous-phase velocity, which is possibly due to the comprehensive effect of lower temperature when the larger droplets are generated and the rapid heat transfer rate between the smaller droplets and the surrounding environment. Within the given continuous-phase velocity range, LMPA droplets with relatively large size are produced when the velocity is less  than or equal to 0.04 m/s (cases 1 and 2), and the temperature of the larger droplets (case 1) is lower than that of the smaller droplets (case 2) just produced. In this case, the influence of the lower temperature on the solidification time exceeds that of the rapid heat transfer between the smaller droplets and the surrounding environment, and thus the time required for a single droplet solidification increases with the increase in continuous-phase temperature. On the contrary, when the continuous-phase velocity is greater than or equal to 0.04 m/s (cases 2, 3, and 4), although the larger droplets just formed have relatively lower temperature, the rapid heat transfer between the smaller droplets and the surrounding environment has more influence on the solidification time than the lower temperature, and thus the time required for a single droplet solidification decreases with the increases in continuous-phase temperature.

Effects of Interfacial Tension
A previous study shows that the interfacial tension has a significant effect on droplet stability and break up time (Zhu and Wang, 2016). In this section, we changed the interfacial tension from 0.05 to 0.16 N/m and all the other factors were kept unchanged as shown in Table 5 (cases 2, 5, 6, and 7). The visual comparison for the LMPA droplet generation and solidification at different interfacial tensions is shown in Figure 9A, which distinctly indicates that the droplet size increases with the increasing interfacial tension. For the frequency, Table 7 shows that it constantly decreases as the interfacial tension grows. This is because the viscous and shear force are suppressed at a higher surface tension value, which leads to the inward contraction of the liquid droplets due to the pressure affected by the interfacial tension. The higher interfacial tension between the two liquids contributes to lower frequency and higher droplet size. On the other hand, when the interfacial tension between liquids is lower, the droplet with higher frequency and size is obtained. Figure 9B shows that the effect of interfacial tension on the temperature distribution of the continuous phase in the microchannel is not obvious. However, as shown in Figure 10, the required time for droplet solidification constantly decreases with the increasing interfacial tension. This is mainly because, within a given interface tension range, the influence of the temperature at the initial formation of droplets on the solidification time exceeds the rapid heat transfer between smaller droplets and the surrounding environment. Therefore, the interface tension and the time required for droplet solidification are linearly related under given conditions.

Effects of Continuous-Phase Temperature
In this section, we changed the continuous-phase temperature from 250 to 265 K, and all the other parameters were kept unchanged as shown in Table 5  ( cases 2, 8, 9, and 10). The comparison of droplet formation and solidification progress at different continuous-phase temperatures is shown in Figure 11A, which clearly indicates that the LMPA droplet size almost keeps unchanged as the continuous-phase temperature increases. In addition, Table 8 shows that the dropletgeneration frequency also almost stays unchanged with the increase in continuous-phase temperature. This is possibly because the influence of the continuous-phase temperature on the shear force or the velocity is not obvious under a given temperature difference and thus barely affects the frequency and droplet size. However, the effect of continuous-phase temperature on the droplet solidification process is relatively obvious, and it is observed from Figure 11A that the droplet that is almost completely solidified in the microchannel at different continuous-phase temperatures is located in different positions in the microchannel. This is because increasing continuous- phase temperature will reduce the temperature gradient along the main microchannel as it is shown in Figure 11B, and thus, it can be observed that the droplets will solidify closer to the inlet at a lower continuous-phase temperature. Figure 12 shows that the time required for a single droplet solidification increases with the increase in continuous-phase temperature, which is mainly because increasing continuous temperature will reduce the heat transfer between droplets and the environment, and further slow the solidification process.

Effects of Cooling Wall Temperature
In this section, the cooling wall temperature was changed from 250 to 265 K, and all the other parameters were kept unchanged as Table 5 (cases 2, 11,12,13). The comparison of LMPA droplet formation and solidification progress at different cooling wall temperatures is shown in Figure 13A, which also clearly indicates that the droplet size almost stays unchanged as the cooling wall temperature increases. Table 9 shows that the droplet-generation frequency almost stays unchanged with the increase in cooling wall temperature. This is possibly because the influence of the cooling wall temperature on the shear force or the phase velocity is not obvious under a given temperature difference and thus  barely affects the frequency and droplet size. However, similar to the effect of continuous-phase temperature on the droplet solidification process, it is observed from Figure 13A that the droplet that is almost completely solidified in the microchannel at different cooling wall temperatures is located in different positions in the microchannel. This is because increasing cooling wall temperature will reduce the temperature gradient along the main microchannel as shown in Figure 13B, and thus, it can be observed that the droplets will solidify closer to the inlet at a lower continuous-phase temperature. Figure 14 shows that the time required for a single droplet solidification increases with the increase in cooling wall temperature, which is mainly because increasing cooling wall temperature will reduce the heat transfer between droplets and the environment, and further slow down the droplet solidification process.

CONCLUSION
In this study, a simplified two-dimensional cross flow-focusing geometry model is given. Then, a combination of the VOF method and enthalpy-porosity method based on the FLUENT code was applied to examine the effects of some key factors on the LMPA droplet formation and solidification process in a microfluidic device. The main results are as follows: • Increasing the continuous-phase velocity reduces the LMPA droplet size and simultaneously increases the droplet-generation frequency. The variation of continuous-phase velocity has little effect on the temperature distribution in the microchannel. The required time from droplet generation to complete solidification first increases and then decreases with increasing continuous-phase velocity.
• Increasing the interfacial tension can increase the LMPA droplet size and simultaneously reduce the dropletgeneration frequency. The interface tension and the time required for droplet solidification are linearly related. The higher the interfacial tension, the shorter the required time from droplet generation to complete solidification. • The continuous-phase temperature has no obvious influence on the size and frequency of the generated droplets. However, it has a significant influence on the temperature distribution and the process of LMPA droplet solidification in the microchannel. With the increase in continuous-phase temperature, the required time for droplet solidification gradually decreases. • Similar to the influence of continuous-phase temperature on the droplet-generation and solidification process, the wall temperature has no obvious effect on the size and frequency of generated LMPA droplets. However, it has a great effect on the temperature distribution and the process of droplet solidification in the microchannel. With the increase in wall temperature, the required time for droplet solidification gradually decreases.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YR proposed the research topic and supervised the research project. YG and TH conducted the numerical modeling and simulation. YG completed the data analysis and drafted the manuscript. JW and CW reviewed the simulation results. All authors reviewed the manuscript.