Ligament break-up simulation through pseudo-potential Lattice Boltzmann Method

TL;DR

This paper presents a numerical study of ligament breakup due to Plateau-Rayleigh instability using a pseudo-potential lattice Boltzmann method, highlighting its potential for engineering jet hydrodynamics applications.

Contribution

It introduces a systematic analysis of ligament fragmentation using SRT pseudo-potential LBM, demonstrating convergence and applicability to engineering problems.

Findings

Simulations show asymptotic behavior consistent with sharp interface hydrodynamics.

Grid convergence achieved with interface thickness smaller than ligament size.

Method effectively models ligament breakup phenomena.

Abstract

The Plateau-Rayleigh instability causes the fragmentation of a liquid ligament into smaller droplets. In this study a numerical study of this phenomenon based on a single relaxation time (SRT) pseudo-potential lattice Boltzmann method (LBM) is proposed. If systematically analysed, this test case allows to design appropriate parameters sets to deal with engineering applications involving the hydrodynamics of a jet. Grid convergence simulations are performed in the limit where the interface thickness is asymptotically smaller than the characteristic size of the ligament. These simulations show a neat asymptotic behaviour, possibly related to the convergence of LBM diffuse-interface physics to sharp interface hydrodynamics.

Figures

Tables

Table 1. Table 1: Densities (liquid and vapour) and surface tension as a function of ligament radius.

$R_0$	7	14	28	42	56	70	84
${\rho }_l$	2.112	2.117	2.125	2.135	2.145	2.156	2.167
${\rho }_v$	0.081	0.084	0.090	0.097	0.106	0.116	0.126
$\sigma$	0.05159	0.05345	0.05775	0.06277	0.06848	0.07491	0.08206

Equations

$$\frac{\partial f_{\alpha }\left(\text{x},t\right)}{\partial t}+{\text{c}}_{\alpha }\nabla f_{\alpha }\left(\text{x},t\right)=-\frac{1}{\tau }\left[f_{\alpha }\left(\text{x},t\right)-f_{\alpha }^{eq}\left(\text{x},t\right)\right]$$

$$f_{\alpha}^{eq}\left(\textbf{x},t\right)=w_{\alpha}\rho\left(\textbf{x},t\right)\bigg{[}\frac{\textbf{c}_{\alpha}\cdot\textbf{u}\left(\textbf{x},t\right)}{c_{s}^{2}}+\frac{\left[\textbf{c}_{\alpha}\cdot\textbf{u}\left(\textbf{x},t\right)\right]^{2}}{2c_{s}^{4}}-\frac{\left[\textbf{u}\left(\textbf{x},t\right)\cdot\textbf{u}\left(\textbf{x},t\right)\right]}{2c_{s}^{2}}\bigg{]}$$

$${\Delta }_{\alpha }{f}_{\alpha }\left(\text{x},t\right)=f_{\alpha }\left(\text{x}+{\text{c}}_{\alpha }\Delta t,t+\Delta t\right)-f_{\alpha }\left(\text{x},t\right)=-\frac{\Delta t}{\tau }\left[f_{\alpha }\left(\text{x},t\right)-f_{\alpha }^{eq}\left(\text{x},t\right)\right]+F_{\alpha }$$

$$\rho \left(\text{x},t\right)=\sum _{\alpha =0}^{N_{pop}-1}{f}_{\alpha }\left(\text{x},t\right)\rho \text{u}\left(\text{x},t\right)=\sum _{\alpha =0}^{N_{pop}-1}{\text{c}}_{\alpha }{f}_{\alpha }\left(\text{x},t\right)$$

$$\text{F}\left(\text{x},t\right)=-G_0\psi \left(\text{x},t\right)\sum _{\alpha }^{N_{pop}}\psi \left(\text{x}+{\text{c}}_{\alpha }\Delta t,t\right){\text{c}}_{\alpha }{w}_{\alpha }$$

$${\text{u}}^{{}^{\prime }}\left(\text{x},t\right)=\text{u}\left(\text{x},t\right)+\frac{\text{F}\left(\text{x},t\right)\tau }{\rho \left(\text{x},t\right)}$$

$$P=p_0+\frac{c_s^2{G}_0}{2}{\psi }^2=\rho c_s^2+\frac{c_s^2{G}_0}{2}{\psi }^2$$

$$Oh=\frac{\nu }{{\rho }_l\sigma \hat{r}}{t}_{cap}=\sqrt{\frac{{\rho }_l{\hat{r}}^3}{\sigma }}$$

Full text

1 INTRODUCTION

Numerical methods coupled with experimental tests are gaining nowadays great importance for solving a large number of problems which involve both innovative applications (from aerodynamics to multiphase) or innovative materials, [1, 2, 3]. In this scenario, during the recent years an increasing attention has been driven towards the lattice Boltzmann method (LBM) for solving many physical problems [4, 5, 6, 7, 8, 9]. The major appeals pertain its simplicity and applicability in a wide range of conditions, the easy handling of complex geometries (e.g. porous media), the possibility to combine it with other methods to design hybrid approaches [10, 11], and even possible extensions to more complex physics phenomena (e.g. non-Newtonian fluids, [12]). In this ballpark, one of the most common (modern) application, where to exploit all the advantages related to such a method, is multiphase fluid dynamics. LBM, in fact, allows to locally solve the equation of state with apparent benefits for understanding physical phenomena, [13, 14, 15, 16, 17, 18, 19, 20, 21]. Many interesting engineering applications have been already analysed with different multiphase approaches, [22, 23, 24, 25, 26, 16, 27]. The aim of this work is to deeply analyse the Plateau-Rayleigh (PR) instability by means of a single-belt pseudo-potential Lattice Boltzmann method. LBM is by definition a diffuse interface method: at the interface between two different phases, the fluid properties (i.e. density, velocity, pressure) change smoothly over a region with assigned thickness. With regard to the PR instability, this sets a compelling case for conducting a series of numerical simulations in order to understand if an asymptotic behaviour may be recognized when the interface thickness is smaller with respect to the characteristic geometrical length-scale of the physical problem at hand.

2 NUMERICAL MODEL

The Lattice Boltzmann Method represents a relatively new computational approach which may be alternative to standard Navier-Stokes based solvers. It starts from a microscopic kinetic approach in order to reconstruct the macroscopic fluid flow solutions by mathematically describing movements and interactions of the particles which constitute the flow. The LBM [4, 8, 9] starts from the Boltzmann Equation and, after a discrete decomposition, it writes the discrete form of the Boltzmann equation itself for a fixed set of speeds. The discrete Boltzmann equation reads as follows:

[TABLE]

where $f\left(\textbf{x},t\right)$ is the particle distribution function which is generally solved along the allowed velocities directions, [9]. In this paper a $D3Q19$ stencil is used for flow discretization coupled with a uniform Cartesian lattice. The single-relaxation time $\tau$ to the local equilibrium $f_{eq}$ , which is a function of the macroscopic flow variables according to the following equation 2, is used:

[TABLE]

where $w_{\alpha}$ represents a set of weights normalized to unity; $\rho\left(\textbf{x},t\right)$ and $\textbf{u}\left(\textbf{x},t\right)$ are respectively global density and speed and $c_{s}$ is the lattice speed of sound.

Equation 1 may be explicitly solved as follows:

[TABLE]

where $F_{\alpha}$ represents the external forcing - in this case the multiphase interaction described below. Macroscopic fluid density and speed may be derived respectively through the $0^{th}$ and the $1^{st}$ population momentum, as in the following equation:

[TABLE]

where $N_{pop}$ denotes the number of discrete velocities (or populations). As previously anticipated, this work is based on the pseudo-potential forcing which, in case of single-belt formulation [13], may be written as:

[TABLE]

being $\psi\left(\textbf{x},t\right)$ a local functional of a density: $\psi\left(\textbf{x},t\right)=\rho_{0}\left[1-\exp\left(-{\rho\left(\textbf{x},t\right)\over{\rho_{0}}}\right)\right]$. In this application the reference density $\rho_{0}$ is set to $\rho_{0}=1$ and $G_{0}$ is the basic parameter which rules the inter-particle interaction. By the assumption on $\rho_{0}$, the inter-particle strength $G_{0}$ is the only free parameter which fixes both the density ratio and the surface tension. Starting from equation 5, the component of the interaction potential along each direction can be evaluated and then used to shift the macroscopic velocities before evaluating the equilibrium distribution functions:

[TABLE]

Finally, the equation of state of the system may be written as follows:

[TABLE]

3 LIGAMENT SET-UP

In this section the test case and the main parameters are presented. First of all it is useful to define the initial ligament configuration in order to reach its breakup. For the Plateau-Rayleigh instability a sinusoidal perturbation has used with a fixed amplitude and a constant frequency. More specifically the perturbation amplitude is a function of the ligament radius, while the period of oscillation is related to the domain length. With referring to Fig. 1, the radius $r$ function of the $x$ coordinate may be written as: $r\left(x\right)=R_{0}+\frac{R_{0}}{10}\sin\left(\frac{4x}{L}\right)$, where $R_{0}$ is the unperturbed ligament radius and $L$ is the domain length.

The phenomenon is driven by some characteristic parameters, between them the Ohnesorge number $Oh$ and the capillary time $t_{cap}$ defined in the following equation:

[TABLE]

where $\rho_{l}$ is the liquid density, $\nu$ the kinematic viscosity, $\sigma$ the surface tension, $\hat{r}$ the average radius of the ligament along $x$ coordinate. Starting from this configuration and with the parameters above defined, the dimensionless time of break up as a function of initial radius will be analysed.

4 RESULTS

In this section results for different ligament radii are presented. More specifically, the inter-particle strength $G_{0}$ has been fixed to 5.3, while the $Oh$ has been fixed to 0.1 and the ligament radius has been varied from 7 LU to 84 LU. With these parameters, the following density ratios and surface tensions have been obtained for the pseudo-potential multiphase model.

With the parameters reported in Table 1 the following dimensionless breakup time (normalized by means of the capillary time) may be obtained, and reported in Fig. 2.

As can be observed from Fig. 2, the breakup time is showing an asymptotic behaviour with increasing the ligament radius. This may lead to strong limitations in terms of grid-size requirement when hydrodynamic time-dependent phenomena has to be analysed with the proposed method. Future work is the comparison of presented results with outcomes of other collision operators, e.g. the Multi Relaxation Time (MRT), which would allow to keep the density ratio constant with varying the kinematic viscosity. Moreover, both SRT and MRT asymptotic results should be compared with predictions coming from sharp interface hydrodynamics.

5 ACKNOWLEDGEMENTS

The numerical simulations were performed on Zeus HPC facility, at the University of Naples “Parthenope”; Zeus HPC has been realized through the Italian Government Grant PAC01$\_$00119 “MITO - Informazioni Multimediali per Oggetti Territoriali”, with Prof. E. Jannelli as the Scientific Responsible.

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 642069.