Fluid-fluid Interaction Study of Different Density and Viscosity Using Smoothed Particle Hydrodynamics Method

In Computational Fluid Dynamics, fluid-fluid interaction is used to model the behaviour of two different fluids as a mixture. The numerical simulations of fluid-fluid interaction were carried out by using Smoothed Particle Hydrodynamics (SPH) method based on Navier-Stokes equation in FORTRAN programming language. Correct assumptions and other suitable conditions are needed to reach best fluid-fluid interaction simulation. The term of viscosity force equation was written to model the multiple-fluid interactions. This research uses two different fluids in density and viscosity within a three-dimensional cube container. In order to get a pure interaction of fluids, other effects and influences like chemical reaction and heat transfer are neglected. The fluids are miscible with each other. The fluid particle collisions to container are regarded as perfectly elastic collisions. The simulations of two different fluids in this research showed different particle movements in each numerical simulation, but similar to its actual behaviour.


Introduction
In Computational Fluid Dynamics (CFD), the simulations can occur between fluid-fluid or solid-fluid interactions. These two ways of interactions will give different results. Fluid-fluid interaction is commonly used to model the behaviour of two different fluids as a mixture. The correct assumption and specific conditions is still researched until today. The one common method used in CFD problems is Smoothed Particle Dynamics (SPH) Method which uses Lagrangian approach.
The SPH method first invented to simulate non-axisymmetric phenomena in astrophysics (Lucy 1977, Gingold & Monaghan 1977 [1]. The SPH method is a simple method and can provide good accuracy. The Lagrangian approach makes the movement of particles which can be modeled well without bounds of grid and particles can move freely. By using this method, modeling of fluid behavior becomes more realistic and applicable for studying fluids with multiple phases.
SPH modeling with varying fluid phases has a few differences where the formulations are affected by the conditions of each phase. Basically, the equations must still satisfy the equation of continuity and momentum equation. In this multi-phase fluid model, the properties of the fluid will not be declared globally, but will be stored in each particle. The characteristic of fluids also affects whether they are miscible or immiscible.
This research is emphasized to analyze the interaction between fluid-fluid. This research uses some simple numerical simulations in a three-dimensional cube container for simulation, with two different 1234567890''"" fluids in density and viscosity. In order to get a pure interaction of fluids, other effects and influences like chemical reaction and heat transfer are neglected. And also, the fluids are miscible with each other.
This study is a further research of SPH to develop the previous numerical program from single-phase fluid to multiple-phase fluids. The previous researches were taken by Marthanty [2] and Lydiana [3] which still used single-phase fluid in their researches. As the development, this research modifies the term of viscosity force equation from single-fluid to multiple-fluid equation used by Müller et al. [4,5] and Kelager [6]. This research was done to make qualitative validation and to study the influence of different density and viscosity to multi-phase fluid behaviour.

Smoothed Particle Hydrodynamics
SPH method was invented to simulate nonaxisymmetric phenomena in astrophysics (Lucy 1977, Gingold & Monaghan 1977) [1]. It is a simple method and can provide good accuracy. Furthermore, this method gives satisfactory results and can be applied to various physics problems.
This method is basically an interpolation method with kernel functions. The integral of any function A(r) is defined by: where r is particle position and W is kernel function (smoothing kernel) with width h (smoothing radius). In numerical form, the function can be written: where m is mass of each particle and ρ is density of each particle.
The SPH method can estimate the properties of fluid and derivation values of a continuous field based on discrete points called smoothed particles. The property of a fluid particle is computed by interpolating the property of surrounding fluid particles in a kernel radius. The derivation value of a fluid's property is calculated by deriving the kernel function W:

SPH Equations for Fluid-Fluid Interaction
The Navier-Stokes equation is used as the governing equation in this method. The Navier-Stokes equation for incompressible and isothermal fluid is: where ρ is density, u is velocity vector, p is pressure, μ is viscosity coefficient, and f is external forces. The total forces consists of internal and external forces, so the equation becomes: Then, the acceleration of particle i can be obtained from the equation:

1234567890''""
International The internal forces consist of forces due to pressure and viscosity effects, while the external forces consist of forces due to gravity, buoyancy, and surface tension [4,5,6]. The equations to compute density, pressure (Tait equation, [1]) and forces are as follows:

Smoothing Kernel
The stability, accuracy, and speed of SPH method are strongly influenced by selection of kernel functions or smoothing kernels. The kernel functions used in this research is based from Müller et al. [4] and Kelager [6]. The equations can be seen on the appendix.

Numerical Parameters
The numerical parameters shown in Table 1 were chosen for all simulations, except for density and viscosity of fluids that are shown in Table 2. Here, each scenario of numerical simulation will consider two types of fluid with different density or viscosity to observe the interactions of different fluid properties. The buoyancy (b) is set to zero because we don't model air particles. Coefficient of restitution (cr) is used to bounce back particles that are move through the boundary, and it is set to 1 to model perfectly elastic collision.

Geometry of Container
The simulation used a three-dimensional cube container as a computational domain with dimension of 0.1 m x 0.1 m x 0.1 m. The model of boundaries (walls and plate) used perfectly elastic collision to prevent all particles leave the container.

Numerical Simulations
In each simulation, two layers of fluid are placed as shown in Figure 1. The combinations of Fluid1-Fluid2 are used for studying the influence of different density of fluid in a container. Then, the combinations of Fluid1-Fluid3 are used for studying the influence of different viscosity of fluid in a container. The two fluids are assumed at rest and the positions of two fluids are interchangeable. The number of fluid particles in a container is 3179 particles with particle spacing equal to 6.25 x 10 -3 m (or 1/16 of dimension of the container). The density was observed by computing the average from all fluid particles. The detail of simulation scenarios is described on Table 3.

Pressure Correction
Some simulations were performed to check the stability of fluid particles, and it is found that the gravity force dominates over the other forces. The gravity force made all particles to move downward and the fluid became compressed. As a result, the density increased to 3.5 times its rest density. And the other forces had only small effects in the simulation, especially pressure force which should be more active to adjust the fluid density. So, the pressure correction is needed by considering the similitude law to adjust the pressure force. The pressure force was compared to the gravity force from the value obtained in simulation and the actual value. The pressure correction coefficient was derived: with ΔpA and ρA are pressure difference and density in actual condition, and ΔpB and ρB are pressure difference and density obtained from the simulation. The pressure correction coefficient makes the pressure force increases and more comparable to the gravity force.

Results
The pressure correction increased the pressure force to achieve the stability of the fluid particles. The pressure force became more dominant than other forces and this pressure force could adjust the particle density to density at rest. Figure 2 shows the average density during simulation time from 0 to 2 s for each simulation. There is some fluctuation happened in density from the start of simulation, because the numerical computation was adjusting the density until reach the stability of computation. After some computation steps, the fluctuation decreased and the density became stable. The difference of density in Figure 2 is proportional to average density in Table 3. Simulation 1A and 1B have higher density than Simulation 2A and 2B that have similar density.

Figure 2. Average density during simulation
Difference of density affects the gravity force which is proportional to the density. The higher gravity force on a particle makes the particle moves downward. As shown in Figure 3 for Simulation 1A, the higher density fluid (lighter color) is positioned above the lower density fluid (darker color) at initial condition of simulation. As the simulation run (t = 0.4s), the higher density particles move through the lower density particles with mixing between the particles. After some time steps (t = 1.2s), most of the higher density particles are on the lower layer and the particles become stable. In the boundary of two layers, some particles mix with the other type particles because of miscible model. In Simulation 1B, the lower density fluid is positioned above the higher density fluid but there is no such position change like the previous one. The different viscosity simulations did not have significant difference for the variables except the viscosity force. Figure 4 shows the average viscosity force for each simulation. Fluid 3 in Simulation 2A and 2B has 1000 times viscosity than Fluid 1, so the effect of viscosity difference is clearly visible in the graph. Higher viscosity forces make effects on particle accelerations and movements. As shown in Figure 5, the effect of higher viscosity makes the particle movements slower and the particles have different behaviour (more viscous) than the other one.  The pressure correction factor makes the simulation could achieve the stable condition after some simulation steps. Table 4 shows the magnitude of the pressure correction factor used in each simulation and the simulation time when most of particles reach stable condition. From the results, we can simulate and model two different density or viscosity fluids but still in miscible model where chemical effect and heat transfer are neglected. All simulations in this research are executed on 1.67 GHz processor and simulate 3179 particles on 1000 time steps to reach simulation time t = 2s with average computer time consumption 12457 s.

Conclusion and Future Work
In this research, a numerical program was developed to simulate different density and viscosity of fluid by SPH method. Pressure correction coefficient was determined in this simulation to increase the pressure which makes the particle movements more stable. The results showed the particle movements and behaviour of fluid-fluid interaction similar to their actual behaviour. These results also showed that the qualitative validation had been achieved by applying the kernel functions and internal forces proposed by [4,5,6]. Further work that needs to be done is to improve the void area near the container found in Simulation 2 by using ghost particle method.
The assumptions used in this program are limited to miscible fluid; so in the future, the program will be improved for immiscible fluids with the effects of interface surface tension and more to develop fluid-solid interaction to simulate a sediment transport phenomenon.