Numerical simulation of spinodal dewetting using single-component multiphase pseudopotential lattice Boltzmann method at high density ratio

Spinodal dewetting is the spontaneous rupture and dewetting of thin liquid film (thickness less than 100 nm) on a solid substrate due to the attractive intermolecular forces between the interfaces of the liquid-bounding fluid and the liquid bounded solid substrate. Most of the reported numerical studies are performed by simplification of the Navier–Stokes (NS) equations using lubrication approximation. The lubrication approximation, however, is not valid for all the types of liquid thin films. For example, in metallic films, contact angles are greater than 30o, and inertial forces are also significant. So, to understand thin film dewetting dynamics, there is a need to solve the NS equation without simplification. The present numerical study investigates the spinodal dewetting of thin liquid films under van der Waals force by indirectly solving NS equations using one of the mesoscopic approaches, the lattice Boltzmann method (LBM). The stability analysis is carried out using a single-component multiphase pseudopotential LBM with a multiple relaxation time collision operator at the density ratio of 98.48, in both 2D (two-dimension) and 3D (three-dimension). D2Q9 (D2 represents two-dimension and Q9 represents nine possible microscopic velocities in which a particle can move) and D3Q15 (D3 represents three-dimension and Q15 represents 15 possible microscopic velocities in which a particle can move) lattice models are used in 2D and 3D, respectively. In-house codes are developed using C language, and the 3D LBM codes have been parallelized using a message-passing interface. The reported dewetting stages, namely, the arrangement of fluctuations on a dominant wavelength, hole formation, and expansion, are successfully observed with the present numerical method in both the dimensions. The breakup of liquid threads between the holes into droplets due to Rayleigh instability has also been captured in 3D. In 2D, at the time of initial rupture of the film, the average dominant wavelength λdom,rupavg is calculated by DFT (discrete Fourier transform), and it was observed that λdom,rupavg was much greater than that predicted by linear stability analysis λdom,LSA, using lubrication approximation. This study reveals that this deviation can be contributed to the shear stresses present at the interface, high contact angles, and diffused interface.