Review Article
Investigate on the Development Direction on Simulating Surface Tension with Smoothed Particle Hydrodynamics
Department of the Computer Science Wuzhou University, No. 82, Fumin 3rd Road, Wuzhou, Guangxi, Postal code 543002, China
Fluids are almost everywhere in peoples daily life. The fluid interface motion caused by surface tension is of great significance in most nature and industrial technologies, including transfer and dispersion of underground pollutant, food processing, cosmetics production, metal casting, glazing and fuel adding, etc. The numerical simulation method is typically required to conduct the detailed analyses on surface tension, in order to assist in understanding the results caused by these nonlinear fluid behaviors (Brackbill et al., 1992).
The smoothed particle hydrodynamics (SPH) is a numerical simulation method. There are some unique advantages for it to simulate surface tension. This article points out the direction on simulating surface tension after studying the development of this simulation with SPH.
Sph is extremely suitable for simulating surface tension: The numerical simulation methods of fluid are divided into grid-based method and meshfree method. Despite the great success, the grid-based numerical methods need to generate and re-generate grids, which limit their applications in many problems. As a meshfree method, SPH is thus widely applied in the area of fluid dynamics because it truly uses no grid (Liang et al., 2012; Liang, 2011a). It is characterized by its own unique advantages when solving such as problems with free surface, deformable boundary, moving interface (for the finite difference method, FDM), etc. (Liu and Liu, 2003). No numerical diffusion is existed in fluid interface transmission, SPH then easily deals with the interfaces containing physical and chemical effects and also those irregular, mobile and even distorted interfaces (Liang, 2011b, c; Song et al., 2009). In the same way, its advantages can be adopted in three-dimensional space (Qiang et al., 2011), the SPH is thus extremely suitable for simulating surface tension.
DEVELOPMENT OF DOMESTIC AND OVERSEAS RESEARCHES ON SIMULATING SURFACE TENSION WITH SPH
As a new and gradually matured hydrodynamics numerical method, SPH has been widely accepted (Liu and Liu, 2003) and research results from applying SPH in surface tension are also continuously published in recent decades. However, a public recognized simulating surface tension method is not yet established which is proven by diversified opinions and practices from scholars. The details are as follows.
CSF model method: Surface tension simulation with SPH is divided into two kinds up to now.
The Continuum Surface Force (CSF) model proposed by Brackbill et al. (1992) is applied in the first kind of method. Based on the definition of color function, surface normal and surface curvature are obtained through such method. The surface tension per unit is translated into a force per unit volume.
Surface tension coefficient is assumed constant, the force per unit volume Fs is then written as follows within the limited range of interface thickness (Qiang et al., 2011; Morris, 2000):
(1) |
where, in the above equation, δs is a normalized function and equals to 1 when located on interface, which is generally |n|. fs is the force per unit area:
(2) |
where, σ is the surface tension coefficient; κ(x) is the curvature of interface at the position of x, ň is the unit normal to the interface. The normal n can be obtained by using:
(3) |
where, c(x) is the color function identifying each fluid in the problem; [c(x)] is the jump in c(x) across the interface; the curvature can be calculated by using:
(4) |
CSF model is firstly introduced into surface tension simulation with SPH by Morris (2000). His test case is the non-linear oscillation of a circular droplet with an initial velocity perturbation. The relative error between his results and the results using other grid methods is within 3%. His method was then followed by various other scholars (Liu et al., 2011; Muller et al., 2003; Fang et al., 2009; Zhang et al., 2009; Liu and Liu, 2005). However, Hu and Adams (2006) found, when the curvature is calculated in accordance with the method of Morris, ∇2c(x) is required to be calculated and it is error-prone in the area where the particles only have a few neighbors. They proposed a modified version of Morris and described the surface tension as a divergence of the stress tensor, which is uniquely defined by the color function c(x):
(5) |
where, α is the surface tension coefficient; I is the unit tensor; d is the coordinate dimension. While the surface tension is given as:
(6) |
Where, in accordance with their models, ∇2c(x) does not need to be calculated. The capillary wave, three-phase interaction, drop deformation in shear flow and the mesoscopic flow in a channel were simulated and sound results were obtained by using their models. In Adami et al. (2010), a color-gradient formulation with density-weighted summation and a formulation (such equation is similar to the method of Chen et al. (1999) but with the calculation improved) to calculate the divergence of interface-normal direction were put forward by them. With these equations adopted, Square-droplet deformation, circular liquid-droplet oscillation and a circular drop in a shear flow were simulated. The simulation results are in conformity with the results from analytical solution or other methods.
It is thus clear that many improvements applying CSF model in SPH method for surface tension simulation are obtained. This approach is straightforward, but the surface tension parameters are user-input parameters. And the surface curvature is also required to be calculated by explicit or implicit expression (Liu et al., 2011; Liu and Liu, 2010), which is so complex that certain great errors are easily made in calculating the front interface curvature under some circumstances. It is not an easy work to conduct those calculations for SPH (Liu and Liu, 2010; Tartakovsky and Meakin, 2005a; b).
IIF model method: Another kind of method is proposed by some scholars to avoid the shortcomings of the first kind of method. This kind of method introduces the inter-particle interaction force (IIF) into the SPH equation, which is also named as the IIF model method. The IIF model includes the short-distance repulsion and long-range attraction and this kind of attractive force between every pair of particles generates the surface tension (Liu et al., 2011). The existing IIF model method has multiple forms, respectively described.
For the first form, one certain state equation is applied. For the long-distance attraction force among atoms is taken into account in the state equation of van der Waals (vdW), the fluid surface tension behaviors can be simulated by such form. In the vdW equation of state, the pressure is given by as follow:
(7) |
where, in the above equation, ρ is the fluid density. is the Boltzmann constant, m is the mass of the particles. = b/m, = a/m2, a and b are the van der Waals constants. In the equation, (-ρ2) item generates the attraction force causing surface tension in pressure gradient calculation process. The vdW drops performing small amplitude oscillations is simulated by Nugent and Posch (2000) using this kind of method. The simulation results are in good agreement with analytical solutions or experimental results. Melean et al. (2004) studied how to remove tensile instability from the simulation of viscous liquid and their study made further improvements for Nugent and Posch (2000) method; later, they used their method to model the formation of a circular van der Waals liquid drop and the results most similar to the experimental ones are obtained (Melean and Sigalotti, 2005).
For the second kind of form, IIF is introduced when applying the van der Waals state equation at the same time. Tartakovsky and Meakin (2005a; b) introduced the following interaction forces:
(8) |
where, in the above equation,|rj-ri|≤h, ri and rj is the position vector of the particles i and j, respectively. h is the smoothing length. sij is the strength of the force acting between particles i and j, which is determined by users in accordance with their experience. It can be seen that this interaction force becomes repulsive force at short-distance and attractive force at long-distance. Accordingly, five different numerical experiments were conducted: (1) Small-amplitude shape fluctuations of a liquid drop, (2) Measurement of the pressure in the SPH liquid as a function of drop radius, (3) Capillary rise of fluid in a small aperture, (4) Gravity-driven flow of droplets between two parallel walls and (5) A flow through an inverted Y-shaped fracture junctions and sound results were achieved.
For the third kind of form, IIF is directly added to simulate surface tension. The method was proposed by Zhou et al. (2008). The van der Waals equation was not adopted, but IIF is directly introduced into simulating surface tension. Such kind of interaction force is only existed in different kinds of particles.
The interaction force of different kinds of liquid-liquid particles is written as:
(9) |
where, in the above equation, rij is the distance between the particles i and j. S is the positive coefficient about surface tension. It is clear that it is just a kind of simple repulsive force.
When the solid can be wetted, the interaction force of different kinds of liquid-solid particles is obtained by:
(10) |
where, r0 is the initial particle spacing. n1 and n2 are, respectively the indexes of repulsive force and attraction force (generally, n1 = 12, n2 = 4). A is the positive coefficient about the wettability of interface. The equation stands for the repulsive force with short distance; and for the attraction force with long distance. In this formulation, the interaction force is repulsive at short range and attractive at long-distance.
While for the solid cannot be wetted, the interaction force of different kinds of liquid-solid particles is obtained by:
(11) |
where, R is the positive coefficient about repulsive strength and it is just a kind of simple repulsive force.
Deformation of a liquid drop in suspension was modeled and the simulation results are in agreement with the ones using the function of volume of flow (VOF). Becker and Teschner (2007) also directly introduced the inter-particle interaction force (IIF), but the forms are slightly different. The force introduced by them is just the attraction force and is only existed between the same kinds of particles. They scaled the attractive forces using the smoothing kernel as a weighting function. The square drop deformation was also modeled by Becker and Teschner (2007) and better results than those obtained by Hu and Adams (2006) were finally achieved.
From the above forms of IIF model, IIF model implicitly simulates the surface tension behaviors through the interaction force between particles and then the fluid surface does not need to be located, the surface curvature does not needed to be calculated and also the surface tension parameters do not needed to be used. Such interaction force is either indirectly generated from state equation or obtained through direct introduction.
COMPARISON ON BOTH ADVANTAGES AND DISADVANTAGES BETWEEN CSF MODEL METHOD AND IIF MODEL METHOD
There are a lot of methods for locating the fluid surface, for example, the reconstructed method adopted by Zhang (2010). In summary, CSF model is straightforward, but it needs to conduct the fluid surface location, calculate the surface curvature and use the surface tension parameters. This approach involves quite complex calculations of front curvatures that may lead to significant errors in some cases. It is not an easy work to conduct those calculations for SPH. IIF model implicitly simulates the surface tension through introducing the interaction force between particles, thus the disadvantages of CSF model are avoided and the advantage of SPH (relatively strong ability to incorporate complicated physical effects into the SPH formulations) is also exerted.
In general, in term of surface tension simulation with SPH, IIF model is therefore one of the main methods for future development. The introduction of certain form of interaction force is most simple and straightforward when applying the IIF model, which is one of the main development directions on simulating surface tension with SPH.
However, certain aspects are also needed to be improved for the IIF model based on the introduction of certain form of interaction force. The main shortcoming of the IIF model is: The introduced interaction force strength coefficient is determined by users. The artificial role in determining the interaction force strength coefficient must be eliminated, in order to make the improvements concerned. Thus, there are many researches needed to be done to achieve the goal.
This study was supported by a Grant-in-Aid for the China National Natural Science Foundation (No. 11162002) and Guangxi Zhuang Autonomous Region Department of Education (No. 2013YB224).