Three-Dimensional Modeling of Interfacial Tension with Smoothed Particle Hydrodynamics by Using Pairwise Potential

Abstract

A new pairwise potential was introduced into the three-dimensional SPH simulation program and the dynamic behavior of small droplets on a solid substrate was calculated. The simulation could reproduce the experimental observations for an aluminum droplet on an alumina substrate as well as those of water on glass and Teflon substrates. The program was applied to an iron droplet on an aluminum substrate, and the droplet’s dynamic behavior was estimated.

1. Introduction

Direct observation of the metallurgical reactions occurring in a steelmaking process is quite difficult because of the extremely high-temperature environment and invisible nature of the materials. Therefore, it is necessary to reproduce and predict the behavior of molten steel and slag by using computer simulations for the development of sophisticated steelmaking processes. The modeling of interfacial tensions is highly important in the estimation of slag and metal deformations. Particle methods are promising simulation models because they can treat large deformations in fluids without the problem of grid tangling. The most popular particle methods used in the metallurgical field are the moving particle semi-implicit (MPS) and smoothed particle hydrodynamics (SPH) methods. Because the former was developed for incompressible fluids, an implicit algorithm is used in the pressure calculation to satisfy the incompressibility.¹⁾ On the other hand, the latter method was developed for compressible fluids, and an explicit algorithm is used, where the incompressibility is approximated using the special equation of state.^2,3) The calculation speed in SPH is usually much faster than that in MPS because of the algorithms employed. Because the number of particles increases dramatically in three-dimensional simulations, the SPH method has the advantage of reduced calculation costs.

In this study, a new pairwise potential was introduced into the SPH method, and a three-dimensional simulation model was developed. The model was compared with experimental observations of a water droplet and an aluminum droplet on solid substrates. The model was applied to an iron droplet on an alumina substrate, and its dynamic behavior was calculated.

2. Simulation Model

2.1. SPH Method^2,3)

When a particle (calculation point) is located in space at a position x, the physical properties of the particle f(x) such as a velocity and a pressure can be described by Eq. (1), where Ω is an integral region, and W is a kernel function.   

$$f\left(x\right)={\int }_{\mathrm{\Omega }}f\left(x^{\prime }\right)W\left(x-x^{\prime },h\right)d{x}^{\prime }$$(1)

If the physical properties of a particle are distributed continuously in a sphere of radius h around the particle, the distribution can be described by the kernel function W(r, h) defined in Eq. (2), where r is the distance between particles, and h is the effective radius of the kernel function.   

$$W\left(r,h\right)=\omega \left(\frac{5}{\pi h^2}\right)\left(1+3\frac{r}{h}\right){\left(1-\frac{r}{h}\right)}^3$$(2)

Equation (1) can be discretized to the summation form expressed in Eq. (3), where n is the number of particles within the effective radius, x_i is the position of the i-th particle and m_i and ρ_i are its mass and density.   

$$f\left(x\right)={\sum }_{i=1}^n\frac{m_i}{{\rho }_i}f\left(x-x_i\right)W\left(x-x_i,h\right)$$(3)

The governing equation is the Navier–Stokes equation, where ν is the velocity, ρ is the density, P is the pressure, μ is the viscosity, and F is an external force. The gravity and interfacial force are included in F.   

$$\frac{Dv}{Dt}=-\frac{1}{\rho }\nabla P+ \frac{\mu }{\rho }\mathrm{\Delta }v+F$$(4)

Since the space derivatives of a physical property can be expressed by using those of the kernel function, Eq. (4) was discretized in terms of Eqs. (2) and (3) for the numerical calculation.

2.2. New Interfacial Tension Model

Tarakovsky et al.⁴⁾ introduced the interparticle potential into the SPH method and reported that the wettability was described reasonably well for a water droplet. They developed a two-dimensional model and proposed the interparticle potential force F, which is repulsive at short distances and attractive at long distances. The force F is given by Eq. (5), where S is the strength coefficient of the force acting between particles.   

$$F=\{\begin{array}{c}S\cdot cos\left(\frac{1.5\pi r}{3h}\right)\hspace{1em}\hspace{1em}\left(r\le h\right)\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(r>h\right)\end{array}$$(5)

However, the theoretical relation between the strength coefficient S and the interfacial tension was not described.

Hongo et al.⁵⁾ proposed a different model for SPH simulation. The interfacial tension Γ was expressed by the surface tensions of both phases, as shown in Eq. (6), where Γ_A is the surface tension of phase A, Γ_B is that of phase B, and α is an adjustable parameter that has a poorly understood physical meaning.   

$$\mathrm{\Gamma }=\alpha {\mathrm{\Gamma }}_A+\left(1-\alpha \right){\mathrm{\Gamma }}_B\hspace{1em}\left(0\le \alpha \le 1\right)$$(6)

Kondo et al.⁶⁾ proposed the pairwise potential given by Eq. (7) and applied it to a two-dimensional simulation using the MPS method, where C is a coefficient, r_min is the initial nearest-particle distance, and r_s is the effective radius of the potential.   

$$P\left(r<r_s\right)=\frac{1}{3}C\left(r-\frac{3}{2}{r}_{min}+\frac{1}{2}{r}_s\right){\left(r-r_s\right)}^2$$(7)

The interparticle force F is obtained as the first derivative of the potential P by r, which is repulsive at short distances and attractive at long distances.   

$$F\left(r<r_s\right)=C\left(r-r_{min}\right)\left(r-r_s\right)$$(8)

The relation between the interparticle force and the interfacial energy was theoretically derived, and the coefficient C can be calculated from the surface tensions and equilibrium contact angle. Natsui et al.⁷⁾ developed Kondo’s method and applied it to a three-dimensional simulator using the MPS method. They compared the simulation results with experiments on the dynamic behavior of a water droplet and reported good agreement.

In this study, the pairwise potential proposed by Kondo et al. [Eq. (7)] was applied to a three-dimensional SPH simulation instead of Tarakovsky et al.’s interparticle potential force F [Eq. (5)], because Eq. (7) rests on the theoretical background of surface thermodynamics, and its much simpler function is preferable for the calculation of a large number of particles to reduce the calculation cost.

When Eqs. (7) and (8) were directly used in the SPH simulation, the accumulated calculation error disturbed the stable computation, especially for small scales such as a droplet having a diameter less than 0.01 m. This is probably because the interparticle force calculated by Eq. (8) becomes a product of very small numbers as the simulation scale decreases. In addition, the effective radius of the potential r_s is independent of the effective radius for the Navier–Stokes equation, and the calculation error might be minimized by adjusting r_s in the MPS method. On the other hand, in the SPH method, the effective radius of the potential should be equal to the effective radius of the kernel function in principle,³⁾ and accumulation of calculation errors is inevitable because an explicit algorithm is used.

To solve these problems, the new pairwise potential φ using a dimensionless distance q (=r/h) was introduced, as shown in Eq. (9), where d is the initial particle distance, and E is the energy coefficient for liquid–liquid particles.   

$$\phi \left(q\le 1\right)=E\left(q-\frac{h-3d}{2h}\right){\left(q-1\right)}^2$$(9)

According to the manner proposed by by Kondo et al.,⁶⁾ The coefficient E is given by Eq. (10).   

$$E=\frac{2\sigma d^2}{h^3\mathrm{\Sigma }\phi \left(\text{q}\right)}$$(10)

For the determination of E, the system composed of the particles within the hemisphere with radius h (region A) and the particles aligned on the normal line with length h at the center of its great circle (region B) was considered. The summation of the pairwise potential between the particle in region A and that in region B in the system was calculated. Then the value of E was determined from the available surface tension data of a liquid using Eq. (10) for the given values of d and h. The coefficient for solid–liquid particles expressed as E_inter is also given by Eq. (11),⁶⁾ where θ is the equilibrium contact angle of a liquid on a solid substrate which is used as the wettability parameter in this model.   

$$E_{inter}=\frac{1}{2}\left(cos\theta +1\right)E$$(11)

By substituting Eq. (9) into Eq. (11), differentiating by q and summing up the interparticle forces, the interfacial force vector F_inter is expressed as Eq. (12), where an arrow on r indicates a vector. Equation (12) is substituted into the external force term in Eq. (4).^2,8)   

$$F_{inter}=\frac{1}{2}\left(cos\theta +1\right)\cdot \mathrm{\Sigma }\frac{\text{d}\phi \left(\text{q}\right)}{\text{dq}}\frac{\overrightarrow{r}}{\left|r\right|}$$(12)

The calculation was stabilized and its error was significantly reduced by using Eqs. (9) and (12) even if the simulation scale was reduced.

3. Three-dimensional Simulation of a Droplet

3.1. Simulation of a Water Droplet on a Glass and a Teflon Substrates

To verify the developed simulation model, the dynamic behavior of a water droplet was experimentally observed under the conditions reported by Natsui et al.⁷⁾ and was compared with the calculation. As mentioned before, an equilibrium contact angle θ is used as the wettability parameter in this model. The simulation conditions are tabulated⁹⁾ in Table 1. The effective radius h was set to 4.0d for the calculation stability. Parallel processing using an Intel Xenon X5675 (3.07 GHz, 6 cores)×2 was performed to reduce the calculation time.

Table 1. Simulation conditions.

	water	Al	Fe
Temperature (K)	298	1373	1873
Density (kg/m3)9)	997	2240	7020
Viscosity (mPa·s)9)	0.089	2.56	4.8
Surface tension (mN/m)9)	73	900	1800
Equilibrium Contact angle (deg)	30 (glass) 90 (Teflon)	90	125
Initial inter-particle distance (m)	1.0×10−4	2.0×10−4	2.0×10−4
The number of liquid particles	4140	4160	8144
Initial diameter of droplet (m)	2.0×10−3	3.92×10−3	5.0×10−3

A water droplet 2 mm in diameter was dropped onto solid substrates of glass (commercial slide glass, θ=30°) and Teflon (commercial Teflon plate, θ=90°) from 2 mm above it using a microsyringe, and its behavior was observed by high-speed video. The surface of the plates were carefully cleaned, however, the surface roughness was not measured in the experiment and its effect on the droplet behavior was not considered in the simulation model. The simulation showed that the water droplet reached the equilibrium contact angle after attenuation of the oscillation, which showed good agreement with the observation. The same calculations have been already reported in detail by Natsui et al.⁷⁾

The motion of the same droplet on a substrate inclined by 60° was observed. Figure 1 shows a snapshot of the moving droplet, where the left panel shows a side view of the observation, the middle panel is that of the simulation, and the right panel is a bird’s-eye view of the simulation. The shapes of the droplet were reproduced reasonably well for both substrates in the simulation. However, the advancing and receding contact angles were slightly overestimated and underestimated. The inaccuracy of the force calculation might be caused by the decrease in the number of particles near the edge of the droplet, which would be solved by increasing the number of particles using the faster calculation method.

Fig. 1.

Comparison of the calculated shapes of a moving water droplet on an inclined substrate with the experimental results: a) on a glass plate, b) on a Teflon plate.

3.2. Simulation of a Liquid Aluminum Droplet on an Alumina Substrate

Fujii et al. studied the dynamic wetting of aluminum on solid substrates by an in-situ observation technique.^10,11) The calculation was performed for their experimental conditions to check the applicability of the developed model to actual molten metal systems. A 70 mg aluminum droplet was dropped onto an alumina substrate from 4 mm above it. The droplet deformed and approached equilibrium at 100 ms after impact, as shown in Fig. 2. The image on the lower right shows a side view of the aluminum droplet in the experiment, which was obtained by tracing the picture at 100 ms observed by Fujii.¹¹⁾ The calculated shape and equilibrium contact angle showed good agreement with those obtained in their experiment.

Fig. 2.

Simulated snapshot of 70 mg aluminum droplet dropped onto an alumina substrate at 1373 K (θ = 90°).

3.3. Simulation of a Liquid Iron Droplet on an Alumina Substrate

The program was applied to an iron droplet on an aluminum substrate, and the droplet’s dynamic behavior was estimated. Figure 3 shows snapshots of a droplet 5 mm in diameter dropped from 5 mm above the substrate at 1873 K. The iron droplet oscillated twice and reached equilibrium at 110 ms after impact. Figure 4 shows the time variation of the dynamic contact angle and height of the droplet. The contact angle was calculated by averaging the right and left ones in the side view. It oscillated between 117° and 145° and reached an equilibrium value of 125°. The droplet height varied between 1.5 and 4.0 mm and also reached an equilibrium value, 3.2 mm.

Fig. 3.

Simulated snapshot of iron droplet 5 mm in diameter dropped onto an alumina substrate at 1873 K (θ = 125°).

Fig. 4.

Variation of the dynamic contact angle and height of the iron droplet with time.

4. Conclusion

A new pairwise potential was introduced into the three-dimensional SPH simulation program, and the dynamic behavior of small droplets on a solid substrate was calculated. The simulation could reproduce the experimental observations for an aluminum droplet on an alumina substrate as well as those of water on glass and Teflon substrates. The program was applied to an iron droplet on an aluminum substrate, and the droplet’s dynamic behavior was estimated. The dynamic contact angle oscillated between 117° and 145° and reached an equilibrium value of 125°.

References

• 1)   M.  Asai,  H.  Nijo and  K.  Ito: ISIJ Int., 49 (2009), 178.
• 2)   K.  Ito: Bull. Iron Steel Inst. Jpn., 19 (2014), 809.
• 3)   G. R.  Liu and  M. B.  Liu: Smoothed Particla Hydrodynamics, World Scientific, Singapore, (2003), 35.
• 4)   A.  Tarakovsky and  P.  Meakin: Phys. Rev. E, 72 (2005), 26301.
• 5)   T.  Hongo,  M.  Shigeta,  S.  Izawa and  Y.  Fukunishi: Proc. 23rd CFD Symp. JSFM, Tokyo, (2009), A8-5.
• 6)   M.  Kondo,  S.  Koshizuka and  M.  Takimoto: Trans. JSCES, 11 (2007), 20070021.
• 7)   S.  Natsui,  R.  Soda,  T.  Kon,  S.  Ueda,  J.  Kano,  R.  Inoue and  T.  Ariyama: Mater. Trans., 53 (2012), 662.
• 8)   M.  Nakano and  K.  Ito: Proc. Asia Steel 2015, ISIJ, Tokyo, (2015), 598.
• 9)   Y.  Kawai and  Y.  Shiraishi: Handbook of Physico-chemical Properties at High Temperatures, ISIJ, Tokyo, (1988).
• 10)   H.  Fujii,  H.  Nakae and  K.  Okada: Acta Metall. Mater., 41 (1993), 2963.
• 11)   H.  Fujii: Ph.D. thesis, Waseda University, (1993), http://dspace.wul.waseda.ac.jp/dspace/handle/2065/41038, (accessed 2015-11-07).