Numerical simulation of oscillating flow field including a droplet

The oscillating flow field including a droplet is simulated numerically by solving the Navier-Stokes equations using the arbitrary Lagrangian-Eulerian and level set coupled method. The oscillating pressure and velocity fields are shown to be simulated well. Parallel computations with a fine grid system are performed, and it is found that the surface-layer flow on the droplet is oscillating with a different timing and reversed before the reversal of the bulk flow. The acoustic streaming with large vortices is shown to appear in the average flow field.


INTRODUCTION
Precise measurement of properties of high-temperature molten material is of practical importance in many engineering fields [1][2][3].Material properties of molten mixtures of nuclear fuel are, for instance, indispensable for evaluation of severe accident scenario of nuclear reactors [1].Measurement techniques using containers are not suitable, since molten mixtures are highly reactive and corrosive, and affected by the container walls.A levitated liquid droplet is used to measure material properties of molten material at high temperature, since the levitated droplet is not in contact with a container, and the effect of container wall is eliminated for precise measurement [4].The levitation of liquid droplet, which is also used for container-less processing of material, is controlled by using electromagnetic [5] electrostatic [2,3,6], ultrasonic [7], or aerodynamic force [1] in the vertical direction.The levitated droplet is affected by an oscillating pressure field when the ultrasonic or acoustic levitator is used.Additionally, rotation is sometimes imposed on the droplet to stabilize its motion [3,8].The rotational motion is controlled by magnetic [3] or acoustic force [8] perpendicular to the vertical axis.It is thus important to know the oscillating flow field including a droplet.Streaming flows associated with ultrasonic levitator have been studied experimentally [9,10] and theoretically [11].Although steady-state flows in and around the droplet were studied, transient flow fields have not yet been discussed well.
The oscillating flow field including a droplet is simulated numerically in this study.The droplet is placed in the center of the simulation region, and the oscillating flow field is modeled with oscillating boundary conditions.The flow field is calculated using the arbitrary Lagrangian-Eulerian (ALE) method, where the computational grid points are moved with the oscillating boundary velocity.Locations of droplet surface are obtained using the level set method.In the level set method, the level set function, which is the distance function from *Corresponding author.E-mail: twata@u-fukui.ac.jp the droplet surface, is calculated by solving the transport equation using the local flow velocity.Incompressible Navier-Stokes equations are solved to obtain the flow field, and the effect of compressibility is taken into account as the density perturbation.The simulation program is parallelized by the domain decomposition technique using the Message Passing Interface (MPI) library.The pressure field in and around the droplet is shown first to see the acoustic standing wave including the droplet.Transient velocity fields in and around the droplet are shown next, and characteristics of the oscillating flow field especially in the thin layer on the droplet surface are discussed.

GOVERNING EQUATIONS
Governing equations for the flow field including a droplet are the equation of continuity and the incompressible Navier-Stokes equation [12]: (1) and (2) where ρ, u, p and µ, respectively, are the density, the velocity, the pressure and the viscosity, D is the viscous stress tensor, and F s is a body force due to the surface tension.External force fields such as the gravity are not simulated in this study.The surface tension force is given by (3) where σ, κ, δ, and φ are the surface tension, the curvature of the interface, the Dirac delta function and the level set function, respectively.The level set function is a distance function defined as the normal distance from the droplet surface: φ = 0 at the surface, φ < 0 in the droplet region, and φ > 0 for the outside region.The curvature of the droplet surface is expressed in terms of φ: The density and viscosity are given, respectively, by (5) and (6) where the subscripts g and l denote gas and liquid phases, respectively, and H is the smeared Heaviside function defined by where ε is a small positive constant for which ∇φ = 1 for φ≤ ε.The time evolution of φ is given by (8) In this study, the ALE method [13] is applied, and the computational grid is assumed to be oscillating in the vertical direction with the same velocity as the boundary velocity.The substantial derivative terms in Eqns ( 2) and ( 8) are thus defined by (9) where U is the velocity of the computational grid or the oscillating boundary velocity.
In order to maintain the level set function as a distance function, reinitialization of the level set function by solving the following equation is proposed [8]: (10) where τ is an artificial time and sign (φ 0 ) indicates the sign of the level set function at the beginning of the reinitialization procedure.The level set function becomes a distance function in the steady-state solution of Eqn (10).The smoothed sign function proposed for numerical treatment of reinitialization [14] is used for Eqn (10); (11) where h is the spatial increment in the finite difference method for solving the governing equations.A smoothed version of the sign function was also used in [12].
The following equation is also solved to preserve the total mass in time [14]; (12) where A 0 denotes the total mass for the initial condition and A denotes the total mass corresponding to the level set function, P is a positive constant for stabilization, and 1.0 was used in [14].The total mass is conserved in the steady-state solution of Eqn (12).The effects of Eqns ( 10) and ( 12) on the nonlinear droplet dynamics have been discussed in [15].
The effect of compressibility is taken into account simply as the density perturbation, since the amplitude of pressure variation and thus the density variation are both small in this study.The density perturbation is defined by (13) where, c is the sound speed.

NUMERICAL SCHEME
The finite difference method is used to solve the governing equations.The staggered mesh system is used for spatial discretization of velocities.The convection terms are discretized using the second order upwind scheme and other terms by the second order central difference scheme.Time integration is performed by the second order Adams-Bashforth method.The SMAC method is used to obtain pressure and velocities [16].The pressure Poisson equation is solved using the Bi-CGSTAB method.The domain decomposition technique is applied in the vertical direction of the simulation region and the MPI library is used for parallel computations.The block Jacobi preconditioner is used for the parallel Bi-CGSTAB method [17,18].

SIMULATION MODELS AND CONDITIONS
The simulation region is a two-dimensional rectangular region shown in Figure 1.The vertical size of the simulation region is 6.4 mm for checking the pressure field and the parallel efficiency, and 12.288 mm for simulating the oscillating flow field.The horizontal size of the simulation region is 3.2 mm for both cases, and the radius of the droplet is 2 mm.A half of the droplet is simulated as shown in Figure 1.Water and air densities and viscosities are assumed for the inside and outside of the droplet, respectively.The slip and symmetric boundary conditions are applied at the right and left sides of the simulation region, respectively, and the oscillating velocity conditions given by the following equation are used at the top and bottom boundaries: (14) where the oscillation amplitude A is 0.076 m and the angular frequency ω is 1.26 × 10 5 rad/s.These values are almost corresponding to the sound pressure of 0.25 kPa and the sound frequency of 20 kHz in air.The sound speed is 340 m/s for air and 1400 m/s for water.The oscillation velocity given by Eqn ( 14) is also used as the grid velocity in Eqn (9) used in the ALE method.This indicates that the simulation region and the calculation grid are oscillated according to (15) where X is the vertical position.The grid size is 3.2 µm both in horizontal and vertical directions, and the number of grid points is 1000 × 2000 for checking the pressure field and the parallel efficiency, and 1000 × 3840 for simulating the oscillating flow field.The time step size is set equal to 0.08 µs, and 625 time steps are corresponding to one oscillation period.

VARIATION OF PRESSURE DISTRIBUTION
Variations of pressure distributions in one oscillation period, T, after about 20 oscillation periods, are shown in Figures 2 and 3.The vertical pressure distributions along the right side of the simulation region are shown in Figure 2, where the pressure variation in the air without the droplet is indicated.The pressure gradient is almost the maximum at 0.2 T, 0.3 T, 0.7 T and 0.8 T, with the pressure node at the center of the simulation region.The maximum pressure gradient is actually at 0.25 T and 0.75 T, though not shown in Figure 2, and almost zero at 0.5 T and 1.0 T. The pressure distributions along the vertical center axis of the droplet are shown in Figure 3.The timings of the maximum and zero pressure gradients are the same as those shown in Figure 2. The pressure distribution is, however, affected by the droplet, and the pressure gradient is larger in the droplet than in the air [19].It is found that the oscillating pressure field including the droplet is simulated well by the present numerical method using the ALE and level set coupled method with a pseudo compressibility.

PARALLEL EFFICIENCY
A fine mesh system is used in this study for simulating small scale flow phenomena around the droplet surface.The simulation program is parallelized by the domain decomposition   method using the MPI library, and parallel computations using 256 processors are performed.The speedup of parallel computation is shown in Figure 4 for 1000 × 2000 grid calculation.
The total calculation time including I/O process is indicated as "Elapsed Time," and the time for matrix calculation in the Bi-CGSTAB method is shown as "Matrix Calculation."The number of time steps is five in one calculation, and the calculation time is obtained as an average of five calculations.The average calculation time using single processor is 3345.6 s for "Elapsed Time" and 3316.8 s for "Matrix Calculation."The time for matrix calculation is thus more than 99% of the total calculation time in this study.The number of iterations for the Bi-CGSTAB method is different in each time step and varied from 668 to 1883 for the single processor case.The speedup or parallel performance is 64.5 for "Elapsed Time," and 76.5 for "Matrix Calculation" using 256 processors.The number of iterations is increased as the number of processors increases.The number of iterations for the 256 processor case is varied from 1459 to 9387.The time for data communication, thus, becomes long as well as the calculation time, as the number of processors increases.Actual calculation time for 625 time steps using 256 processors is about 4300 s.The domain decomposition method is applied in the vertical direction in this study.The number of vertical grid points is 2000 in the test case shown in Figure 4, and the number of grid points in each processor is not always the same.For instance, 250 grid points are calculated in one processor if 8 processors are used, but if 32 processors are used, 62 and 63 grid points are calculated in each of 16 and other 16 processors, respectively.It is indicated in Figure 4 that the parallel efficiency is not so good for the cases with 32, 128, and 256 processors.For 128 processors, 15 grid points are calculated in each of 48 processors and 16 grid points are calculated in each of other 80 processors.The calculation time is thus relatively shorter for 48 processors, and longer for 80 processors.In other words, 48 processors have to wait for other 80 processors.The parallel efficiency is thus not good for the case in which the number of grid points cannot be divided by the number of processors.This is also the case with 256 processors.For studying the oscillating flow field in the following, large-scale simulations are performed with 3840 vertical grid points using 256 processors, and 15 grid points are calculated in each processor.

OSCILLATING FLOW FIELD
Simulated flow fields in one oscillation period in the simulation region of 3.2 mm × 12.288 mm are shown in Figures 5 and 6.The pressure field is in Figure 5, where the high and low pressure regions are indicated by red and blue, respectively.The pressure field is almost flat at 0.000 T, 0.499 T and 1.000 T. The pressure in the lower part is, however, higher at 0.248 T and lower at 0.752 T, while that in the upper part is lower at 0.248 T and higher at 0.752 T. It is seen that the pressure field is oscillated vertically with the pressure node at the mid elevation in Figure 5.The pressure distribution along the droplet center is gradually varied in the air to the pressure distribution along the right side.It is noted that the pressure field with the maximum gradient at 0.248 T or 0.752 T is qualitatively the same as the sound pressure field reported in [19].
Corresponding velocity field is shown in Figure 6, where the distribution of vertical velocity component is shown to see the vertical oscillating flow field clearly.Upward and downward velocity regions are indicated by red and blue, respectively.The bulk flow direction is downward at 0.000 T and 1.000 T, while upward at 0.499 T. The maximum velocity thus corresponds to the flat pressure distribution shown in Figure 5.This is because the pressure gradient is related to the time derivative of the flow velocity as given by the Navier-Stokes equation, Eqn (2).The 0.000 T 0.248 T 0.499 T 0.752 T 1.000 T flow velocity in the droplet is calculated to be very small in comparison with the outer flow velocity, and seems to be zero in Figure 6.It is found that the velocity field is much affected by the droplet.The flow velocity becomes larger at the side of the droplet regardless of the flow direction, since the flow area in the air is small due to the droplet.The flow velocity is, on the contrary, small at the top and bottom of the droplet due to the flow stagnation.It is found in Figure 6 that the velocity field on the droplet surface is different from the bulk velocity field especially at 0.248 T and 0.752 T. The bulk flow velocity is almost zero at these timings, and the pressure gradient is the maximum as shown in Figure 5.The thin red region is, however, seen on the side surface of the droplet at 0.248 T. This indicates that the upward flow appears in the surface layer at 0.248 T, before the bulk flow becomes upward at 0.499 T. The blue downward flow region is appears next at 0.752 T, before the bulk flow becomes downward at 1.000 T. The velocity vectors at the side surface of the droplet are shown in Figure 7, where the small region with several grid points is enlarged to see the surface layer clearly.It is noted that the magnitude of the velocity vector in the droplet, which is shown in the left part of each figure, is 100 times larger than the outside vector in the right part.It is seen at 0.000 T that the inside and outside flows are both downward and the flow velocity varies gradually across the inside and outside boundary layers at the droplet surface.The outside flow is slightly upward at 0.248 T, while the inside flow is still downward, and the flow velocity is larger near the surface as indicated in Figure 6.The outside flow velocity becomes large at 0.499 T, and the inside flow is also upward.The outside flow is slightly downward at 0.752 T, while the inside flow is still upward.The downward velocity is larger near the droplet surface as shown in Figure 6.The outside velocity is large and the inside velocity is downward at the end of one oscillation period at 1.000 T. The surface-layer flow seems to change its direction before the bulk flow.The oscillating flow velocity in the surface layer is smaller than the bulk flow velocity, and sensitive to the outside pressure field.The surface-layer flow is, thus, reversed first according to the pressure field.

FORMATION OF ACOUSTIC STREAMING
Average flow field is shown in Figure 8, where the time average is taken over one oscillation period.The average flow velocity is almost zero in one oscillation period, since the magnitude of upward and downward flow velocities are exactly the same at the top and bottom boundaries.It is, however, shown in the average flow field that the steady vortex flow appears.The steady flow in the oscillating pressure field is known as the acoustic streaming, and observed in the experiments [9,10] and in the theoretical analysis [11].The recirculation flow region on the droplet was also predicted by the theoretical treatment, Numerical simulation of oscillating flow field including a droplet 0.000 T 0.248 T 0.499 T 0.752 T 1.000 T  though it is not clear in the experiments and in this study.The acoustic streaming is due to the interaction between inertia and viscosity [20], and originated in the surface layer.The average flow field is obtained after about 20 oscillation periods in this study, and the flow velocity is still very small.The magnitude of average flow velocity is about 100 times smaller than the oscillating flow velocity.Longer time calculations might be necessary to obtain the recirculation flow region.

CONCLUSION
In this study, the oscillating flow field including a droplet has been simulated numerically by solving the Navier-Stokes equations using the arbitrary Lagrangian-Eulerian (ALE) and level set coupled method.Detailed transient flow fields in and around the droplet were calculated by performing parallel computations with 256 processors, and the flow phenomena in the thin layer around the droplet surface were discussed.The pseudo compressibility was assumed for the incompressible governing equations, and the oscillating pressure field was simulated well as well as the oscillating velocity field.It was found that the oscillating flow also appeared in the surface layer, where the oscillation timing was different from that in the bulk flow region.This surface-layer flow was formed by the pressure field, and reversed before the reversal of the bulk flow in each oscillation.The acoustic streaming with large vortices was shown in the average flow field, though the recirculation flow region predicted by the theory was not clear.The relation between the large-scale acoustic streaming and the small-scale surface-layer flow would be made clear in the near future.

Figure 1
Figure 1 Schematic diagram of droplet and simulation region.

Figure 2
Figure 2 Vertical pressure distribution along the right side of the simulation region.

Figure 3
Figure 3 Vertical pressure distribution along the center axis of the droplet.

Figure 4
Figure 4 Speedup of parallel computation.

Figure 7
Figure 7 Velocity vectors at the droplet surface.