The case of a sudden fluid impact on a vertical wall (Peregrine, 2003) has been examinated on a geometrically simple set up. (Viccione et al., 2009) shown how such kind of phenomenon is strongly affected by fluid compressibility, especially during the first stages. A fluid mass, 0.50m high and 4.00m long, moving with an initial velocity v0 = 10m/s is discretized into a collection of 20.000 particles whith an interparticle distance d0 = 0.01m. The resulting mass is at a close distance to the vertical wall, so the impact process takes place after few timesteps (Fig.

Fig. 6. Initial conditions with fluid particles (blue dots) approaching the wall (green dots)

The following Fig. 7 shows the results in terms of pressure at different times.

6). Timestep is automatically adjusted to satisfy the Courant limit of stability.

Fig. 5. Comparison among different cell sizes.

0,10 0,20 0,30 0,40 0,50 0,60 0,70 0,80 0,90 1,00

 /s s,b

practically of no use whatsoever.

time and memory requirements.

surfaces flows.

**5. Applications of the SPH method**

**5.1 Slamming loads on a vertical structure**

Fig. 4.d

Fig. 4.a. refers to the case with no fixed cell structure. Here obviously the "Verlet list" procedure is highly beneficial, even though it appears that the size of the list must be carefully chosen, in order to fully exploit it effects. Figs. 4.b to 4.d show different computational times, depending on the cell size. Figure 5, gives a better insight of the results, showing the non-dimensional running cost trend s/s,b. As can be seen, minimum is achieved for a certain grid size.

78 Hydrodynamics – Optimizing Methods and Tools

0,70 0,75 0,80 0,85 0,90 0,95 1,00

> 0,40 0,45 0,50 0,55 0,60 0,65 0,70

0,15 0,20 0,25 0,30 0,35 0,40

0,15 0,20 0,25 0,30

 /s s,b

tot/tot,b

 /s s,b

tot/tot,b

 /s s,b

tot/tot,b

1,00 1,10 1,20 1,30 1,40 1,50

Relative partial and total computational time Domain not partitioned into cells

=rv/rc

Relative partial and total computational time Domain partitioned into vertical slices: Dx,cell=7.50m; Dy,cell=7.09m

1,00 1,10 1,20 1,30 1,40 1,50

=rv/rc

Relative partial and total computational time Domain partitioned in cells whose dimensions are: Dx,cell=0,50m; Dy,cell=0,50m

1,00 1,10 1,20 1,30 1,40 1,50

=rv/rc

Relative partial and total computational time Domain partitioned in cells whose dimensions are: Dx,cell=0,25m; Dy,cell=0,25m

1,00 1,10 1,20 1,30 1,40 1,50

=rv/rc

 /s s,b

tot/tot,b

Fig. 4.a. refers to the case with no fixed cell structure. Here obviously the "Verlet list" procedure is highly beneficial, even though it appears that the size of the list must be carefully chosen, in order to fully exploit it effects. Figs. 4.b to 4.d show different computational times, depending on the cell size. Figure 5, gives a better insight of the results, showing the non-dimensional running cost trend s/s,b.

Fig. 4.a

Fig. 4.b

Fig. 4.c

Fig. 4.d

Fig. 5. Comparison among different cell sizes.

It appears that while both the linked cell list and the Verlet list do relieve the computational time, the comparative advantage of the linked cell increases to the point where it is practically of no use whatsoever.

Later on, (Domínguez et. al., 2010) proposed an innovative searching algorithm based on a dynamic updating of the Verlet list yielding more satisfying results in term of computational time and memory requirements.

#### **5. Applications of the SPH method**

Smoothed Particle Hydrodynamics has been applied to a number of cases involving free surfaces flows.

#### **5.1 Slamming loads on a vertical structure**

The case of a sudden fluid impact on a vertical wall (Peregrine, 2003) has been examinated on a geometrically simple set up. (Viccione et al., 2009) shown how such kind of phenomenon is strongly affected by fluid compressibility, especially during the first stages. A fluid mass, 0.50m high and 4.00m long, moving with an initial velocity v0 = 10m/s is discretized into a collection of 20.000 particles whith an interparticle distance d0 = 0.01m. The resulting mass is at a close distance to the vertical wall, so the impact process takes place after few timesteps (Fig. 6). Timestep is automatically adjusted to satisfy the Courant limit of stability.

Fig. 6. Initial conditions with fluid particles (blue dots) approaching the wall (green dots) The following Fig. 7 shows the results in terms of pressure at different times.

Simulating Flows with SPH: Recent Developments and Applications 81

point for which the relative distance yields a pressure greater than the threshold value. Once reached such point, those neighbouring particles, previously fixed, are then set free to move. Runout velocity is instead controlled by handling the shear stress bed with the fixed bed.

PT3

PT1

PT2

Fig. 8. Spatial discretization. Red circles represent the area where local triggering is imposed.

Fig. 9. Neighbour particle destabilization. a) Particle "i" is approaching the neighbour particle "j". b) Despite the relative distance "|rij|" is decreased, particle "j" is still fixed because pij< plim. c) Particle "j" is set free to move because the pressure "pij" has reached the

Next Figures show three instants for each SPH based simulation, with the indication of the

t = 50 secs t = 100 secs t = 150 secs t = 100 secs t = 50 secs t = 150 secs

Fig. 10. PT1 Particle triggered zone, limit pressure plim = 300kgf /cm2(left side), plim = 200kgf

threshold value "plim".

/cm2 (right side), viscosity coefficient bed=0.1.

volume mobilized.

Fig. 7. Pressure contour as the impact progress takes place.

The rising and the following evolution of high pressure values is clearly evident. The order of magnitude is about 106 Pa, as it would be expected according to the Jokowski formula p = ρ C0 Δv, with v = v0 = 10m/s. After about 1/100 seconds most of the Jokowsky like pressure peak, generated by the sudden impact with the surface, disappeared, following that, the pressure starts building up again at a slower rate.

#### **5.2 Simulating triggering and evolution of debris-flows with SPH**

The capability into simulating debris-flow initiation and following movement with the Smoothed Particle Hydrodynamics is here investigated. The available domain taken from an existing slope, has been discretized with a reference distance being d0=2.5m and particles forming triangles as equilateral as possible. A single layer of moving particles has been laid on the upper part of the slope (blue region in Fig. 8).

Triggering is here settled randomly, releasing a particle located in the upper part of a slope, while all the remaining ones are initially frozen. Motion is then related to the achievement of a pressure threshold plim (Fig. 9). The resulting process is like a domino effect or a cascading failure. While some particles are moving, they may approach others initially still, to the 80 Hydrodynamics – Optimizing Methods and Tools

t = 0.0005sec t = 0.001sec

t = 0.002sec t = 0.005sec

t = 0.009sec t = 0.020sec

The rising and the following evolution of high pressure values is clearly evident. The order of magnitude is about 106 Pa, as it would be expected according to the Jokowski formula p = ρ C0 Δv, with v = v0 = 10m/s. After about 1/100 seconds most of the Jokowsky like pressure peak, generated by the sudden impact with the surface, disappeared, following

The capability into simulating debris-flow initiation and following movement with the Smoothed Particle Hydrodynamics is here investigated. The available domain taken from an existing slope, has been discretized with a reference distance being d0=2.5m and particles forming triangles as equilateral as possible. A single layer of moving particles has been laid

Triggering is here settled randomly, releasing a particle located in the upper part of a slope, while all the remaining ones are initially frozen. Motion is then related to the achievement of a pressure threshold plim (Fig. 9). The resulting process is like a domino effect or a cascading failure. While some particles are moving, they may approach others initially still, to the

Fig. 7. Pressure contour as the impact progress takes place.

that, the pressure starts building up again at a slower rate.

on the upper part of the slope (blue region in Fig. 8).

**5.2 Simulating triggering and evolution of debris-flows with SPH**

point for which the relative distance yields a pressure greater than the threshold value. Once reached such point, those neighbouring particles, previously fixed, are then set free to move. Runout velocity is instead controlled by handling the shear stress bed with the fixed bed.

Fig. 8. Spatial discretization. Red circles represent the area where local triggering is imposed.

Fig. 9. Neighbour particle destabilization. a) Particle "i" is approaching the neighbour particle "j". b) Despite the relative distance "|rij|" is decreased, particle "j" is still fixed because pij< plim. c) Particle "j" is set free to move because the pressure "pij" has reached the threshold value "plim".

Next Figures show three instants for each SPH based simulation, with the indication of the volume mobilized.

Fig. 10. PT1 Particle triggered zone, limit pressure plim = 300kgf /cm2(left side), plim = 200kgf /cm2 (right side), viscosity coefficient bed=0.1.

Simulating Flows with SPH: Recent Developments and Applications 83

Also, computation aspects related to the neighbourhood definition have been discussed. Field variables, such as particle velocity or density, have been evaluated by smoothing interpolation of the corresponding values over the nearest neighbour particles located inside a cut-off radius "rc" The generation of a neighbour list at each time step takes a considerable portion of CPU time. Straightforward determination of which particles are inside the interaction range requires the computation of all pair-wise distances, a procedure whose computational time

Lastly, applications of SPH in fluid hydrodynamics concerning wave slamming and propagation of debris flows have been discussed. These phenomena – involving complex geometries and

The work has been equally shared among the authors. Special thanks to the C.U.G.Ri. (University Centre for the Prediction and Prevention of Great Hazards), center, for allowing

Allen M.P. & Tildesley D.J. (1987). Computer Simulation of Liquids; *Clarendon Press*; Oxford. Belytschko T.; Lu Y.Y. & Gu L. (1994). Element-free Galerkin methods. *International Journal*

Belytschko, T. & Xiao, S. (2002). Stability analysis of particle methods with corrected derivatives, *Computers and Mathematics with Applications*, Vol. 43, pp. 329-350. Benz W. (1990). Smoothed Particle Hydrodynamics: a review, in numerical modellying of Non-Linear Stellar Pulsation: Problems and Prospects, Kluwer Academic, Boston. Blink J.A. & Hoover WG. (1985). Fragmentation of suddenly heated liquids, *Phys Rev A*; Vol.

Bonet, J.; Kulasegaram S.; Rodriguez-Paz M.X. & Profit M. (2004). Variational formulation

Chialvo A.A. & Debenedetti P.G. (1983). On the use of the Verlet neighbour list in molecular

Cleary, P.W. (1998). Modelling confined multi-material heat and mass flows using SPH,

Dilts G. A. (1999). Moving –least squares-particle hydrodynamics I, consistency and

Domínguez, J. M.; Crespo, A. J. C. ; Gómez-Gesteira, M. & Marongiu, J. C. (2011). Neighbour

Dymond, J. H. & Malhotra, R. (1988). The Tait equation: 100 years on, International Journal of Thermophysics, *Vol*. 9, No. 6, pp. 941-951, doi: 10.1007/bf01133262. Gingold, R.A. & Monaghan, J.J. (1977). Smoothed Particle hydrodynamics: theory and application to non-spherical stars. *Mon. Not. R. Astr. Soc.*, Vol. 181, pp. 375-389.

for the smooth particle hydrodynamics (SPH) simulation of fluid and solid problems, *Computer Methods in Applied Mechanics and Engineering*, Vol. 193, No. 12,

stability. *International Journal for Numerical Methods in Engineering*, Vol. 44, No. 8,

lists in smoothed particle hydrodynamics. International Journal for Numerical

would be of the order O(N2), and therefore unpractical for large domains.

**7. Acknowledgment**

32, No. 2, pp. 1027-1035.

pp. 1245–1257.

pp. 1115–1155.

**8. References**

rapidly-varied free surfaces - are of great importance in science and technology.

all the computations here presented on the Opteron quad processor machine.

*for Numerical methods in Engineering*, Vol. 37, pp. 229 - 256.

dynamics, *Comp. Ph. Comm*, Vol. 60, pp. 215-224.

Methods in Fluids, 66: n/a. doi: 10.1002/fld.2481.

*Appl. Math. Modelling*, Vol. 22, pp. 981–993.

Fig. 11. PT2 Particle triggered zone, limit pressure plim = 300kgf /cm2(left side), plim = 200kgf /cm2 (right side), viscosity coefficient bed=0.1.

Fig. 12. PT3 Particle triggered zone, limit pressure plim = 300kgf /cm2(left side), plim = 200kgf /cm2 (right side), viscosity coefficient bed=0.1.

As can been seen from the above Figures 10 to 12, by varying the location of the triggering area and the limit pressure plim, the condition of motion are quite different. More specifically, the mobilized area increases when the isotropic pressure plim decreases.

#### **6. Conclusion**

Recent theoretical developments and practical applications of the Smoothed Particle Hydrodynamics (SPH) method have been discussed, with specific concern to liquids. The main advantage is the capability of simulating the computational domain with large deformations and high discontinuities, bearing no numerical diffusion because advection terms are directly evaluated.

Recent achievements of SPH have been presented, concerning (1) numerical schemes for approximating Navier Stokes governing equations, (2) smoothing or kernel function properties needed to perform the function approximation to the Nth order, (3) restoring consistency of kernel and particle approximation, yielding the SPH approximation accuracy. Also, computation aspects related to the neighbourhood definition have been discussed. Field variables, such as particle velocity or density, have been evaluated by smoothing interpolation of the corresponding values over the nearest neighbour particles located inside a cut-off radius "rc" The generation of a neighbour list at each time step takes a considerable portion of CPU time. Straightforward determination of which particles are inside the interaction range requires the computation of all pair-wise distances, a procedure whose computational time would be of the order O(N2), and therefore unpractical for large domains.

Lastly, applications of SPH in fluid hydrodynamics concerning wave slamming and propagation of debris flows have been discussed. These phenomena – involving complex geometries and rapidly-varied free surfaces - are of great importance in science and technology.

#### **7. Acknowledgment**

The work has been equally shared among the authors. Special thanks to the C.U.G.Ri. (University Centre for the Prediction and Prevention of Great Hazards), center, for allowing all the computations here presented on the Opteron quad processor machine.

#### **8. References**

82 Hydrodynamics – Optimizing Methods and Tools

t = 50 secs t = 100 secs t = 150 secs t = 50 secs t = 100 secs t = 150 secs t

Fig. 11. PT2 Particle triggered zone, limit pressure plim = 300kgf /cm2(left side), plim = 200kgf

t = 50 secs t = 100 secs t = 150 secs t = 50 secs t = 100 secs t = 150 secs

Fig. 12. PT3 Particle triggered zone, limit pressure plim = 300kgf /cm2(left side), plim = 200kgf

As can been seen from the above Figures 10 to 12, by varying the location of the triggering area and the limit pressure plim, the condition of motion are quite different. More

Recent theoretical developments and practical applications of the Smoothed Particle Hydrodynamics (SPH) method have been discussed, with specific concern to liquids. The main advantage is the capability of simulating the computational domain with large deformations and high discontinuities, bearing no numerical diffusion because advection

Recent achievements of SPH have been presented, concerning (1) numerical schemes for approximating Navier Stokes governing equations, (2) smoothing or kernel function properties needed to perform the function approximation to the Nth order, (3) restoring consistency of kernel and particle approximation, yielding the SPH approximation accuracy.

specifically, the mobilized area increases when the isotropic pressure plim decreases.

**6. Conclusion**

terms are directly evaluated.

/cm2 (right side), viscosity coefficient bed=0.1.

/cm2 (right side), viscosity coefficient bed=0.1.

Allen M.P. & Tildesley D.J. (1987). Computer Simulation of Liquids; *Clarendon Press*; Oxford.

**5**

*Venezuela*

**3D Coalescence Collision of Liquid Drops Using**

The importance of modeling liquid drops collisions (see figure 1) is due to the existence of natural and engineering process where it is useful to understand the droplets dynamics in specific phenomena. Examples of applications are the combustion of fuel sprays, spray coating, emulsification, waste treatment and raindrop formation (Bozzano & Dente, 2010;

In this study we apply the Smoothed Particle Hydrodynamics method (SPH) to simulate for the first time the hydrodynamic collision of liquid drops on a vacuum environment in a three-dimensional space. When two drops collide a circular flat film is formed, and for sufficiently energetic collisions the evolution of the dynamics leads to a broken interface and to a bigger drop as a result of coalescence. We have shown that the SPH method can be useful to simulate in 3D this kind of process. As a result of the collision between the droplets the formation of a circular flat film is observed and depending on the approach velocity between the droplets different scenarios may arise: (i) if the film formed on the droplets collision is stable, then flocks of attached drops can appear; (ii) if the attractive interaction across the interfacial film is predominant, then the film is unstable and ruptures may occur leading to the formation of a bigger drop (permanent coalescence); (iii) under certain conditions the drops can rebound and the emulsion will be stable. Another possible scenario

Bradley & Stow, 1978;Park & Blair, 1975; Rourke & Bracco, 1980; Shah et al., 1972).

when two drops collide in a vacuum environment is the fragmentation of the drops.

Many studies has been proposed for the numerical simulation of the coalescence and break up of droplets (Azizi & Al Taweel, 2010; Cristini et al., 2001; Decent et al., 2006; Eggers et al., 1999; Foote ,1974; Jia et al., 2006; Mashayek et al., 2003; Narsimhan, 2004; Nobari et al., 1996; Pan & Suga, 2005; Roisman, 2004; Roisman et al., 2009; Sun et al., 2009; Xing et al., 2007; Yoon et al., 2007). In these studies, the authors propose different methods to approach the dynamics of liquid drops by a numerical integration of the Navier-Stokes equations. These examine the motion of droplets and the dynamics that it follows in time and study the liquid bridge that arises when two drops collide. The effects of parameters such as Reynolds number, impact velocity, drop size ratio and internal circulation are investigated and different regimes for droplets collisions are simulated. In some cases, those calculations yield results corresponding to four regimes of binary collisions: bouncing, coalescence, reflexive separation and stretching separation. These numerical simulations suggest that the collisions that lead to rebound between the drops are governed by macroscopic dynamics. In these simulations the mechanism of formation of satellite drops was also studied,

**1. Introduction**

**Smoothed Particle Hydrodynamics**

Alejandro Acevedo-Malavé and Máximo García-Sucre

*Venezuelan Institute for Scientific Research (IVIC)*

## **3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics**

Alejandro Acevedo-Malavé and Máximo García-Sucre *Venezuelan Institute for Scientific Research (IVIC) Venezuela*

#### **1. Introduction**

84 Hydrodynamics – Optimizing Methods and Tools

Krongauz Y. & Belytschko T. (1997). Consistent pseudo derivatives in meshless methods. *Computer methods in applied mechanics and engineering*, Vol. 146, pp. 371-386. Lee, E.S.; Violeau, D.; Issa, R. & Ploix, S. (2010) Application of weakly compressible and

Liu M. B.; Liu G. R. & Lam K. Y. (2003a). A one dimensional meshfree particle formulation for simulating shock waves, *Shock Waves*, Vol. 13, No. 3, pp.201 – 211. Liu M. B.; Liu G. R. & Lam K. Y. (2003b). Constructing smoothing functions in smoothed

Liu, G. R. & Liu, M. B. (2003). *Smoothed particle hydrodynamics: a meshfree particle method*,

Liu W. L.; Jun S., Li S. ; Adee J. & Belytschko T. (1995). Reproducing kernel particle methods

Chen, J. S.; Yoon, S.; Wang, H. P. & Liu, W. K. (2000). An Improved Reproducing Kernel

Lucy, L.B. (1977). A numerical approach to the testing of the fission hypothesis. *Astronomical*

Monaghan, J.J. (1988). Introduction to SPH. *Computer Physics Communication*, Vol.48, pp. 89 –96. Monaghan, J.J. (1994). Smoothed particle hydrodynamics, *Annual Review of Astronomy and*

Monaghan, J.J. & Gingold, R.A. (1983). Shock simulation by the particle method SPH, *Journ.*

Monaghan, J.J. & Lattanzio J.C. (1985). A refined particle method for astrophysical problems.

Oger, G.; Doring, M.; Alessandrini, B.; & Ferrant, P. (2007). An improved SPH method:

Peregrine, D.H. (2003). Water wave impact on walls. *Ann. Rev. Fluid Mech*, Vol. 35, pp. 23-43. Randles, P. W.; Libersky, L. D. & Petschek, A. G. (1999). On neighbors, derivatives, and viscosity in particle codes, *Proceedings of ECCM Conference*, Munich, Germany. Swegle, J. W.; Attaway, S. W.; Heinstein, M. W.; Mello, F. J. & Hicks, D. L. (1994). An analysis of smooth particle hydrodynamics. *Sandia Report* SAND93-2513. Verlet L. (1967). Computer Experiments on Classical Fluids. *Phys. Rev.* Vol. 159, No. 1, pp.

Vila, J.P. (1999). On particle weighted methods and smooth particle hydrodynamics. *Mathematical Models and Methods in Applied Sciences*, Vol.9, No.2, pp. 191–209. Viccione, G., Bovolin, V. & Carratelli, E. P. (2008). Defining and optimizing algorithms for

*for Numerical Methods in Fluids*, Vol. 58 pp. 625–638. doi: 10.1002/fld.1761. Viccione G.; Bovolin, V. & Pugliese Carratelli E. (2009). Influence of the compressibility in

neighbouring particle identification in SPH fluid simulations, *International Journal*

Fluid - Structure interaction Using Weakly Compressible SPH. 4rd ERCOFTAC

Towards higher order convergence. *Journal of Computational Physics*, Vol. 225, No.2,

*in Applied Mechanics and Engineering*, Vol. 181, No. 1-3, pp. 117-145.

*Research*, Vol. 48(Extra Issue), pp. 50–60, doi:10.3826/jhr.2010.0003.

*Mathematics*, Vol. 155, No. 2, pp. 263-284.

Vol. 38, pp. 1655-1679.

pp. 1472-1492.

98-103.

*Journal*, Vol. 82, pp. 1013–1024.

*Astrophysics*, Vol. 30, pp. 543–574.

*of Comp. Physics*, Vol. 82, pp. 374-389.

*Astronomy and Astrophysics*, Vol. 149, pp. 135-143.

SPHERIC workshop on SPH applications. Nantes.

World Scientific, ISBN 981-238-456-1, Singapore.

truly incompressible SPH to 3-D water collapse in waterworks, *Journal of Hydraulic*

particle hydrodynamics with applications, *Journal of Computational and Applied*

for structural dynamics. *International Journal for Numerical Methods in Engineering*,

Particle Method for Nearly Incompressible Hyperelastic Solids, *Computer Methods*

The importance of modeling liquid drops collisions (see figure 1) is due to the existence of natural and engineering process where it is useful to understand the droplets dynamics in specific phenomena. Examples of applications are the combustion of fuel sprays, spray coating, emulsification, waste treatment and raindrop formation (Bozzano & Dente, 2010; Bradley & Stow, 1978;Park & Blair, 1975; Rourke & Bracco, 1980; Shah et al., 1972).

In this study we apply the Smoothed Particle Hydrodynamics method (SPH) to simulate for the first time the hydrodynamic collision of liquid drops on a vacuum environment in a three-dimensional space. When two drops collide a circular flat film is formed, and for sufficiently energetic collisions the evolution of the dynamics leads to a broken interface and to a bigger drop as a result of coalescence. We have shown that the SPH method can be useful to simulate in 3D this kind of process. As a result of the collision between the droplets the formation of a circular flat film is observed and depending on the approach velocity between the droplets different scenarios may arise: (i) if the film formed on the droplets collision is stable, then flocks of attached drops can appear; (ii) if the attractive interaction across the interfacial film is predominant, then the film is unstable and ruptures may occur leading to the formation of a bigger drop (permanent coalescence); (iii) under certain conditions the drops can rebound and the emulsion will be stable. Another possible scenario when two drops collide in a vacuum environment is the fragmentation of the drops.

Many studies has been proposed for the numerical simulation of the coalescence and break up of droplets (Azizi & Al Taweel, 2010; Cristini et al., 2001; Decent et al., 2006; Eggers et al., 1999; Foote ,1974; Jia et al., 2006; Mashayek et al., 2003; Narsimhan, 2004; Nobari et al., 1996; Pan & Suga, 2005; Roisman, 2004; Roisman et al., 2009; Sun et al., 2009; Xing et al., 2007; Yoon et al., 2007). In these studies, the authors propose different methods to approach the dynamics of liquid drops by a numerical integration of the Navier-Stokes equations. These examine the motion of droplets and the dynamics that it follows in time and study the liquid bridge that arises when two drops collide. The effects of parameters such as Reynolds number, impact velocity, drop size ratio and internal circulation are investigated and different regimes for droplets collisions are simulated. In some cases, those calculations yield results corresponding to four regimes of binary collisions: bouncing, coalescence, reflexive separation and stretching separation. These numerical simulations suggest that the collisions that lead to rebound between the drops are governed by macroscopic dynamics. In these simulations the mechanism of formation of satellite drops was also studied,

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 87

Tartakovsky & Meakin (2005) have shown that the artificial surface tension that emerge from the standard formulation of the Smoothed Particle Hydrodynamics (SPH) method (Gingold & Monaghan, 1977) could be eliminated by using SPH equations based on the number density of particles instead of the density of particles in the fluid. The contribution of Tartakovsky & Meakin (2005) could be very useful when modeling the hydrodynamic interaction of drops in liquid emulsions. Combining these schemes with some continuousdiscrete hybrid approach (Cui et al., 2006; Koumoutsakos, 2005; Li et al., 1998; Nie et al., 2004; O'Connell & Thompson, 1995) it could be constructed an interesting model to discuss the collapse and disappearance of the interfacial film in emulsion media (Bibette et al., 1992; Ivanov & Dimitrov, 1988; Ivanov & Kralchevsky, 1997; Kabalnov & Wennerström, 1996; Sharma & Ruckenstein, 1987). Ivanov & Kralchevsky (1997) conducted a study on the possible outcomes for the collision of liquid droplets in emulsions. According to this study, when the collision between two drops occurs, an interfacial film of flat circular section is formed, and coalescence or flocculation may arise (Ivanov & Kralchevsky, 1997). These authors did not carry out the hydrodynamical modeling of collision between drops. Instead, they discuss thermodynamics and hydrodynamics aspects of the problem and raise some

In this work we apply the SPH method to simulate for the first time in three-dimensional space the hydrodynamic coalescence collision of liquid drops in a vacuum environment. This method is employed in order to obtaining approximate numerical solutions of the equations of fluid dynamics by replacing the fluid with a set of particles. These particles may be interpreted as corresponding to interpolation points from which properties of the fluid can be determined. Each SPH particle can be considered as a system of smaller particles. The SPH method is particularly useful when the fluid motion produce big

All our calculations were performed defining inside the SPH code two drops composed by 4700 SPH particles, running on a Dell Work Station with 8 processors Intel Xeon of 3.33 Ghz

The SPH method was invented first and simultaneously by Lucy, (1977) and Gingold & Monaghan (1977) to solve astrophysical problems. This method has been used to study a range of astrophysical topics including formation of galaxies, formation of stars, supernovas, stellar collisions, and so on. This method has the advantage that if you want to model more than one material, the interface problems arising can be modeled easily, while they are hard to model using other methods based on finite differences. An additional advantage is that SPH method can be considered as a bridge between continuous and fragmented material, which makes it one of the best method to study problems of fragmentation in solids (Benz & Asphaug, 1994, 1995). Another feature that makes the SPH method attractive is that it yields solutions depending on space and time, making it versatile for treating a wide variety of problems in physics. Furthermore, given the similarity between SPH and molecular dynamics, combination of these two methods can be used to treat complex problems in systems that differ considerably in their length scales. The easiness of the method to be adaptable and their Lagrangian character make of SPH one of the most popular among existing numerical methods used for modeling fluids. On the other hand, the SPH method can be used to describe the dynamics of deformable bodies (Desbrun

possible outcomes when two liquid droplets collide.

deformations and a large velocity of the whole fluid.

**2. Smoothed particle hydrodynamics method**

with 32.0 GB of RAM memory.

confirming that the principal cause of the formation of satellite drops is the "end pinching" while the capillary wave instabilities are the dominant feature in cases where a large value of the parameter impact is employed.

Experimental studies on the coalescence process involving the production of satellite droplets has been reported in the literature (Ashgriz & Givi, 1987, 1989; Brenn & Frohn, 1989; Brenn & Kolobaric, 2006; Zhang et al., 2009). These authors found out that when the Weber number increases, the collision takes the form of a high-energy one and results of different type may arise. In these references the results show that the collision of the droplets can be bouncing, grazing and generating satellite drops. Based on data from experiments on the formation and breaking up of ligaments, the process of satellite droplets formation is modeled by these authors and the experiments are carried out using various liquid streams. On the other hand, for Weber numbers corresponding to a high-energy collision, permanent coalescence occurs and the bigger drop is deformed producing satellite drops. Experimental studies on the binary collision of droplets for a wide range of Weber numbers and impact parameters have been carried out and reported in the literature (Ashgriz & Poo, 1990; Gotaas et al., 2007b; Menchaca-Rocha et al., 1997; Qian & Law, 1997). These authors identified two types of collisions leading to drops separation, which can be reflexive or stretching separation. It was found that the reflexive separation occurs for headon collisions, while stretching separation occurs for high values of the impact parameter. Carrying out Experiments, the authors reported the transition between two types of separation, and also collisions that lead to coalescence. In these references experimental investigations of the transition between different regimes of collisions were reported. The authors analyzed the results using photographic images, which showed the evolution of the dynamics exhibited by the droplets. As a result of these experiments were proposed five different regimes governing the collision between droplets: (i) coalescence after a small deformation, (ii) bouncing, (iii) coalescence after substantial deformation, (iv) coalescence followed by separation for head-on collisions, and (v) coalescence followed by separation for off-center collisions.

Li (1994) and Chen (1985) studied the coalescence of two small bubbles or drops using a model for the dynamics of the thinning film in which both, London-van der Waals and electrostatic double layer forces, are taken into account. Li (1994) proposes a general expression for the coalescence time in the absence of the electrostatic double layer forces. The model proposed by Chen (1985), depending on the radius of the drops and the physical properties of the fluids and surfaces, describes the film profile evolution and predicts the film stability, time scale and film thickness.

The dynamics of collision between equal-sized liquid drops of organic substances has also been reported in the literature (Ashgriz & Givi, 1987, 1989; Gotaas et al., 2007a; Jiang et al., 1992; Podgorska, 2007). They reported the experimental results of the collision of water and normal-alkane droplets in the radius range of 150 m. These results showed that for the studied range of Weber numbers, the behavior of hydrocarbon droplets is more complex than the observed for water droplets. For water droplets head-on collisions, permanent coalescence always result. Experimental studies on the different ways in which may occur the coalescence of drops, have been performed by different authors (Gokhale et al., 2004; Leal, 2004; Menchaca-Rocha et al., 2001; Mohamed-Kassim & Longmire, 2004; Thoroddsen et al., 2007; Wang et al., 2009; Wu et al., 2004). In these studies are reported the evolution in time of the surface shape as well as a broad view of the contact region between the droplets.

86 Hydrodynamics – Optimizing Methods and Tools

confirming that the principal cause of the formation of satellite drops is the "end pinching" while the capillary wave instabilities are the dominant feature in cases where a large value

Experimental studies on the coalescence process involving the production of satellite droplets has been reported in the literature (Ashgriz & Givi, 1987, 1989; Brenn & Frohn, 1989; Brenn & Kolobaric, 2006; Zhang et al., 2009). These authors found out that when the Weber number increases, the collision takes the form of a high-energy one and results of different type may arise. In these references the results show that the collision of the droplets can be bouncing, grazing and generating satellite drops. Based on data from experiments on the formation and breaking up of ligaments, the process of satellite droplets formation is modeled by these authors and the experiments are carried out using various liquid streams. On the other hand, for Weber numbers corresponding to a high-energy collision, permanent coalescence occurs and the bigger drop is deformed producing satellite drops. Experimental studies on the binary collision of droplets for a wide range of Weber numbers and impact parameters have been carried out and reported in the literature (Ashgriz & Poo, 1990; Gotaas et al., 2007b; Menchaca-Rocha et al., 1997; Qian & Law, 1997). These authors identified two types of collisions leading to drops separation, which can be reflexive or stretching separation. It was found that the reflexive separation occurs for headon collisions, while stretching separation occurs for high values of the impact parameter. Carrying out Experiments, the authors reported the transition between two types of separation, and also collisions that lead to coalescence. In these references experimental investigations of the transition between different regimes of collisions were reported. The authors analyzed the results using photographic images, which showed the evolution of the dynamics exhibited by the droplets. As a result of these experiments were proposed five different regimes governing the collision between droplets: (i) coalescence after a small deformation, (ii) bouncing, (iii) coalescence after substantial deformation, (iv) coalescence followed by separation for head-on collisions, and (v) coalescence followed by separation

Li (1994) and Chen (1985) studied the coalescence of two small bubbles or drops using a model for the dynamics of the thinning film in which both, London-van der Waals and electrostatic double layer forces, are taken into account. Li (1994) proposes a general expression for the coalescence time in the absence of the electrostatic double layer forces. The model proposed by Chen (1985), depending on the radius of the drops and the physical properties of the fluids and surfaces, describes the film profile evolution and predicts the

The dynamics of collision between equal-sized liquid drops of organic substances has also been reported in the literature (Ashgriz & Givi, 1987, 1989; Gotaas et al., 2007a; Jiang et al., 1992; Podgorska, 2007). They reported the experimental results of the collision of water and normal-alkane droplets in the radius range of 150 m. These results showed that for the studied range of Weber numbers, the behavior of hydrocarbon droplets is more complex than the observed for water droplets. For water droplets head-on collisions, permanent coalescence always result. Experimental studies on the different ways in which may occur the coalescence of drops, have been performed by different authors (Gokhale et al., 2004; Leal, 2004; Menchaca-Rocha et al., 2001; Mohamed-Kassim & Longmire, 2004; Thoroddsen et al., 2007; Wang et al., 2009; Wu et al., 2004). In these studies are reported the evolution in time of the surface shape as well as a broad view of the contact region

of the parameter impact is employed.

for off-center collisions.

between the droplets.

film stability, time scale and film thickness.

Tartakovsky & Meakin (2005) have shown that the artificial surface tension that emerge from the standard formulation of the Smoothed Particle Hydrodynamics (SPH) method (Gingold & Monaghan, 1977) could be eliminated by using SPH equations based on the number density of particles instead of the density of particles in the fluid. The contribution of Tartakovsky & Meakin (2005) could be very useful when modeling the hydrodynamic interaction of drops in liquid emulsions. Combining these schemes with some continuousdiscrete hybrid approach (Cui et al., 2006; Koumoutsakos, 2005; Li et al., 1998; Nie et al., 2004; O'Connell & Thompson, 1995) it could be constructed an interesting model to discuss the collapse and disappearance of the interfacial film in emulsion media (Bibette et al., 1992; Ivanov & Dimitrov, 1988; Ivanov & Kralchevsky, 1997; Kabalnov & Wennerström, 1996; Sharma & Ruckenstein, 1987). Ivanov & Kralchevsky (1997) conducted a study on the possible outcomes for the collision of liquid droplets in emulsions. According to this study, when the collision between two drops occurs, an interfacial film of flat circular section is formed, and coalescence or flocculation may arise (Ivanov & Kralchevsky, 1997). These authors did not carry out the hydrodynamical modeling of collision between drops. Instead, they discuss thermodynamics and hydrodynamics aspects of the problem and raise some possible outcomes when two liquid droplets collide.

In this work we apply the SPH method to simulate for the first time in three-dimensional space the hydrodynamic coalescence collision of liquid drops in a vacuum environment. This method is employed in order to obtaining approximate numerical solutions of the equations of fluid dynamics by replacing the fluid with a set of particles. These particles may be interpreted as corresponding to interpolation points from which properties of the fluid can be determined. Each SPH particle can be considered as a system of smaller particles. The SPH method is particularly useful when the fluid motion produce big deformations and a large velocity of the whole fluid.

All our calculations were performed defining inside the SPH code two drops composed by 4700 SPH particles, running on a Dell Work Station with 8 processors Intel Xeon of 3.33 Ghz with 32.0 GB of RAM memory.

#### **2. Smoothed particle hydrodynamics method**

The SPH method was invented first and simultaneously by Lucy, (1977) and Gingold & Monaghan (1977) to solve astrophysical problems. This method has been used to study a range of astrophysical topics including formation of galaxies, formation of stars, supernovas, stellar collisions, and so on. This method has the advantage that if you want to model more than one material, the interface problems arising can be modeled easily, while they are hard to model using other methods based on finite differences. An additional advantage is that SPH method can be considered as a bridge between continuous and fragmented material, which makes it one of the best method to study problems of fragmentation in solids (Benz & Asphaug, 1994, 1995). Another feature that makes the SPH method attractive is that it yields solutions depending on space and time, making it versatile for treating a wide variety of problems in physics. Furthermore, given the similarity between SPH and molecular dynamics, combination of these two methods can be used to treat complex problems in systems that differ considerably in their length scales. The easiness of the method to be adaptable and their Lagrangian character make of SPH one of the most popular among existing numerical methods used for modeling fluids. On the other hand, the SPH method can be used to describe the dynamics of deformable bodies (Desbrun

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 89

*W*(**r**

The SPH discretization reduces the Navier-Stokes equation to a system of ordinary differential equations having the form of Newton's second law of motion for each particle. This simplicity allows taking into account a variety of chemical effects with relatively little effort in the development of computational codes. Also, since the number of particles remains constant in the simulation and the interactions are symmetrical, the mass, momentum and energy are conserved exactly, and the systems like dynamic boundaries and interfaces can be modeled without too much difficulty. Hoover (1998), and Colagrossi & Landrini (2003), used the SPH method to model immiscible flows and found that the standard formulation of SPH proposed by Gingold & Monaghan (1977) creates an artificial surface tension on the border between the two fluids. Colagrossi & Landrini (2003) put forward an SPH formulation for the simulation of interfacial flows, that is, flow fields of different fluids separated by interfaces. The scheme proposed for the simulation of interfacial flows starts considering that the fluid field is represented by a collection of N particles interacting with each other according to evolution equations of the general form

*<sup>i</sup>* **<sup>r</sup>** *<sup>j</sup>* ,*h*)

*<sup>i</sup>* **<sup>r</sup>** *<sup>j</sup>* ,*h*)

. (7)

. (8)

(9)

(10)

*<sup>i</sup> <sup>m</sup> <sup>j</sup>*

*Ai*

*d i dt <sup>i</sup> Mij <sup>j</sup>* ,

*d***u***i dt* <sup>1</sup> *i* F*ij j f i* ,

*d***x***i dt* **<sup>u</sup>***<sup>i</sup>* .

difficulty can be circumvented using the following discrete approximations

*div*(**u***<sup>i</sup>* ) (**<sup>u</sup>** *<sup>j</sup>* **<sup>u</sup>***<sup>i</sup>* ) *j*

*Ai* (*<sup>A</sup> <sup>j</sup> Ai* ) *j*

The terms Mij and Fij arise from the mass and momentum conservation equations. In the equations (9) appear the density i, the velocity ui of the particles, and the force fi can be any body force. When there are fluid regions with a sharp density gradient (interfaces), the SPH standard formulations must be modified in order to be applied to treat such systems. This

*<sup>W</sup> ji*

*<sup>W</sup> ji*

Here W is the Kernel or Smoothing Function and A can be any scalar field or continuous function. The small difference between the equation (10) and the standard equation that uses

*m j j* ,

*m j j* .

Similarly, the gradient can be calculated using the expression

*j*

*A j n j i W* (**r**

*j*

The mass density is given by

& Gascuel, 1996). Currently there are several applications of SPH in different areas related to fluid dynamics, such as: incompressible flows, elastic flows, multiphase flows, supersonic flows, shock wave simulation, heat transfer, explosive phenomena, and so on (Liu & Liu, 2003; Monaghan, 1992). A major advantage of SPH is that their physical interpretation is relatively simple.

In the SPH model, the fluid is represented by a discrete set of N particles. The position of the ith particle is denoted by the vector ri, i=1,…, N. We start introducing the function As(r), that is the smoothed representation of any arbitrary function A(r) (the function A(r) is any physical quantity of the hydrodynamical model and As(r) is the smoothed version of this quantity). The SPH scheme is based on the idea of a smoothed representation As(r) of the continuous function A(r) that can be obtained from the convolution integral

$$A\_S(\mathbf{r}) = \left[ A(\mathbf{r}')W(\mathbf{r} - \mathbf{r}', h)d\mathbf{r}'.\right].\tag{1}$$

Here h is the smoothing length, and the smoothing function W satisfies the normalization condition

$$\int W(\mathbf{r} - \mathbf{r}', h)d\mathbf{r}' = 1.\tag{2}$$

The integration is performed over the whole space. In the limit of h tending to zero, the smoothing function W becomes a Dirac delta function, and the smoothed representation As(r) tends to A(r).

In the SPH scheme, the properties associated with particle i, are calculated by approximating the integral in eq. (1) by the sum

$$\begin{aligned} A\_{\bar{i}} &= \sum\_{\bar{j}} \Delta V\_{\bar{j}} A\_{\bar{j}} W(\mathbf{r}\_{\bar{i}} - \mathbf{r}\_{\bar{j}}, h) \\ &= \sum\_{\bar{j}} m\_{\bar{j}} \frac{A\_{\bar{j}}}{\rho\_{\bar{j}}} W(\mathbf{r}\_{\bar{i}} - \mathbf{r}\_{\bar{j}}, h). \end{aligned} \tag{3}$$

Here ∆Vj is the fluid volume associated with particle j, and mj and j are the mass and density of the jth particle, respectively. In equation (3), Aj is the value of a physical field A(r) on the particle j, and the sum is performed over all particles. Furthermore, the gradient of A is calculated using the expression

$$
\nabla A\_{\dot{i}} = \sum\_{\dot{j}} m\_{\dot{j}} \frac{A\_{\dot{j}}}{\rho\_{\dot{j}}} \nabla\_{\dot{j}} W(\mathbf{r}\_{\dot{i}} - \mathbf{r}\_{\dot{j}}, h). \tag{4}
$$

In the equation (3), i/mi can be replaced by the particle number density ni=i/mi, so that

$$A\_{\dot{l}} = \sum\_{\dot{j}} \frac{A\_{\dot{j}}}{n\_{\dot{j}}} W(\mathbf{r}\_{\dot{l}} - \mathbf{r}\_{\dot{j}}, h). \tag{5}$$

The particle number density can be calculated using the expression

$$m\_{\vec{i}} = \sum\_{\vec{j}} W(\mathbf{r}\_{\vec{i}} - \mathbf{r}\_{\vec{j}}, h). \tag{6}$$

The mass density is given by

88 Hydrodynamics – Optimizing Methods and Tools

& Gascuel, 1996). Currently there are several applications of SPH in different areas related to fluid dynamics, such as: incompressible flows, elastic flows, multiphase flows, supersonic flows, shock wave simulation, heat transfer, explosive phenomena, and so on (Liu & Liu, 2003; Monaghan, 1992). A major advantage of SPH is that their physical interpretation is

In the SPH model, the fluid is represented by a discrete set of N particles. The position of the ith particle is denoted by the vector ri, i=1,…, N. We start introducing the function As(r), that is the smoothed representation of any arbitrary function A(r) (the function A(r) is any physical quantity of the hydrodynamical model and As(r) is the smoothed version of this quantity). The SPH scheme is based on the idea of a smoothed representation As(r) of the

Here h is the smoothing length, and the smoothing function W satisfies the normalization

The integration is performed over the whole space. In the limit of h tending to zero, the smoothing function W becomes a Dirac delta function, and the smoothed representation

In the SPH scheme, the properties associated with particle i, are calculated by

*A j W*(**r** *<sup>i</sup>* **<sup>r</sup>** *<sup>j</sup>* ,*h*)

Here ∆Vj is the fluid volume associated with particle j, and mj and j are the mass and density of the jth particle, respectively. In equation (3), Aj is the value of a physical field A(r) on the particle j, and the sum is performed over all particles. Furthermore, the gradient of A

> *A j j j <sup>i</sup>*

In the equation (3), i/mi can be replaced by the particle number density ni=i/mi, so that

*A j n j W*(**r**

*j*

*ni <sup>W</sup>* (**<sup>r</sup>**

*j*

*<sup>i</sup>* **<sup>r</sup>** *<sup>j</sup>* ,*h*).

*W*(**r**

*<sup>i</sup>* **<sup>r</sup>** *<sup>j</sup>* ,*h*).

*<sup>i</sup>* **<sup>r</sup>** *<sup>j</sup>* ,*h*).

(**r**) *<sup>A</sup>*(**r**)*<sup>W</sup>* (**<sup>r</sup> <sup>r</sup>**, *<sup>h</sup>*)*d***r**. *<sup>s</sup> <sup>A</sup>* (1)

( , ) 1. *W* **r r** *h d***r** (2)

(3)

*<sup>i</sup>* **<sup>r</sup>** *<sup>j</sup>*,*h*). (4)

(5)

(6)

continuous function A(r) that can be obtained from the convolution integral

*Ai <sup>V</sup> <sup>j</sup>*

*j*

*A j j j W* (**r**

*<sup>m</sup> <sup>j</sup>*

*Ai <sup>m</sup> <sup>j</sup>*

*Ai*

The particle number density can be calculated using the expression

relatively simple.

condition

As(r) tends to A(r).

approximating the integral in eq. (1) by the sum

is calculated using the expression

$$\rho\_{\dot{i}} = \sum\_{\dot{j}} m\_{\,\,\dot{j}} W(\mathbf{r}\_{\dot{i}} - \mathbf{r}\_{\dot{j}}, h). \tag{7}$$

Similarly, the gradient can be calculated using the expression

$$
\nabla A\_{\dot{i}} = \sum\_{j} \frac{A\_{\dot{j}}}{n\_{\dot{j}}} \nabla\_{\dot{j}} W(\mathbf{r}\_{\dot{i}} - \mathbf{r}\_{\dot{j}}, h). \tag{8}
$$

The SPH discretization reduces the Navier-Stokes equation to a system of ordinary differential equations having the form of Newton's second law of motion for each particle. This simplicity allows taking into account a variety of chemical effects with relatively little effort in the development of computational codes. Also, since the number of particles remains constant in the simulation and the interactions are symmetrical, the mass, momentum and energy are conserved exactly, and the systems like dynamic boundaries and interfaces can be modeled without too much difficulty. Hoover (1998), and Colagrossi & Landrini (2003), used the SPH method to model immiscible flows and found that the standard formulation of SPH proposed by Gingold & Monaghan (1977) creates an artificial surface tension on the border between the two fluids. Colagrossi & Landrini (2003) put forward an SPH formulation for the simulation of interfacial flows, that is, flow fields of different fluids separated by interfaces. The scheme proposed for the simulation of interfacial flows starts considering that the fluid field is represented by a collection of N particles interacting with each other according to evolution equations of the general form

$$\begin{aligned} \frac{d\rho\_{\bar{i}}}{dt} &= -\rho\_{\bar{i}} \sum\_{\bar{j}} M\_{\bar{i}\bar{j}}, \\ \frac{d\mathbf{u}\_{\bar{i}}}{dt} &= -\frac{1}{\rho\_{\bar{i}}} \sum\_{\bar{j}} \mathbf{F}\_{\bar{i}\bar{j}} + f\_{\bar{i}}, \\ \frac{d\mathbf{x}\_{\bar{i}}}{dt} &= \mathbf{u}\_{\bar{i}}. \end{aligned} \tag{9}$$

The terms Mij and Fij arise from the mass and momentum conservation equations. In the equations (9) appear the density i, the velocity ui of the particles, and the force fi can be any body force. When there are fluid regions with a sharp density gradient (interfaces), the SPH standard formulations must be modified in order to be applied to treat such systems. This difficulty can be circumvented using the following discrete approximations

$$\begin{split}d\boldsymbol{\dot{m}}(\mathbf{u}\_{\dot{I}}) &= \boldsymbol{\Sigma}(\mathbf{u}\_{\ \dot{J}} - \mathbf{u}\_{\ \dot{I}}) \cdot \nabla W\_{j\dot{I}} \frac{m\_{\ \ \dot{J}}}{\rho\_{\ \ \dot{J}}}, \\ \nabla A\_{\dot{I}} &= \boldsymbol{\Sigma}(A\_{\ \dot{J}} - A\_{\dot{I}}) \nabla W\_{j\dot{I}} \frac{m\_{\ \ \dot{J}}}{\rho\_{\ \ \dot{J}}}.\end{split} \tag{10}$$

Here W is the Kernel or Smoothing Function and A can be any scalar field or continuous function. The small difference between the equation (10) and the standard equation that uses

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 91

In this formulation special attention must be paid to the kernel. In fact depending on which kernel is used, eq. (13) could introduce additional errors. For this reason a first-order interpolation scheme is suitable to re-initialize the density field by using the equation

The XSPH (Extended Smoothed Particle Hydrodynamics, which is a variant of the SPH method for the modeling of free surface flows (Monaghan, 1994)) velocity correction ui is introduced to prevent particles inter-penetration (Colagrossi & Landrini, 2003), which takes into account the velocity of the neighbor particles using a mean value of the velocity,

> '

*m j j ij*

> *j* 2

*xi i* 2 *i i i*

*h* ,

2

*ij* is the mean value of density between the ith and jth particle, and ' is the relative

*W j*

*j*

*MLS* (*xi* )

, (14)

(*<sup>u</sup> <sup>j</sup> ui* )*<sup>W</sup> ji*, (15)

(16)

, (17)

<sup>0</sup> (18)

*MLS* (*xi* )

*dV <sup>j</sup> <sup>m</sup> <sup>j</sup>*

, *ui*

change of an arbitrary quantity between simulations (Colagrossi & Landrini, 2003).

 *i*

<sup>v</sup> *<sup>j</sup>*

 *Wij*

In the above equation p is the pressure, is the dynamic viscosity and is the shear strain

In the present work, our calculations are performed in three dimensions and we use the cubic B-spline kernel (Monaghan, 1985). We consider water drops, and the equation of state that we use in the hydrodynamical code was a general Mie-Gruneisen form of equation of state with different analytic forms for states of compression (/0-1)>0 and tension (/0- 1)<0 (Liu & Liu, 2003). This equation has several parameters, namely the density , the

 *i* <sup>2</sup> *j*

*j* 1

*Wij*

*<sup>i</sup>*

Where *MLS Wj* is the moving-least-square kernel.

according to the equations

where is the total stress tensor.

*dEi dt* <sup>1</sup>

*<sup>P</sup> <sup>A</sup>*<sup>1</sup>

 *<sup>A</sup>*<sup>2</sup>

 0 1

where

rate.

 *j W j*

*j*

**<sup>u</sup>***<sup>i</sup> ui ui*

The velocity and acceleration fields are (Liu & Liu, 2003)

*d***r** *i dt* **<sup>v</sup>***<sup>i</sup>* ,

*d***v***i*

<sup>2</sup> *<sup>m</sup> <sup>j</sup>*

*N*

*dt <sup>m</sup> <sup>j</sup>*

*N*

The internal energy evolution is given by the expression (Liu & Liu, 2003) :

*pi i* 2 *p j*

vi

reference density 0, and the constants A1, A2, A3, C1 and C2 . The pressure P is

 2

*<sup>A</sup>*<sup>3</sup>

 3

*if*

 0 1

 0 1

 0 1

*<sup>j</sup>* <sup>1</sup>

 *j* 2

mj/i instead mj/j is important for the treatment of the case of small density ratios. On the other hand, it can be shown that the pressure gradient can be written as

$$
\nabla p\_{\dot{j}} = \sum\_{\dot{j}} (p\_{\ j} + p\_{\dot{j}}) \nabla W\_{j\dot{i}} dV\_{\dot{j}}.\tag{11}
$$

The equation (11) is variationally consistent with eq. (10). In this scheme the terms Mij and Fij appearing in eq. (9) are given by the expressions

Fig. 1. Definition of the problem: head-on coalescence collision in three dimensions between two drops of equal size approaching with a velocity of collision Vcol and radius R in empty space. Each drop is composed by 4700 SPH particles.

A density re-initialization is needed when each particle has a fixed mass, and when the number of particles is constant the mass conservation is satisfied. Yet if one uses eq. (9) for the density, the consistency between mass, density and occupied area is not satisfied. To solve this problem, the density is periodically re-initialized applying the expression

$$
\rho\_{\dot{l}} = \sum\_{\dot{j}} m\_{\dot{j}} W\_{\dot{l}\dot{j}}.\tag{13}
$$

In this formulation special attention must be paid to the kernel. In fact depending on which kernel is used, eq. (13) could introduce additional errors. For this reason a first-order interpolation scheme is suitable to re-initialize the density field by using the equation

$$\left\langle \rho\_{\dot{I}} \right\rangle = \sum\_{\dot{J}} \rho\_{\dot{J}} W\_{\dot{J}}^{MLS}(\mathbf{x}\_{\dot{I}}) dV\_{\dot{J}} = \sum\_{\dot{J}} m\_{\dot{J}} W\_{\dot{J}}^{MLS}(\mathbf{x}\_{\dot{I}}),\tag{14}$$

Where *MLS Wj* is the moving-least-square kernel.

90 Hydrodynamics – Optimizing Methods and Tools

mj/i instead mj/j is important for the treatment of the case of small density ratios. On the

The equation (11) is variationally consistent with eq. (10). In this scheme the terms Mij and Fij

*m j j* ,

*m j j* .

*Mij* (*<sup>u</sup> <sup>j</sup> ui* ) *<sup>W</sup> ji*

*Fij* ( *<sup>p</sup> <sup>j</sup> pi* )*<sup>W</sup> ji*

Fig. 1. Definition of the problem: head-on coalescence collision in three dimensions between two drops of equal size approaching with a velocity of collision Vcol and radius R in empty

A density re-initialization is needed when each particle has a fixed mass, and when the number of particles is constant the mass conservation is satisfied. Yet if one uses eq. (9) for the density, the consistency between mass, density and occupied area is not satisfied. To

. (13)

solve this problem, the density is periodically re-initialized applying the expression

 *<sup>i</sup> <sup>m</sup> <sup>j</sup> Wij <sup>j</sup>*

space. Each drop is composed by 4700 SPH particles.

*<sup>W</sup> jidV <sup>j</sup>* . (11)

(12)

*pi* (*<sup>p</sup> <sup>j</sup> pi* ) *j*

other hand, it can be shown that the pressure gradient can be written as

appearing in eq. (9) are given by the expressions

The XSPH (Extended Smoothed Particle Hydrodynamics, which is a variant of the SPH method for the modeling of free surface flows (Monaghan, 1994)) velocity correction ui is introduced to prevent particles inter-penetration (Colagrossi & Landrini, 2003), which takes into account the velocity of the neighbor particles using a mean value of the velocity, according to the equations

$$
\left< \mathbf{u}\_{\dot{i}} \right> = \boldsymbol{u}\_{\dot{i}} + \Delta \boldsymbol{u}\_{\dot{i}}, \quad \Delta \boldsymbol{u}\_{\dot{i}} = \frac{\varepsilon'}{2} \sum\_{\dot{j}} \frac{m}{\overline{\rho}\_{\dot{i}\dot{j}}} (\boldsymbol{u}\_{\\_j} - \boldsymbol{u}\_{\\_i}) \boldsymbol{W}\_{j\dot{i}}, \tag{15}
$$

where *ij* is the mean value of density between the ith and jth particle, and ' is the relative change of an arbitrary quantity between simulations (Colagrossi & Landrini, 2003). The velocity and acceleration fields are (Liu & Liu, 2003)

$$\begin{aligned} \frac{d\mathbf{r}\_i}{dt} &= \mathbf{v}\_i, \\ \frac{d\mathbf{v}\_i^{\alpha}}{dt} &= \sum\_{j=1}^N m\_j \left( \frac{\sigma\_i^{\alpha\beta}}{\rho\_i^2} + \frac{\sigma\_j^{\alpha\beta}}{\rho\_j^2} \right) \cdot \nabla W\_{ij}^{\mathbf{h}}, \end{aligned} \tag{16}$$

where is the total stress tensor.

The internal energy evolution is given by the expression (Liu & Liu, 2003) :

$$\frac{dE\_{\dot{i}}}{dt} = \frac{1}{2} \sum\_{j=1}^{N} m\_{j} \left( \frac{p\_{\dot{i}}}{\rho\_{\dot{i}}^{2}} + \frac{p\_{\dot{j}}}{\rho\_{\dot{j}}^{2}} \right) \mathbf{v}\_{\dot{i}}^{\beta} - \mathbf{v}\_{\dot{j}}^{\beta} \Big| \frac{\partial \mathbf{W}\_{\dot{i}\dot{j}}}{\partial \mathbf{x}\_{\dot{i}}^{\beta}} + \frac{\mu\_{\dot{i}}}{2\rho\_{\dot{j}}} \mathbf{s}\_{\dot{i}}^{\alpha\beta} \mathbf{s}\_{\dot{i}}^{\alpha\beta},\tag{17}$$

In the above equation p is the pressure, is the dynamic viscosity and is the shear strain rate.

In the present work, our calculations are performed in three dimensions and we use the cubic B-spline kernel (Monaghan, 1985). We consider water drops, and the equation of state that we use in the hydrodynamical code was a general Mie-Gruneisen form of equation of state with different analytic forms for states of compression (/0-1)>0 and tension (/0- 1)<0 (Liu & Liu, 2003). This equation has several parameters, namely the density , the reference density 0, and the constants A1, A2, A3, C1 and C2 . The pressure P is

$$P = A\_1 \left(\frac{\rho}{\rho\_0} - 1\right) + A\_2 \left(\frac{\rho}{\rho\_0} - 1\right)^2 + A\_3 \left(\frac{\rho}{\rho\_0} - 1\right)^3 \quad \text{if} \quad \left(\frac{\rho}{\rho\_0} - 1\right) > 0 \tag{18}$$

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 93

Fig. 2. Sequence of times showing the evolution of the collision between two drops (permanent coalescence) with Vcol = 1.0 mm/ms and We = 4.5. The time scale is given in

milliseconds.

and

$$P = C\_1 \left(\frac{\rho}{\rho\_0} - 1\right) + C\_2 \left(\frac{\rho}{\rho\_0} - 1\right) \quad \text{if} \quad \left(\frac{\rho}{\rho\_0} - 1\right) < 0. \tag{19}$$

In all our calculations we use the following values for the constants: A1=2.20x106 kPa, A2=9.54x106 kPa, A3=1.46x107 kPa, C1=2.20x106 kPa, C2=0.00 kPa, and 0=1000.0 Kg/m3.

#### **3. Coalescence, fragmentation and flocculation of liquid drops in three dimensions**

In order to model the collision of liquid drops several calculations were carried out. We have varied the velocity of collision for modeling the permanent coalescence of droplets in the three dimensional space (3D) in a vacuum environment using the SPH method. In order to proceed we have defined drops with diameter of 30μm and 4700 SPH particles for each drop with a collision velocity of 1.0 mm/ms.

In figure 2 is illustrated a sequence of times showing the evolution of the collision between two drops (permanent coalescence) with Vcol = 1.0 mm/ms and We= 4.5. The evolution of time is shown in milliseconds. It can be seen in this figure that at t=0.0009 ms a flat circular section appears (Ivanov & Kralchevsky, 1997), which increases its diameter as dynamics progresses. The appearance of this flat circular section has been reported for the case of collision of drops in emulsion media (Ivanov & Kralchevsky, 1997), yet in this reference the hydrodynamic modeling of the collision of liquid drops is not considered. In the dynamics we observe that at t = 0.0053 ms a bridge structure between the two drops appears in the region of contact (this bridge is the structure that joins the drops through their flat circular interfaces placed in the center of the droplets coalescence), which disappears at a later time due to the penetration of particles of one drop into the other. After that, a process of coalescence occurs (see figure 2 at t=0.0069 ms) and a bigger drop is formed (see figure 2 at t=0.0077ms).

Figure 3 shows the velocity vector field inside the droplets and in the region of contact between them at t=0.004ms. Notice that inside the drops, the fluid tends to a velocity value lower than the initial velocity of 1.00 mm/ms, while in the area of contact between the drops we observe an increase in the fluid velocity to a value of 1.436 mm/ms. Once the coalescence process occurs, the velocity of the fluid inside the drops tends to zero, i.e. the largest drop size that is formed after some time tends to equilibrium. This can be seen in figure 4 where it is illustrated the time evolution of the kinetic and internal energy of the bigger drop. When two drops collide the possible results of the collision (coalescence, flocculation or fragmentation of drops) depends only on the kinetic energy and the Weber number (We) (Foote, 1974), which is given by

$$
\mathcal{W}e = \frac{\rho V\_r^2 d}{\sigma}.\tag{20}
$$

Here Vr is the difference between the velocities of the drops, d the diameter of the drop, and the surface tension.

92 Hydrodynamics – Optimizing Methods and Tools

0. (19)

 0 1

 *if*

In all our calculations we use the following values for the constants: A1=2.20x106 kPa, A2=9.54x106 kPa, A3=1.46x107 kPa, C1=2.20x106 kPa, C2=0.00 kPa, and 0=1000.0 Kg/m3.

In order to model the collision of liquid drops several calculations were carried out. We have varied the velocity of collision for modeling the permanent coalescence of droplets in the three dimensional space (3D) in a vacuum environment using the SPH method. In order to proceed we have defined drops with diameter of 30μm and 4700 SPH particles for each

In figure 2 is illustrated a sequence of times showing the evolution of the collision between two drops (permanent coalescence) with Vcol = 1.0 mm/ms and We= 4.5. The evolution of time is shown in milliseconds. It can be seen in this figure that at t=0.0009 ms a flat circular section appears (Ivanov & Kralchevsky, 1997), which increases its diameter as dynamics progresses. The appearance of this flat circular section has been reported for the case of collision of drops in emulsion media (Ivanov & Kralchevsky, 1997), yet in this reference the hydrodynamic modeling of the collision of liquid drops is not considered. In the dynamics we observe that at t = 0.0053 ms a bridge structure between the two drops appears in the region of contact (this bridge is the structure that joins the drops through their flat circular interfaces placed in the center of the droplets coalescence), which disappears at a later time due to the penetration of particles of one drop into the other. After that, a process of coalescence occurs (see figure 2 at t=0.0069 ms) and a bigger drop is formed (see figure 2 at

Figure 3 shows the velocity vector field inside the droplets and in the region of contact between them at t=0.004ms. Notice that inside the drops, the fluid tends to a velocity value lower than the initial velocity of 1.00 mm/ms, while in the area of contact between the drops we observe an increase in the fluid velocity to a value of 1.436 mm/ms. Once the coalescence process occurs, the velocity of the fluid inside the drops tends to zero, i.e. the largest drop size that is formed after some time tends to equilibrium. This can be seen in figure 4 where it is illustrated the time evolution of the kinetic and internal energy of the bigger drop. When two drops collide the possible results of the collision (coalescence, flocculation or fragmentation of drops) depends only on the kinetic energy and the Weber

> 2 . *V d*

(20)

Here Vr is the difference between the velocities of the drops, d the diameter of the drop, and

*We r*

 0 1

**3. Coalescence, fragmentation and flocculation of liquid drops in three**

*<sup>P</sup> <sup>C</sup>*<sup>1</sup>

drop with a collision velocity of 1.0 mm/ms.

number (We) (Foote, 1974), which is given by

 *C*<sup>2</sup>

 0 1

and

**dimensions**

t=0.0077ms).

the surface tension.

Fig. 2. Sequence of times showing the evolution of the collision between two drops (permanent coalescence) with Vcol = 1.0 mm/ms and We = 4.5. The time scale is given in milliseconds.

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 95

permanent coalescence occurs. The outcomes reported by Qian & Law (1997) are in good agreement with our results. In our SPH calculation, the relative velocity is not enough to produce fragmentation of the bigger drop and subsequently to produce small satellite droplets. In this calculation, the coalescence is permanent and the bigger drop that is formed reaches the equilibrium (see figure 4). On the other hand, the experiments of Qian & Law (1997) do not have a sufficient resolution to show in detail the deformation of the drops just before the formation of the bridge. However, the appearance of the flat circular section shown in figure 2 is in good agreement with the experimental and theoretical outcomes reported in the literature (Bibette et al., 1992; Ivanov & Dimitrov, 1988; Ivanov &

Kralchevsky, 1997; Kabalnov & Wennerström, 1996; Sharma & Ruckenstein, 1987).

Fig. 4. Evolution of the Kinetic and Internal energy for the collision between two equal-sized

On the other hand, it is observed that if we choose a Weber number for the collision greater than the range of values producing permanent coalescence, the phenomenon of fragmentation arises, i.e. the regime 2 reported by Qian & Law (1997) occurs giving rise to coalescence followed by separation into small satellite drops. The following calculations were performed for droplets with 30μm of diameter, 4700 SPH particles for each drop, and a collision velocity of 10.0 mm/ms (We=450) which is a characteristic velocity for the elements of a liquid spray (Choo & Kang, 2003). In the first stage of the calculation at t=2.0x10-4ms the collision of the two droplets is shown in figure 5. It can be seen the formation of a flat circular section between the drops (Ivanov & Kralchevsky, 1997). This circular section vanishes completely at t = 5.6 x10-4ms. At this time a portion of fluid appears to form a wave

drops with Vcol = 1.0 mm/ms and We = 4.5.

front propagating in the plane x = 0.

From the values of density, relative velocity, droplet diameter and surface tension we obtain the Weber number. The Surface tension is determined using the Laplace equation

$$\mathbb{E}\left[p(r=0) - p(r\to\infty)\right] = \frac{\sigma}{R} \,. \tag{21}$$

The first term p (r = 0) on the left side of the equation (21) is determined at the drop center and the second term p (r → ∞) is taken as the vanishing pressure far away from the drop. The calculations were made in a vacuum environment and only head-on collisions were considered. The value of the pressure at the drop center is 1.78kPa and the Weber number for the coalescence collision is We=4.5. Values of the Weber number in the range 1We19 have been chosen, which corresponds to the range reported by Ashgriz & Poo (1990) for experimental head-on collisions and coalescence of water drops.

Fig. 3. Velocity vector field for the collision between two drops at t=0.004 ms (permanent coalescence process) with Vcol = 1.0 mm/ms and We = 4.5. The time scale is given in milliseconds.

The coalescence processes occurring in droplets collision at low Weber numbers illustrated in the figures 1, 2 and 3 are governed by the competition between stretching and drop drainage. This drainage occurs when the liquid flows from the high-pressure drop region toward the point of contact to form a liquid bridge.

For larger values of We, the initially merged droplets would subsequently split apart, with the simultaneous production of smaller satellite droplets. It can be observed in our simulations (see figure 2) that the circular flat section disappears. Also, due to the effect of the surface tension a structure having the form of bridge is observed at t=0.0053ms and the 94 Hydrodynamics – Optimizing Methods and Tools

From the values of density, relative velocity, droplet diameter and surface tension we obtain

*<sup>R</sup>* . (21)

*p r* <sup>0</sup> *p r*

Fig. 3. Velocity vector field for the collision between two drops at t=0.004 ms (permanent coalescence process) with Vcol = 1.0 mm/ms and We = 4.5. The time scale is given in

The coalescence processes occurring in droplets collision at low Weber numbers illustrated in the figures 1, 2 and 3 are governed by the competition between stretching and drop drainage. This drainage occurs when the liquid flows from the high-pressure drop region

For larger values of We, the initially merged droplets would subsequently split apart, with the simultaneous production of smaller satellite droplets. It can be observed in our simulations (see figure 2) that the circular flat section disappears. Also, due to the effect of the surface tension a structure having the form of bridge is observed at t=0.0053ms and the

milliseconds.

toward the point of contact to form a liquid bridge.

The first term p (r = 0) on the left side of the equation (21) is determined at the drop center and the second term p (r → ∞) is taken as the vanishing pressure far away from the drop. The calculations were made in a vacuum environment and only head-on collisions were considered. The value of the pressure at the drop center is 1.78kPa and the Weber number for the coalescence collision is We=4.5. Values of the Weber number in the range 1We19 have been chosen, which corresponds to the range reported by Ashgriz & Poo (1990) for

the Weber number. The Surface tension is determined using the Laplace equation

experimental head-on collisions and coalescence of water drops.

permanent coalescence occurs. The outcomes reported by Qian & Law (1997) are in good agreement with our results. In our SPH calculation, the relative velocity is not enough to produce fragmentation of the bigger drop and subsequently to produce small satellite droplets. In this calculation, the coalescence is permanent and the bigger drop that is formed reaches the equilibrium (see figure 4). On the other hand, the experiments of Qian & Law (1997) do not have a sufficient resolution to show in detail the deformation of the drops just before the formation of the bridge. However, the appearance of the flat circular section shown in figure 2 is in good agreement with the experimental and theoretical outcomes reported in the literature (Bibette et al., 1992; Ivanov & Dimitrov, 1988; Ivanov & Kralchevsky, 1997; Kabalnov & Wennerström, 1996; Sharma & Ruckenstein, 1987).

Fig. 4. Evolution of the Kinetic and Internal energy for the collision between two equal-sized drops with Vcol = 1.0 mm/ms and We = 4.5.

On the other hand, it is observed that if we choose a Weber number for the collision greater than the range of values producing permanent coalescence, the phenomenon of fragmentation arises, i.e. the regime 2 reported by Qian & Law (1997) occurs giving rise to coalescence followed by separation into small satellite drops. The following calculations were performed for droplets with 30μm of diameter, 4700 SPH particles for each drop, and a collision velocity of 10.0 mm/ms (We=450) which is a characteristic velocity for the elements of a liquid spray (Choo & Kang, 2003). In the first stage of the calculation at t=2.0x10-4ms the collision of the two droplets is shown in figure 5. It can be seen the formation of a flat circular section between the drops (Ivanov & Kralchevsky, 1997). This circular section vanishes completely at t = 5.6 x10-4ms. At this time a portion of fluid appears to form a wave front propagating in the plane x = 0.

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 97

Fig. 5. Sequence of times showing the evolution of the collision between two drops with Vcol = 10.0 mm/ms and We = 450. This figure illustrates the formation of small satellite

droplets. The time scale is given in milliseconds.

This wave front begins to form little satellite drops and increases its amplitude until t=1.8x10-3ms. There is no substantial growth of these satellite drops and the structure tends to a flatten form as the dynamic runs. Figure 6 shows the velocity vector field (seen from the plane y = 0) after the fragmentation of the drops has taken place. As shown in figure 6, the fluid velocity at the center of the structure is 8.7 mm/ms, which is less than the initial rate of collision, while the fluid that is spread to the edges is accelerated reaching a speed of 15.0 mm/ms. A longer stretched ligament is produced and the amount of satellite droplets increases with the evolution of dynamics. Figure 5 illustrates that a portion of the fluid begins to separate, stretching away from the bigger drop, and a non-uniform pressure field is created inside the ligament. This is related to the value in the velocity vector field differences, and due to this pressure differential a flow parallel to the plane x=0 is produced. The fragmentation phenomena and the subsequent formation of satellites drops may be analyzed following the conjectures made by Qian & Law (1997).

Once the ligament begins to form (see figure 5), a flow is generated directed in the opposite direction to the vector field shown in figure 6. This motion transform this portion of the ligament in a bulbous due to the accumulation of mass in this volume, and this change in geometry implies the appearance of a local minimum in the pressure field which is located between the bulbous and its neighboring region.

As a result of this pressure difference a local flow is generated through the point of minimum pressure that opposes to the flow that is coming from the bulbous. This fluid motion causes a local reduction of mass and therefore the ligament between the bulbous and the neighboring region starts to decrease its radius (at the point of minimal pressure). Because of this local decrease of the ligament radius the pressure rises, which creates a flow with the same direction of the flow that comes from the end of the bulbous and other flow in the opposite direction from the point of local reduction of mass. Given these opposing flows emerging from this point, the radius of the ligament decreases even more. Then the system tends to relax this unstable situation reducing the radius of this region to zero, giving rise to a division of the fluid and so producing a satellite drop (see figure 5). Subsequently, this process is repeated in the other regions of the ligament, producing more satellites drops. These shattering collisions occur only at high velocities making the surface tension forces of secondary importance (the phenomenon is inertial dominated).

When the Weber number for the collision is decreased below the range corresponding to the permanent coalescence regime, then flocculation occurs. These calculations were performed for droplets with 30μm of diameter, 4700 SPH particles for each drop, and a collision velocity of 0.2mm/ms (We=0.18). At the beginning of the calculation one observes at t=0.29ms (see figure 7) that a flat circular section appears between the two droplets (Ivanov & Kralchevsky, 1997), which has already increased in diameter at t=1.0ms. Then, there is a stretching of the surface of the drop as can be seen at t=1.77ms. This stretch is deforming the drops until t=3.76 ms, and after that the drop shape remains constant. The chosen collision velocity cannot produce coalescence between the droplets. In fact no penetration was observed through the plane x=0. In this case, only the drops stay together, interacting through their surfaces, giving rise to flocs (Ivanov & Kralchevsky, 1997).

It has been reported that these flocs are formed in emulsions when the interfacial film between drops is very stable or the drops approach each other with a very small kinetic energy (Ivanov & Dimitrov, 1988; and Ivanov & Kralchevsky, 1997). In this case, the wave front that appears in the plane x=0 of figure 5 was not observed. Figure 8 shows the velocity vector field at t = 1.78ms. It can be seen in this figure that the value of the fluid velocity 96 Hydrodynamics – Optimizing Methods and Tools

This wave front begins to form little satellite drops and increases its amplitude until t=1.8x10-3ms. There is no substantial growth of these satellite drops and the structure tends to a flatten form as the dynamic runs. Figure 6 shows the velocity vector field (seen from the plane y = 0) after the fragmentation of the drops has taken place. As shown in figure 6, the fluid velocity at the center of the structure is 8.7 mm/ms, which is less than the initial rate of collision, while the fluid that is spread to the edges is accelerated reaching a speed of 15.0 mm/ms. A longer stretched ligament is produced and the amount of satellite droplets increases with the evolution of dynamics. Figure 5 illustrates that a portion of the fluid begins to separate, stretching away from the bigger drop, and a non-uniform pressure field is created inside the ligament. This is related to the value in the velocity vector field differences, and due to this pressure differential a flow parallel to the plane x=0 is produced. The fragmentation phenomena and the subsequent formation of satellites drops may be

Once the ligament begins to form (see figure 5), a flow is generated directed in the opposite direction to the vector field shown in figure 6. This motion transform this portion of the ligament in a bulbous due to the accumulation of mass in this volume, and this change in geometry implies the appearance of a local minimum in the pressure field which is located

As a result of this pressure difference a local flow is generated through the point of minimum pressure that opposes to the flow that is coming from the bulbous. This fluid motion causes a local reduction of mass and therefore the ligament between the bulbous and the neighboring region starts to decrease its radius (at the point of minimal pressure). Because of this local decrease of the ligament radius the pressure rises, which creates a flow with the same direction of the flow that comes from the end of the bulbous and other flow in the opposite direction from the point of local reduction of mass. Given these opposing flows emerging from this point, the radius of the ligament decreases even more. Then the system tends to relax this unstable situation reducing the radius of this region to zero, giving rise to a division of the fluid and so producing a satellite drop (see figure 5). Subsequently, this process is repeated in the other regions of the ligament, producing more satellites drops. These shattering collisions occur only at high velocities making the surface tension forces of

When the Weber number for the collision is decreased below the range corresponding to the permanent coalescence regime, then flocculation occurs. These calculations were performed for droplets with 30μm of diameter, 4700 SPH particles for each drop, and a collision velocity of 0.2mm/ms (We=0.18). At the beginning of the calculation one observes at t=0.29ms (see figure 7) that a flat circular section appears between the two droplets (Ivanov & Kralchevsky, 1997), which has already increased in diameter at t=1.0ms. Then, there is a stretching of the surface of the drop as can be seen at t=1.77ms. This stretch is deforming the drops until t=3.76 ms, and after that the drop shape remains constant. The chosen collision velocity cannot produce coalescence between the droplets. In fact no penetration was observed through the plane x=0. In this case, only the drops stay together, interacting

It has been reported that these flocs are formed in emulsions when the interfacial film between drops is very stable or the drops approach each other with a very small kinetic energy (Ivanov & Dimitrov, 1988; and Ivanov & Kralchevsky, 1997). In this case, the wave front that appears in the plane x=0 of figure 5 was not observed. Figure 8 shows the velocity vector field at t = 1.78ms. It can be seen in this figure that the value of the fluid velocity

analyzed following the conjectures made by Qian & Law (1997).

secondary importance (the phenomenon is inertial dominated).

through their surfaces, giving rise to flocs (Ivanov & Kralchevsky, 1997).

between the bulbous and its neighboring region.

Fig. 5. Sequence of times showing the evolution of the collision between two drops with Vcol = 10.0 mm/ms and We = 450. This figure illustrates the formation of small satellite droplets. The time scale is given in milliseconds.

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 99

Fig. 7. Sequence of times showing the evolution of the collision between two drops

(flocculation) with Vcol = 0.2 mm/ms and We = 0.18. The time scale is given in milliseconds.

decreases from the border to the center of the drop. In this case, the velocity has decreased below its initial value, which is 0.2mm/ms in all zones of the fluid. As is shown in figure 9, after an elapsed time of 3.76ms, the fluid velocities decrease even more, reaching a value of 2.24x10-2 mm/ms in the zone of interaction between the two drops and 4.48x10-2 mm/ms near the border.

Fig. 6. Velocity vector field showing the fragmentation of two colliding drops at t=1.2x10-3ms (see from the plane z-x) with Vcol=10.0 mm/ms and We = 450. The time scale is given in milliseconds.

98 Hydrodynamics – Optimizing Methods and Tools

decreases from the border to the center of the drop. In this case, the velocity has decreased below its initial value, which is 0.2mm/ms in all zones of the fluid. As is shown in figure 9, after an elapsed time of 3.76ms, the fluid velocities decrease even more, reaching a value of 2.24x10-2 mm/ms in the zone of interaction between the two drops and 4.48x10-2 mm/ms

Fig. 6. Velocity vector field showing the fragmentation of two colliding drops at t=1.2x10-3ms (see from the plane z-x) with Vcol=10.0 mm/ms and We = 450. The time

scale is given in milliseconds.

near the border.

Fig. 7. Sequence of times showing the evolution of the collision between two drops (flocculation) with Vcol = 0.2 mm/ms and We = 0.18. The time scale is given in milliseconds.

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 101

Fig. 9. Velocity vector field showing the flocculation of two liquid drops at t=3.76 ms with

Vcol = 0.2mm/ms and We = 0.18. The time scale is given in milliseconds.

Fig. 8. Velocity vector field showing the flocculation of two liquid drops at t=1.78 ms with Vcol = 0.2mm/ms and We = 0.18. The time scale is given in milliseconds.

100 Hydrodynamics – Optimizing Methods and Tools

Fig. 8. Velocity vector field showing the flocculation of two liquid drops at t=1.78 ms with

Vcol = 0.2mm/ms and We = 0.18. The time scale is given in milliseconds.

Fig. 9. Velocity vector field showing the flocculation of two liquid drops at t=3.76 ms with Vcol = 0.2mm/ms and We = 0.18. The time scale is given in milliseconds.

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 103

Bozzano, G. & Dente, M. (2010). Mechanism of drop coalescence at the interface of two

Bradley, S.G. & Stow, C.D. (1978). Collision between liquid drops. *Philos. Trans. R. Soc.*

Brenn, G. & Frohn, A. (1989). Collision and merging of two equal droplets of propanol. *Exp.*

Brenn, G. & Kolobaric, V. (2006). Satellite droplet formation by unstable binary drop

Chen, J.D. (1985). A model of coalescence between two equal-sized spherical drops or

Choo, Y.J. & Kang, B.S. (2003). A study on the velocity characteristics of the liquid elements produced by two impinging jets. *Exp. in Fluids*, 34, 2003, pp. 655-661; Colagrossi, A. & Landrini, M. (2003). Numerical simulation of interfacial flows by smoothed

Cristini, V., Bawzdziewicz, J. & Loewenberg, M. (2001). An Adaptive Mesh Algorithm for

Cui, J., He, G.W. & Qi, D.W. (2006). A constrained particle dynamics for continuum-particle

Decent, S.P., Sharpe, G., Shaw, A.J. & Suckling, P.M. (2006). The formation of a liquid bridge during the coalescence of drops. Int. J. *Multi-phase Flow*, 32, 2006, pp. 717-738; Desbrun, M. & Gascuel, M.P. (1996). Smoothed Particles: a new paradigm for animating

Duchemin, L., Eggers, J. & Josserand, C. (2003). Inviscid coalescence of drops. *J. Fluid Mech*.,

Eggers, J., Lister, J.R. & Stone, H.A. (1999). Coalescence of liquid drops. *J. Fluid Mech.,* 401,

Foote, G.B. (1974). The water drop rebound problem: Dynamics of collision. *J. Atmos. Sci*., 32,

Gingold, R.A. & Monaghan, J.J. (1977). Smoothed particle hydrodynamics: theory and application to non-spherical stars. *Roy. Astronom. Soc*., 181, 1977, pp. 375-389; Gokhale, S.J., Dasgupta, S., Plawsky, J.L. & Wayner, P.C. (2004). Reflectivity-based

Gotaas, C., Havelka, P., Jakobsen, H., Svendsen, H., Hase, M., Roth, N. & Weigand, B. (2007,

Gotaas, C., Havelka, P., Jakobsen, H. & Svendsen, H. (2007, b). Evaluation of the impact

coalesced drop. *Phys. Rev. E.*, 70, 2004, pp. 1-12;

numerical comparison. *Phys. Fluids*, 19, 2007, pp. 1-17;

evaluation of the coalescence of two condensing drops and shape evolution of the

a). Effect of viscosity on droplet-droplet collision outcome: Experimental study and

parameter in droplet-droplet collision experiments by the aliasing method. *Phys.*

Evolving Surfaces: Simulations of Drop Breakup and Coalescence. *J. Comput. Phys*.,

hybrid method in micro- and nano fluidics. Acta Mechanica Sinica, 22, 2006, pp.

highly deformable bodies. *Proceedings of Eurographics Workshop on animation and*

particle hydrodynamics. *J. Comput. Phys*., 191, 2003, pp. 448-475;

Engineering – ESCAPE20, 28, 2010, pp. 55-60;

*London Ser. A,* 287, 1978, pp. 635-642;

collisions. *Phys. Fluids*, 18, 2006, pp. 1-18;

bubbles. *J. Colloid Interface Sci.*, 107, 1985, pp. 209-220;

*Fluids*, 7, 1989, pp. 441-446;

168, 2001, pp. 445-463;

503-508;

*simulation*, 1996.

1999, pp. 293-310;

1974, pp. 390-401;

*Fluids*, 19, 2007, pp. 1-14;

487, 2003, pp. 167-180;

immiscible liquids. 20th European Symposium on Computer Aided Process

#### **4. Conclusions**

An adequate methodology using the SPH method in three-dimensional space was presented for the calculation of hydrodynamic collisions of liquid drops in a vacuum environment. Several features of binary collisions in three dimensions have been observed in our calculations. As a result of the collision between droplets the formation of a circular flat film was obtained for a range of values of the Weber number. We have found three possible outcomes for the collision: permanent coalescence, fragmentation and flocculation of drops. We have ascribed a range of Weber number values for the three possible outcomes of the collisions.

Inside the drops, the fluid tends to a velocity value lower than the initial velocity of collision, while in the area of contact between the drops, we observed an increase in the fluid velocity. When the Weber number of the drops is larger than the range of values corresponding to coalescence, a fragmentation phenomenon is observed. Otherwise, when the Weber number is below this range, we obtained drop flocculation. In the case of drop fragmentation, we have observed the formation of satellite droplets emerging from the contact zone between the two drops. The mechanism of formation of these satellite droplets has been discussed.

We can see in the velocity vector field that after the fragmentation of the drops has taken place, the fluid velocity at the center of the structure is less than the initial rate of collision, while the fluid that is spread to the edges is accelerated producing a long stretched ligament and several satellite droplets. In the SPH dynamics that the drops follow when flocculation occurs, we have obtained that the system tends to a state of equilibrium in which the kinetic energy remains constant and the drops interact mainly through their deformed surfaces. On the velocity vector fields can be seen that in this case the value of the fluid velocity decreases from the border to the center of the drops. In all zones of the fluid the velocity has decreased below its initial value. After an elapsed time the fluid velocities decrease even more in the zone of interaction between the two drops near the border.

#### **5. References**

102 Hydrodynamics – Optimizing Methods and Tools

An adequate methodology using the SPH method in three-dimensional space was presented for the calculation of hydrodynamic collisions of liquid drops in a vacuum environment. Several features of binary collisions in three dimensions have been observed in our calculations. As a result of the collision between droplets the formation of a circular flat film was obtained for a range of values of the Weber number. We have found three possible outcomes for the collision: permanent coalescence, fragmentation and flocculation of drops. We have ascribed a range of Weber number values for the three possible outcomes of the

Inside the drops, the fluid tends to a velocity value lower than the initial velocity of collision, while in the area of contact between the drops, we observed an increase in the fluid velocity. When the Weber number of the drops is larger than the range of values corresponding to coalescence, a fragmentation phenomenon is observed. Otherwise, when the Weber number is below this range, we obtained drop flocculation. In the case of drop fragmentation, we have observed the formation of satellite droplets emerging from the contact zone between the two drops. The mechanism of formation of these satellite droplets

We can see in the velocity vector field that after the fragmentation of the drops has taken place, the fluid velocity at the center of the structure is less than the initial rate of collision, while the fluid that is spread to the edges is accelerated producing a long stretched ligament and several satellite droplets. In the SPH dynamics that the drops follow when flocculation occurs, we have obtained that the system tends to a state of equilibrium in which the kinetic energy remains constant and the drops interact mainly through their deformed surfaces. On the velocity vector fields can be seen that in this case the value of the fluid velocity decreases from the border to the center of the drops. In all zones of the fluid the velocity has decreased below its initial value. After an elapsed time the fluid velocities decrease even more in the

Ashgriz, N. & Givi, P. (1987). Binary collision dynamics of fuel droplets. *Int. J. Heat Fluid*

Ashgriz, N. & Givi, P. (1989). Coalescence efficiencies of fuel droplets in binary collisions.

Ashgriz, N. & Poo, J.Y. (1990). Coalescence and separation of binary collisions of liquid

Azizi, F. & Al Taweel, A.M. (2010). Algorithm for the accurate numerical solution of PBE for

Benz, W. & Asphaug. (1994). Impact simulations with fracture: I. Method and tests. *Icarus*,

Benz, W. & Asphaug. (1995). Simulations of brittle solids using smoothed particle

Bibette, J., Morse, D.C., Witten, T.A. & Weitz, D.A. (1992). Stability criteria for emulsions.

hydrodynamics. *Comput. Phys. Commun.,* 87, 1995, pp. 253-265;

drop breakup and coalescence under high shear rates. *Chem. Eng. Sci*., 65, 2010, pp.

zone of interaction between the two drops near the border.

*Int. Commun. Heat Mass Transfer,* 16, 1989, pp. 11-17;

drops. *J. Fluid Mech.,* 221, 1990, pp. 183-204;

*Phys. Rev. Lett*., 69, 1992, pp. 2439-2443;

*Flow,* 8, 1987, pp. 205-208;

**4. Conclusions**

collisions.

has been discussed.

**5. References**

6112-6127;

1233, 1994, pp. 98-116.

3D Coalescence Collision of Liquid Drops Using Smoothed Particle Hydrodynamics 105

Nobari, M.R., Jan, Y.J. & Tryggvason, G. (1996). Head-on collision of drops-A numerical

O'Connell, S.T. & Thompson, P.A. (1995). Molecular dynamics-continuum hybrid

Pan, Y. & Suga, K. (2005). Numerical simulation of binary liquid droplet collision. *Phys.*

Park, J.Y. & Blair, L.M. (1975). The effect of coalescence on drop size distribution in an agitated liquid–liquid dispersion. *Chem. Eng. Sci.,* 30, 1975, pp. 1057-1064; Podgorska, W. (2007). Influence of Dispersed Phase Viscosity on Drop Coalescence in

Qian, J. & Law, C.K. (1997). Regimes of coalescence and separation in droplet collision. *J.*

Roisman, L. (2004). Dynamics of inertia dominated binary drops collisions. *Phys. Fluids,* 16,

Roisman, L., Berberovic, E. & Tropea, C. (2009). Inertia dominated drops collisions. I. On the

Rourke, P.J. & Bracco, F.V. (1980). Modelling of drop interactions in thick sprays and a

Shah, P.S., Fan, L.T., Kao, I.C. & Erickson, L.E. (1972). Modeling of growth processes with

Sharma, A. & Ruckenstein, E. (1987). Stability, critical thickness, and the time of rupture of

Sun, Z., Xi, G. & Chen, X. (2009). Mechanism study of deformation and mass transfer for binary droplet collisions with particle method. *Phys. Fluids*, 21, 2009, pp. 1-13; Tartakovsky, A.M. & Meakin, P.A. (2005). Smoothed particle hydrodynamics model for

Thoroddsen, S.T., Qian, B., Etoh, T.G. & Takehara, K. (2007). The initial coalescence of

Wang, W., Gong, J., Ngan, K.H. & Angeli, P. (2009). Effect of glycerol on the binary

Wu, M., Cubaud, T. & Ho, C. (2004). Scaling law in liquid drop coalescence driven by

Xing, X.Q., Butler, D.L., Ng, S.H., Wang, Z., Danyluk, S. & Yang, C. (2007). Simulation of

Yiantsios, S.G. & Davis, R.H. (1991). Close approach and deformation of two viscous drops

comparison with experiments, in: Stratified Charged Auto Engineering Conference.

two liquid phases: a review of drop phenomena, mixing, and growth. *Adv. Appl.*

miscible flow in three-dimensional fractures and the two-dimensional Rayleigh–

coalescence of water drops in stagnant oil phase. *Chem. Eng. Res. Design*, 87, 2009,

droplet formation and coalescence using lattice Boltzmann-based single-phase

due to gravity and van der Waals forces. *J. Colloid Interface Sci*., 144, 1991, pp. 412-

Turbulent Flow. *Chem. Eng. Res. Design,* 85, 2007, pp. 721-729;

universal flow in the lamella. *Phys. Fluids*, 21, 2009, pp. 1-10;

*Institute of Mechanical Engineering Publications*, 1980, pp. 101-116;

thinning foam and emulsion films. *Langmuir*, 3, 1987, pp. 760-768;

Taylor instability. *J. Comput. Phys*., 207, 2005, pp. 610-624.

miscible drops. *Phys. Fluids*, 19, 2007, 1-21;

surface tension. *Phys. Fluids*, 16, 2004, pp. 51-54;

model. *J. Coll. Interf. Sci*., 311, 2007, pp. 609-618;

computations: a tool for studying complex fluid flows. *Phys. Rev. E,* 52, 1995, pp.

investigation. Phys. Fluids, 8, 1996, pp. 29-42;

5792-5795.

*Fluids,* 17, 2005, pp. 1-14;

2004, pp. 3438-3449;

pp. 1640-1648;

433;

*Fluid Mech*., 331, 1997, pp. 59-80.

*Microbiol.,* 15, 1972, pp. 367-414;

104 Hydrodynamics – Optimizing Methods and Tools

Hoover, W.G. (1998). Isomorphism linking smooth particles and embedded atoms. *Physica*

Ivanov, I.B. & Dimitrov, D.S. (1988). Thin Liquid Films Fundamentals and Applications.

Ivanov, I.B. & Kralchevsky, P.A. (1997). Stability of emulsions under equilibrium and

Jia, X., McLaughlin, J.B. & Kontomaris, K. (2006). Lattice Boltzmann simulations of drop

Jiang, Y.J., Umemura, A. & Law, C.K. (1992). An experimental investigation on the collision behaviour of hydrocarbon droplets. *J. Fluid Mech*., 234, 1992, pp. 171-190; Kabalnov, A. & Wennerström, H. (1996). Macroemulsion stability: The oriented wedge

Koumoutsakos, P. (2005). Multiscale flow simulation using particle. *Ann. Rev. Fluid Mech.*,

Leal, L.G. (2004). Flow induced coalescence of drops in a viscous fluid. *Phys. Fluids*, 16, 2004,

Li, D. (1994). Coalescence between two small bubbles or drops. *J. Colloid Interface Sci*., 163,

Li, J., Liao, D. & Yip, S. (1998). Coupling continuum to molecular-dynamics simulation:

Liu, G.R. & Liu, M.B. (2003). Smoothed Particle Hydrodynamics – A Meshfree Particle

Lucy, L.B. (1977). A numerical approach to the testing of the fission hypothesis. *Astron. J.,* 82,

Mashayek, F., Ashgriz, N., Minkowycz, W.J. & Shotorban, B. (2003). Coalescence collision of

Menchaca-Rocha, A., Huidobro, F., Martinez-Davalos, A., Michaelian, K., Perez, A.,

Menchaca-Rocha, A., Martinez-Davalos, A., Nuñez, R. (2001). Coalescence of liquid drops

Mohamed-Kassim, Z. & Longmire, E.K. (2004). Drop coalescence through a liquid/liquid

Monaghan, J.J. (1994). Simulating Free Surface Flows with SPH. *J. Comput. Phys.,* 110, 1994,

Monaghan, J.J. (1992). Smoothed particle hydrodynamics. *Annu. Rev. Astron. Astrophys,* 1992,

Monaghan, J.J. (1985). Extrapolating B splines for interpolation. *J. Comput. Phys.,* 60, 1985,

Narsimhan, G. (2004). Model for drop coalescence in a locally isotropic turbulent flow field.

Nie, X.B., Chen, S.Y., WN, E. & Robbins, M.O. (2004). A continuum and molecular dynamics hybrid method for micro- and nano-fluid flow. *J. Fluid Mech*., 500, 2004, pp. 55-64.

Rodriguez, V. & Carjan, N. (1997). Coalescence and fragmentation of colliding

liquid drops. *Int. J. Heat Mass Trans.,* 46, 2003, pp. 77-89;

mercury drops. *J. Fluid Mech.,* 346, 1997, pp. 291-318;

by surface tension. *Phys. Rev. E.,* 63, 2001, pp. 1-5;

interface. *Phys. Fluids*, 2004, pp. 1-47;

*J. Coll. Interf. Sci*., 272, 2004, pp. 197-209;

reflecting particle method and the field estimator. *Phys. Rev. E,* 57, 1998, pp. 7259-

dynamic conditions. *Coll. Surf. A*, 128, 1997, pp. 155-175;

theory revisited. *Langmuir*, 12, 1996, pp. 276-292;

coalescence and chemical mixing. *Physica A*, 362, 2006, pp. 62-67;

*A,* 260, 1998, pp. 244-255;

37, 2005, pp. 457-487.

pp. 1833-1851;

7267;

1994, pp. 108-119;

1977, pp. 1013-1024;

pp. 399-406;

pp. 543-574;

pp. 253-262;

Method. World Scientific, 2003;

Dekker, 1988.

**Part 2**

**Models and Codes in Fluid Dynamics**

## **Part 2**

## **Models and Codes in Fluid Dynamics**

106 Hydrodynamics – Optimizing Methods and Tools

Yoon, Y., Baldessari, F., Ceniceros, H.D. & Leal, L.G. (2007). Coalescence of two equal-sized deformable drops in an axisymmetric flow. *Phys*. *Fluids*, 19, 2007, pp. 1-24; Zhang, F.H., Li, E.Q. & Thoroddsen, S.T. (2009). Satellite formation during coalescence of

unequal size drops. *Phys. Rev. Lett.*, 102, 2009, pp. 1-4;

**0**

**6**

*Chile*

**Eulerian-Lagrangian Formulation for**

Carlos Cartes and Orazio Descalzi

*Complex Systems Group, Universidad de los Andes*

**Compressible Navier-Stokes Equations**

The Eulerian-Lagrangian formulation of the (inviscid) Euler dynamics in terms of advected Weber-Clebsch potentials (Lamb, 1932), was extended by Constantin to cover the viscous Navier-Stokes dynamics (Constantin, 2001). Numerical studies (Ohkitani & Constantin, 2003), of this formulation of the Navier-Stokes equations concluded that the diffusive Lagrangian map becomes non-invertible under time evolution and requires resetting for its calculation. They proposed the observed sharp increase of the frequency of resettings as a new diagnostic

In previous work we were able (Cartes et al., 2007; 2009) to complement these results, using an approach that is based on a generalised set of equations of motion for the Weber-Clebsch potentials, that turned out to depend on a parameter *τ*, which has the unit of time for the Navier-Stokes case. Also to extend our formulation to magnetohydrodynamics, and thereby

In this work we present a generalisation of the Weber-Clebsch variables in order to describe the compressible Navier-Stokes dynamics. Our main result is a good agreement between the dynamics for the velocity and density fields that come from the dynamics of Weber-Clebsch variables and direct numerical simulations of the compressible Navier-Stokes equations. We first present the inviscid Eulerian-Lagrangian theory, then Constantin's extension to viscous fluids and derive our equations of motion for the Weber-Clebsch potentials that describe the compressible Navier-Stokes dynamics. Then, performing direct numerical simulations of the Taylor-Green vortex, we check that our formulation reproduces the

Let us consider the incompressible Euler equations with constant density, fixed to one, for the

∇ · **u** = 0 ,

*∂t***u** + **u** · ∇**u** = −∇*p* (1)

**1. Introduction**

of vortex reconnection.

compressible dynamics.

velocity field **u**

**2. Eulerian-Lagrangian theory**

**2.1 Euler equations and Clebsch variables**

obtain a new diagnostic for magnetic reconnection.

## **Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations**

Carlos Cartes and Orazio Descalzi *Complex Systems Group, Universidad de los Andes Chile*

#### **1. Introduction**

The Eulerian-Lagrangian formulation of the (inviscid) Euler dynamics in terms of advected Weber-Clebsch potentials (Lamb, 1932), was extended by Constantin to cover the viscous Navier-Stokes dynamics (Constantin, 2001). Numerical studies (Ohkitani & Constantin, 2003), of this formulation of the Navier-Stokes equations concluded that the diffusive Lagrangian map becomes non-invertible under time evolution and requires resetting for its calculation. They proposed the observed sharp increase of the frequency of resettings as a new diagnostic of vortex reconnection.

In previous work we were able (Cartes et al., 2007; 2009) to complement these results, using an approach that is based on a generalised set of equations of motion for the Weber-Clebsch potentials, that turned out to depend on a parameter *τ*, which has the unit of time for the Navier-Stokes case. Also to extend our formulation to magnetohydrodynamics, and thereby obtain a new diagnostic for magnetic reconnection.

In this work we present a generalisation of the Weber-Clebsch variables in order to describe the compressible Navier-Stokes dynamics. Our main result is a good agreement between the dynamics for the velocity and density fields that come from the dynamics of Weber-Clebsch variables and direct numerical simulations of the compressible Navier-Stokes equations.

We first present the inviscid Eulerian-Lagrangian theory, then Constantin's extension to viscous fluids and derive our equations of motion for the Weber-Clebsch potentials that describe the compressible Navier-Stokes dynamics. Then, performing direct numerical simulations of the Taylor-Green vortex, we check that our formulation reproduces the compressible dynamics.

#### **2. Eulerian-Lagrangian theory**

#### **2.1 Euler equations and Clebsch variables**

Let us consider the incompressible Euler equations with constant density, fixed to one, for the velocity field **u**

$$\begin{aligned} \partial\_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} &= -\nabla p \\ \nabla \cdot \mathbf{u} &= 0 \end{aligned} \tag{1}$$

[∇, *Dt*] ≡ (∇**u**) · ∇, (10)

. (11)

(*a*, *t*) its position at

(*a*, *t*) then the Lagrangian

we can obtain the evolution equation for **u**

Newton equations for the fluid element are

*<sup>X</sup>* (*a*, *t*) are given by

and *p* (*X* (*a*, *t*)) is the pressure field in Eulerian coordinates. Therefore the movement equations for the fluid elements are

*Fi*

*∂*2*X<sup>i</sup>*

systems. We also note that this transformation is always invertible.

**2.2 Weber transformation**

velocity of a fluid element is

time *t* and note *A<sup>i</sup>*

and its acceleration

where the forces *F<sup>i</sup>*

coordinates, verifies

*Dt***<sup>u</sup>** <sup>=</sup> −∇

Let us note *a<sup>i</sup>* as the initial coordinate (at *t* = 0) of a fluid element and *X<sup>i</sup>*

(*x*, *<sup>t</sup>*) the inverse application: *<sup>a</sup><sup>i</sup>* <sup>≡</sup> *<sup>A</sup><sup>i</sup>*

*u*˜ *i*

*∂u*˜*<sup>i</sup>*

*∂*2*X<sup>i</sup>*

(*a*, *<sup>t</sup>*) = *<sup>∂</sup>X<sup>i</sup>*

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) = *<sup>∂</sup>*2*X<sup>i</sup>*

*<sup>∂</sup>t*<sup>2</sup> (*a*, *<sup>t</sup>*) <sup>=</sup> *<sup>F</sup><sup>i</sup>*

*<sup>X</sup>* (*a*, *<sup>t</sup>*) <sup>=</sup> <sup>−</sup> *<sup>∂</sup><sup>p</sup>*

*<sup>∂</sup>t*<sup>2</sup> (*a*, *<sup>t</sup>*) <sup>=</sup> <sup>−</sup> *<sup>∂</sup><sup>p</sup>*

det

From Eq. (16) we perform a coordinate transformation for the derivatives of *p* using

*∂ <sup>∂</sup>x<sup>i</sup>* <sup>=</sup> *<sup>∂</sup>A<sup>j</sup> ∂x<sup>i</sup>*

*∂ <sup>∂</sup>a<sup>i</sup>* <sup>=</sup> *<sup>∂</sup>X<sup>j</sup> ∂a<sup>i</sup>*

For an incompressible fluid, the transformation matrix, between Lagrangian and Eulerian

*∂X<sup>i</sup> ∂a<sup>j</sup>*

this value is fixed from the relation between the volume elements in the two coordinate

*∂*

*∂ ∂x<sup>j</sup>*

At time *t* Eulerian coordinates are by definition the variables *x<sup>i</sup>* = *X<sup>i</sup>*

*Dtφ* + 1 2 **u**2

Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations 111

(*X<sup>i</sup>*

(*a*, *t*), *t*).

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) (12)

*<sup>∂</sup>t*<sup>2</sup> (*a*, *<sup>t</sup>*) . (13)

*<sup>X</sup>* (*a*, *t*) , (14)

*<sup>∂</sup>x<sup>i</sup>* (*<sup>X</sup>* (*a*, *<sup>t</sup>*), *<sup>t</sup>*) (15)

*<sup>∂</sup>x<sup>i</sup>* (*<sup>X</sup>* (*a*, *<sup>t</sup>*), *<sup>t</sup>*) . (16)

= 1 , (17)

*<sup>∂</sup>a<sup>j</sup>* (18)

here *p* is the pressure field. Now the equations for evolution of the vorticity *ω* = ∇ × **u** field

$$D\_l \omega = \omega \cdot \nabla \mathbf{u}\_{\prime} \tag{2}$$

where *Dt* is the convective derivative

$$D\_t = \partial\_t + \mathbf{u} \cdot \nabla \,. \tag{3}$$

A well known consequence of this equation is the preservation of vorticity lines (Helmholtz's theorem).

Here we introduce Clebsch variables (Lamb, 1932). They can be considered as a representation of vorticity lines. In fact from this transformation, which defines the velocity field in terms of scalar variables (*λ*, *μ*, *φ*)

$$\mathbf{u} = \lambda \nabla \mu - \nabla \phi,\tag{4}$$

we can write the vorticity field as

$$
\omega = \nabla \times \mathbf{u} = \nabla \lambda \times \nabla \mu. \tag{5}
$$

Vorticity lines **r**(*s*) are defined as the solutions of

$$\frac{d\mathbf{r}}{ds} = \omega(\mathbf{r}(\mathbf{s})) \,\,\,\,\,\,\tag{6}$$

which admits integrals

$$
\lambda(\mathbf{r}(\mathbf{s})) = \text{const.} \tag{7}
$$

$$
\mu(\mathbf{r}(\mathbf{s})) = \text{const.}
$$

In other words the intersections of surfaces *λ* =const. and *μ* =const. are the vorticity lines. If vorticity lines follow Euler equations and are preserved, then the fields *λ* and *μ* follow the fluid.

Clebsch variables can also be used to find a variational principle for Euler equations. We can write a Lagrangian density for Euler equations

<sup>L</sup> <sup>=</sup> <sup>|</sup>**u**<sup>|</sup> 2 <sup>2</sup> <sup>+</sup> *λ∂t<sup>μ</sup>* , (8)

and the variations of L in function of the fields *λ*, *μ* and *φ* give us the system equations

$$\begin{aligned} \frac{\delta \mathcal{L}}{\delta \mu} &= -D\_t \lambda = 0 \\ \frac{\delta \mathcal{L}}{\delta \lambda} &= \begin{array}{c} D\_t \mu \end{array} = 0 \\ \frac{\delta \mathcal{L}}{\delta \phi} &= \nabla \cdot \mathbf{u} \ = 0. \end{aligned}$$

From this system and the identity

$$\mathbb{E}\left[\nabla\_{\prime}D\_{t}\right] \equiv (\nabla \mathbf{u}) \cdot \nabla\_{\prime} \tag{10}$$

we can obtain the evolution equation for **u**

$$D\_l \mathbf{u} = -\nabla \left( D\_l \phi + \frac{1}{2} \mathbf{u}^2 \right) . \tag{11}$$

#### **2.2 Weber transformation**

2 Will-be-set-by-IN-TECH

here *p* is the pressure field. Now the equations for evolution of the vorticity *ω* = ∇ × **u** field

A well known consequence of this equation is the preservation of vorticity lines (Helmholtz's

Here we introduce Clebsch variables (Lamb, 1932). They can be considered as a representation of vorticity lines. In fact from this transformation, which defines the velocity field in terms of

*d***r**

*μ*(**r**(*s*)) = const.

In other words the intersections of surfaces *λ* =const. and *μ* =const. are the vorticity lines. If vorticity lines follow Euler equations and are preserved, then the fields *λ* and *μ* follow the

Clebsch variables can also be used to find a variational principle for Euler equations. We can

2

*δλ* <sup>=</sup> *Dt<sup>μ</sup>* <sup>=</sup> <sup>0</sup>

*δφ* <sup>=</sup> ∇ · **<sup>u</sup>** <sup>=</sup> 0 .

<sup>L</sup> <sup>=</sup> <sup>|</sup>**u**<sup>|</sup>

*δ*L

*δ*L

*δ*L

and the variations of L in function of the fields *λ*, *μ* and *φ* give us the system equations

where *Dt* is the convective derivative

theorem).

scalar variables (*λ*, *μ*, *φ*)

which admits integrals

fluid.

we can write the vorticity field as

Vorticity lines **r**(*s*) are defined as the solutions of

write a Lagrangian density for Euler equations

From this system and the identity

*Dtω* = *ω* · ∇**u** , (2)

*Dt* = *∂<sup>t</sup>* + **u** · ∇ . (3)

**u** = *λ*∇*μ* − ∇*φ*, (4)

*ds* <sup>=</sup> *<sup>ω</sup>*(**r**(*s*)), (6)

*λ*(**r**(*s*)) = const. (7)

<sup>2</sup> <sup>+</sup> *λ∂t<sup>μ</sup>* , (8)

*δμ* <sup>=</sup> <sup>−</sup>*Dt<sup>λ</sup>* <sup>=</sup> <sup>0</sup> (9)

*ω* = ∇ × **u** = ∇*λ* × ∇*μ*. (5)

Let us note *a<sup>i</sup>* as the initial coordinate (at *t* = 0) of a fluid element and *X<sup>i</sup>* (*a*, *t*) its position at time *t* and note *A<sup>i</sup>* (*x*, *<sup>t</sup>*) the inverse application: *<sup>a</sup><sup>i</sup>* <sup>≡</sup> *<sup>A</sup><sup>i</sup>* (*X<sup>i</sup>* (*a*, *t*), *t*). At time *t* Eulerian coordinates are by definition the variables *x<sup>i</sup>* = *X<sup>i</sup>* (*a*, *t*) then the Lagrangian velocity of a fluid element is

$$
\tilde{\mu}^i(a,t) = \frac{\partial X^i}{\partial t} \left(a\_\prime t\right) \tag{12}
$$

and its acceleration

$$\frac{\partial \vec{u}^i}{\partial t}(a, t) = \frac{\partial^2 X^i}{\partial t^2}(a, t) \; . \tag{13}$$

Newton equations for the fluid element are

$$\frac{\partial^2 X^i}{\partial t^2} \left( a, t \right) = F\_X^i \left( a, t \right) \,, \tag{14}$$

where the forces *F<sup>i</sup> <sup>X</sup>* (*a*, *t*) are given by

$$F\_X^i{}\left(a,t\right) = -\frac{\partial p}{\partial \mathbf{x}^i}\left(X\left(a,t\right),t\right) \tag{15}$$

and *p* (*X* (*a*, *t*)) is the pressure field in Eulerian coordinates. Therefore the movement equations for the fluid elements are

$$\frac{\partial^2 X^i}{\partial t^2} \left( a, t \right) = -\frac{\partial p}{\partial x^i} \left( X \left( a, t \right), t \right) \,. \tag{16}$$

For an incompressible fluid, the transformation matrix, between Lagrangian and Eulerian coordinates, verifies

$$\det\left(\frac{\partial X^i}{\partial a^j}\right) = 1 \,, \tag{17}$$

this value is fixed from the relation between the volume elements in the two coordinate systems. We also note that this transformation is always invertible.

From Eq. (16) we perform a coordinate transformation for the derivatives of *p* using

$$\begin{aligned} \frac{\partial}{\partial \mathbf{x}^i} &= \frac{\partial A^j}{\partial \mathbf{x}^i} \frac{\partial}{\partial a^j} \\ \frac{\partial}{\partial a^i} &= \frac{\partial X^j}{\partial a^i} \frac{\partial}{\partial \mathbf{x}^j} \end{aligned} \tag{18}$$

*∂λ*˜ *<sup>i</sup>*

*∂μ*˜*<sup>i</sup>*

*∂φ*˜ *∂t*

If we now go to the Eulerian coordinates, identifying *μ<sup>i</sup>*

<sup>0</sup> (*μ*(*x*,*t*)), we obtain the *Weber-Clebsch transformation*

be written in Eulerian coordinates as

*ui*

potential pairs.

*h* =

terms, we integrate by parts

*V*

∇ · *ω* = 0 and therefore we have

*V*

<sup>∇</sup>*<sup>φ</sup>* · (∇*<sup>λ</sup>* × ∇*μ*) *<sup>d</sup>*3*<sup>x</sup>* <sup>=</sup>

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) = <sup>0</sup>

*ui*

*Dtλ<sup>i</sup>*

*Dtμ<sup>i</sup>*

If we use Clebsch variables Eq. (4) to represent the velocity field **u**

we have the problem that **u** is restricted to fields with mean helicity

(*x*, *t*) =

(*x*, *t*) = 0

*Dtφ*(*x*, *<sup>t</sup>*) = *<sup>p</sup>*(*x*, *<sup>t</sup>*) <sup>−</sup> <sup>1</sup>

*h* =

of value zero (Grossmann, 1975). In fact, writing *h* in terms of Clebsch variables

*∂V*

(*λ*∇*<sup>μ</sup>* − ∇*φ*) · (∇*<sup>λ</sup>* × ∇*μ*) *<sup>d</sup>*3*<sup>x</sup>* <sup>=</sup> <sup>−</sup>

*V*

the term *λ*∇*μ* is perpendicular to ∇*λ* × ∇*μ* and then their scalar product is zero. For the other

but in a periodic domain the first term in the right hand side is zero. We also know that

*φ* (∇*λ* × ∇*μ*) · *d***s** −

(*a*, *<sup>t</sup>*) = *<sup>p</sup>*˜ (*a*, *<sup>t</sup>*) <sup>−</sup> <sup>1</sup>

3 ∑ *j*=1

Using the convective derivative, the dynamic equations for the Clebsch variables Eq. (27) can

The Weber-Clebsch transformation Eq. (28) and its evolution laws Eq. (29) are very similar to Clebsch variables Eq. (4) and the system (9). An important difference is the number of

2

Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations 113

*<sup>λ</sup><sup>j</sup> ∂μ<sup>j</sup> <sup>∂</sup>x<sup>i</sup>* <sup>−</sup> *∂φ*

> 2 *ui* (*x*, *t*) 2 .

*V*

*V*

*∂X<sup>i</sup> <sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*)

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) = <sup>0</sup> (27)

  2 .

(*x*, *t*) = *A<sup>i</sup>*

(*x*, *t*) = 0 (29)

**u** = *λ*∇*μ* − ∇*φ* , (30)

**<sup>u</sup>** · *<sup>ω</sup>d*3*<sup>x</sup>* (31)

<sup>∇</sup>*<sup>φ</sup>* · (∇*<sup>λ</sup>* × ∇*μ*) *<sup>d</sup>*3*<sup>x</sup>* , (32)

*φ*∇ · (∇*λ* × ∇*μ*) (33)

(*x*, *t*) and *λ<sup>i</sup>*

*<sup>∂</sup>x<sup>i</sup>* . (28)

(*x*, *t*) =

to obtain

$$\frac{\partial^2 X^i}{\partial t^2} \left( a, t \right) = -\frac{\partial A^j}{\partial x^i} \frac{\partial \tilde{p}}{\partial a^j} \left( a, t \right) \tag{19}$$

where *p*˜ (*a*, *t*) is the pressure field in Lagrangian coordinates.

We multiply Eq. (19) with the inverse coordinate transformation *<sup>∂</sup>X<sup>i</sup> <sup>∂</sup>aj* in order to obtain

$$\frac{\partial^2 X^i}{\partial t^2} \left( a, t \right) \frac{\partial X^i}{\partial a^j} \left( a, t \right) = -\frac{\partial \vec{p}}{\partial a^j} \left( a, t \right) \tag{20}$$

which is the *Lagrangian form* for the dynamic equations. The left hand side of this equation can be written as

$$\frac{\partial^2 X^i}{\partial t^2} \left( a, t \right) \frac{\partial X^i}{\partial a^j} \left( a, t \right) = \frac{\partial}{\partial t} \left( \frac{\partial X^i}{\partial t} \left( a, t \right) \frac{\partial X^i}{\partial a^j} \left( a, t \right) \right) - \frac{1}{2} \frac{\partial}{\partial a^j} \left| \frac{\partial X^i}{\partial t} \left( a, t \right) \right|^2 \tag{21}$$

and Eq. (20) becomes

$$\frac{\partial}{\partial t} \left( \frac{\partial X^i}{\partial t} \left( a, t \right) \frac{\partial X^i}{\partial a^j} \left( a, t \right) \right) = -\frac{\partial \vec{q}}{\partial a^j} \left( a, t \right) \tag{22}$$

where the term *q*˜(*a*, *t*) is given by

$$
\tilde{q}\left(a,t\right) = \tilde{p}\left(a,t\right) - \frac{1}{2} \left|\frac{\partial X^i}{\partial t}\left(a,t\right)\right|^2. \tag{23}
$$

Now let us integrate Eq. (22) over *t*, maintaining *a<sup>i</sup>* fixed

$$\left[\frac{\partial X^i}{\partial t}\left(a,t\right)\frac{\partial X^i}{\partial a^j}\left(a,t\right)\right]\_0^t = \frac{\partial X^i}{\partial t}\left(a,t\right)\frac{\partial X^i}{\partial a^j}\left(a,t\right) - \tilde{u}\_0^j(a)$$

to obtain

$$\frac{\partial X^{i}}{\partial t}\left(a\_{\prime}t\right)\frac{\partial X^{i}}{\partial a^{j}}\left(a\_{\prime}t\right) - \tilde{u}\_{0}^{j}\left(a\right) = -\frac{\partial \tilde{\phi}}{\partial a^{j}}\left(a\_{\prime}t\right) \,, \tag{24}$$

where *φ*˜ is written as

$$\tilde{\phi}\left(a,t\right) = \int\_0^t \left(\tilde{\rho}\left(a,s\right) - \frac{1}{2} \left|\frac{\partial X^i}{\partial t}\left(a,s\right)\right|^2\right)ds\,. \tag{25}$$

This equation system (24) is called *Weber transformation* (Lamb, 1932). Now we perform a coordinate transformation

$$\frac{\partial X^{i}}{\partial t}\left(a\_{\prime}t\right) = \mathfrak{u}\_{0}^{j}(a)\frac{\partial A^{j}}{\partial x^{i}} - \frac{\partial A^{j}}{\partial x^{i}}\frac{\partial \tilde{\Phi}}{\partial a^{j}}\left(a\_{\prime}t\right) \,. \tag{26}$$

identifying *μ*˜*<sup>i</sup>* (*a*, *t*) = *a<sup>i</sup>* and the initial velocity *λ*˜ *<sup>i</sup>* (*a*, *t*) = *u*˜*<sup>i</sup>* <sup>0</sup> (*a*) we obtain the evolution equations for the fields in Lagrangian coordinates

4 Will-be-set-by-IN-TECH

*∂x<sup>i</sup>*

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*) <sup>=</sup> <sup>−</sup> *<sup>∂</sup>p*˜

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) *<sup>∂</sup>X<sup>i</sup>*

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*)

2 � � � � �

<sup>=</sup> *<sup>∂</sup>X<sup>i</sup>*

*j*

2 � � � � � *∂p*˜

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*)

<sup>=</sup> <sup>−</sup> *<sup>∂</sup>q*˜

�

*∂X<sup>i</sup> <sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*)

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) *<sup>∂</sup>X<sup>i</sup>*

<sup>0</sup>(*a*) = <sup>−</sup> *∂φ*˜

*∂X<sup>i</sup> <sup>∂</sup><sup>t</sup>* (*a*,*s*)

� − 1 2 *∂ ∂a<sup>j</sup>*

> � � � � �

2

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*) <sup>−</sup> *<sup>u</sup>*˜

� � � � �

*∂φ*˜

(*a*, *t*) = *u*˜*<sup>i</sup>*

2 ⎞ *j* <sup>0</sup>(*a*)

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*) (19)

*<sup>∂</sup>aj* in order to obtain

� � � � � 2

(21)

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*) (20)

� � � � �

*∂X<sup>i</sup> <sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*)

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*) (22)

. (23)

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*) , (24)

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*) . (26)

<sup>0</sup> (*a*) we obtain the evolution

⎠ *ds*. (25)

*<sup>∂</sup>t*<sup>2</sup> (*a*, *<sup>t</sup>*) <sup>=</sup> <sup>−</sup> *<sup>∂</sup>A<sup>j</sup>*

*∂*2*X<sup>i</sup>*

where *p*˜ (*a*, *t*) is the pressure field in Lagrangian coordinates. We multiply Eq. (19) with the inverse coordinate transformation *<sup>∂</sup>X<sup>i</sup>*

which is the *Lagrangian form* for the dynamic equations. The left hand side of this equation can be written as

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*) <sup>=</sup> *<sup>∂</sup>*

� *∂X<sup>i</sup>*

*∂ ∂t*

Now let us integrate Eq. (22) over *t*, maintaining *a<sup>i</sup>* fixed

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) *<sup>∂</sup>X<sup>i</sup>*

*∂X<sup>i</sup>*

*φ*˜ (*a*, *t*) =

*∂X<sup>i</sup>*

Now we perform a coordinate transformation

equations for the fields in Lagrangian coordinates

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) *<sup>∂</sup>X<sup>i</sup>*

� *t* 0

This equation system (24) is called *Weber transformation* (Lamb, 1932).

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) <sup>=</sup> *<sup>u</sup>*˜

(*a*, *t*) = *a<sup>i</sup>* and the initial velocity *λ*˜ *<sup>i</sup>*

⎛

*∂*2*X<sup>i</sup>*

where the term *q*˜(*a*, *t*) is given by

� *∂X<sup>i</sup>*

and Eq. (20) becomes

to obtain

where *φ*˜ is written as

identifying *μ*˜*<sup>i</sup>*

*<sup>∂</sup>t*<sup>2</sup> (*a*, *<sup>t</sup>*) *<sup>∂</sup>X<sup>i</sup>*

*∂*2*X<sup>i</sup>*

*<sup>∂</sup>t*<sup>2</sup> (*a*, *<sup>t</sup>*) *<sup>∂</sup>X<sup>i</sup>*

*∂t*

� *∂X<sup>i</sup>*

*<sup>∂</sup><sup>t</sup>* (*a*, *<sup>t</sup>*) *<sup>∂</sup>X<sup>i</sup>*

*<sup>q</sup>*˜(*a*, *<sup>t</sup>*) <sup>=</sup> *<sup>p</sup>*˜ (*a*, *<sup>t</sup>*) <sup>−</sup> <sup>1</sup>

�*t*

0

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*) <sup>−</sup> *<sup>u</sup>*˜

<sup>⎝</sup>*p*˜ (*a*,*s*) <sup>−</sup> <sup>1</sup>

*j* <sup>0</sup>(*a*) *∂A<sup>j</sup> <sup>∂</sup>x<sup>i</sup>* <sup>−</sup> *<sup>∂</sup>A<sup>j</sup> ∂x<sup>i</sup>*

*<sup>∂</sup>a<sup>j</sup>* (*a*, *<sup>t</sup>*)

to obtain

$$\begin{aligned} \frac{\partial \bar{X}^i}{\partial t}(a, t) &= 0\\ \frac{\partial \bar{\mu}^i}{\partial t}(a, t) &= 0 \\ \frac{\partial \tilde{\phi}}{\partial t}(a, t) &= \vec{\mu} \left(a, t\right) - \frac{1}{2} \left| \frac{\partial X^i}{\partial t} \left(a, t\right) \right|^2 . \end{aligned} \tag{27}$$

If we now go to the Eulerian coordinates, identifying *μ<sup>i</sup>* (*x*, *t*) = *A<sup>i</sup>* (*x*, *t*) and *λ<sup>i</sup>* (*x*, *t*) = *ui* <sup>0</sup> (*μ*(*x*,*t*)), we obtain the *Weber-Clebsch transformation*

$$u^{i}(\mathbf{x},t) = \sum\_{j=1}^{3} \lambda^{j} \frac{\partial \mu^{j}}{\partial \mathbf{x}^{i}} - \frac{\partial \phi}{\partial \mathbf{x}^{i}} \,. \tag{28}$$

Using the convective derivative, the dynamic equations for the Clebsch variables Eq. (27) can be written in Eulerian coordinates as

$$\begin{aligned} D\_l \lambda^i(\mathbf{x}, t) &= 0 \\ D\_l \mu^i(\mathbf{x}, t) &= 0 \\ D\_l \phi(\mathbf{x}, t) &= p(\mathbf{x}, t) - \frac{1}{2} \left| \mu^i(\mathbf{x}, t) \right|^2. \end{aligned} \tag{29}$$

The Weber-Clebsch transformation Eq. (28) and its evolution laws Eq. (29) are very similar to Clebsch variables Eq. (4) and the system (9). An important difference is the number of potential pairs.

If we use Clebsch variables Eq. (4) to represent the velocity field **u**

$$\mathbf{u} = \lambda \nabla \mu - \nabla \phi \,\tag{30}$$

we have the problem that **u** is restricted to fields with mean helicity

$$h = \int\_{V} \mathbf{u} \cdot \boldsymbol{\omega} d^3 \mathbf{x} \tag{31}$$

of value zero (Grossmann, 1975). In fact, writing *h* in terms of Clebsch variables

$$h = \int\_{V} \left(\lambda \nabla \mu - \nabla \phi\right) \cdot \left(\nabla \lambda \times \nabla \mu\right) d^3 \mathbf{x} = -\int\_{V} \nabla \phi \cdot \left(\nabla \lambda \times \nabla \mu\right) d^3 \mathbf{x} \,, \tag{32}$$

the term *λ*∇*μ* is perpendicular to ∇*λ* × ∇*μ* and then their scalar product is zero. For the other terms, we integrate by parts

$$\int\_{V} \nabla \phi \cdot (\nabla \lambda \times \nabla \mu) \, d^3 \mathbf{x} = \oint\_{\partial V} \phi \left( \nabla \lambda \times \nabla \mu \right) \cdot d\mathbf{s} - \int\_{V} \phi \nabla \cdot (\nabla \lambda \times \nabla \mu) \tag{33}$$

but in a periodic domain the first term in the right hand side is zero. We also know that ∇ · *ω* = 0 and therefore we have

*∂ <sup>∂</sup>a<sup>i</sup>* <sup>=</sup>

*∂ <sup>∂</sup>x<sup>i</sup>* <sup>=</sup>

> *∂ <sup>∂</sup>x<sup>i</sup>* , *<sup>∂</sup> ∂x<sup>k</sup>*

 *∂ <sup>∂</sup>a<sup>i</sup>* , *<sup>∂</sup> ∂a<sup>k</sup>* = 0 .

the original Lagrangian position *μ*, we can express the commutator Eq. (43) as

Γ*m*

over a vector or a matrix each component is taken in an independent way.

Lagrangian coordinates

and *∂a*

formula

operator

We also have the relation for the commutators

 *∂ <sup>∂</sup>a<sup>i</sup>* , *<sup>∂</sup> ∂x<sup>k</sup>* = *∂ <sup>∂</sup>a<sup>i</sup>* ,

> *∂ <sup>∂</sup>a<sup>i</sup>* , *<sup>∂</sup> ∂x<sup>k</sup>* <sup>=</sup> *<sup>∂</sup> ∂a<sup>i</sup>*

The term *Cm*,*k*;*<sup>i</sup>* is related to the Christoffel coefficients Γ*<sup>m</sup>*

3 ∑ *m*=1

In the same way it is possible to write the derivatives of the Eulerian coordinate in terms of

Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations 115

Using relations Eq. (40), Eq. (41) and Eq. (42) we can compute the commutators between *∂x*

Introducing the displacement vector *<sup>m</sup>* <sup>=</sup> *<sup>μ</sup><sup>m</sup>* <sup>−</sup> *<sup>x</sup><sup>m</sup>* which relates the Eulerian position *<sup>x</sup>* to

 *∂<sup>m</sup> ∂x<sup>k</sup>*

*ij* <sup>=</sup> <sup>−</sup> *<sup>∂</sup>X<sup>k</sup>*

We consider now the diffusive evolution of our fields, with that goal in mind we define the

where *ν* is the viscosity and **u** is the Eulerian velocity. When the operator Eq. (46) is applied

We also need a coordinate transformation that can be invertible at any time *t*, that condition is always satisfied when the diffusion is zero (*ν* = 0) and the fluid is incompressible, because the fluid element volume is preserved by the coordinate transformation, and therefore

Constantin imposes that the coordinates *μ<sup>i</sup>* are advected and diffused so they follow

*∂*

*<sup>∂</sup>a<sup>m</sup>* <sup>=</sup> *Cm*,*k*;*<sup>i</sup>*

*∂ ∂a<sup>m</sup>* <sup>=</sup> *<sup>∂</sup> ∂a<sup>i</sup>*  *∂μ<sup>m</sup> ∂x<sup>k</sup>*

*∂*

*∂*

*<sup>∂</sup>a<sup>j</sup> Cm*,*k*;*<sup>i</sup>* . (45)

Γ = *∂<sup>t</sup>* + **u** ·∇− *ν*� , (46)

Γ*μ<sup>i</sup>* = 0. (47)

Det(∇*μ*) = 1 , (48)

*∂μ<sup>m</sup> ∂x<sup>k</sup>*

3 ∑ *m*=1 *∂X<sup>m</sup> ∂a<sup>i</sup>*

*∂μ<sup>m</sup> ∂x<sup>i</sup>*

*∂*

*∂*

*<sup>∂</sup>x<sup>m</sup>* . (40)

*<sup>∂</sup>a<sup>m</sup>* . (41)

= 0 (42)

*<sup>∂</sup>a<sup>m</sup>* . (43)

*<sup>∂</sup>a<sup>m</sup>* . (44)

*ij* of the flat connection in **<sup>R</sup>**<sup>3</sup> by the

$$\nabla \cdot (\nabla \lambda \times \nabla \mu) = 0 \tag{34}$$

and we finally get