Mesoscopic modeling of a two-phase flow in the presence of boundaries: The contact angle

We present a mesoscopic model, based on the Boltzmann equation, for the interaction between a solid wall and a nonideal fluid. We present an analytic derivation of the contact angle in terms of the surface tension between the liquid-gas, the liquid-solid, and the gas-solid phases. We study the dependency of the contact angle on the two free parameters of the model, which determine the interaction between the fluid and the boundaries, i.e. the equivalent of the wall density and of the wall-fluid potential in molecular dynamics studies. We compare the analytical results obtained in the hydrodynamical limit for the density profile and for the surface tension expression with the numerical simulations. We compare also our two-phase approach with some exact results obtained by E. Lauga and H. Stone [J. Fluid. Mech. 489, 55 (2003)] and J. Philip [Z. Angew. Math. Phys. 23, 960 (1972)] for a pure hydrodynamical incompressible fluid based on Navier-Stokes equations with boundary conditions made up of alternating slip and no-slip strips. Finally, we show how to overcome some theoretical limitations connected with the discretized Boltzmann scheme proposed by X. Shan and H. Chen [Phys. Rev. E 49, 2941 (1994)] and we discuss the equivalence between the surface tension defined in terms of the mechanical equilibrium and in terms of the Maxwell construction.

Figure 1

The stationary state configuration for the density field as a function of the relative distance from the wall. We present the density field normalized to the center channel density for two different values of the average density of the system: $⟨\rho ⟩=1.1{\rho }_l$ (◻) and $⟨\rho ⟩=0.7{\rho }_l$ (엯). In both cases we integrate numerically the lattice Boltzmann equation with $\tau =1$ in a $2d$ channel with two walls (bottom and lower) and a periodic boundary condition in the stream-wise direction. The dimensions of the channel are $L_x\times L_y=80\times 45$ and ${\mathcal{G}}_b=-6.0$, being ${\mathcal{G}}_c=-4.0$ the critical point of the system. The liquid and gas density for this value are, respectively, ${\rho }_l=2.65$ and ${\rho }_g=0.07$ in lattice Boltzmann units. Note that the density tends to match a given value at the wall (for $y∕L_y=0$) which is different from both liquid and gas values; this is because of the chosen boundary condition, $\psi ({\rho }_w=0.5)$.

Figure 2

Surface tensions in a three phase system for ${\mathcal{G}}_b=-6.0$ as a function of ${\rho }_w$. The horizontal line represents the surface tension between liquid and gas $\left({\gamma }_{lg}\right)$. The surface tension between solid and gas $\left({\gamma }_{sg}\right)$ is zero for ${\rho }_w={\rho }_g=0.07$ and equal to ${\gamma }_{lg}$ for ${\rho }_w={\rho }_l=2.65$. The surface tension between liquid and solid $\left({\gamma }_{ls}\right)$ reflects the opposite behavior.

Figure 3

Contact angle for the case ${\mathcal{G}}_b=-6.0$ as a function of ${\rho }_w$. The analytical estimate is compared to the numerical results of lattice Boltzmann simulations (+). The lattice Boltzmann simulations are carried out in a $2d$ domain periodic along the stream-wise direction and with two walls (upper and lower wall). The grid mesh used is $L_x\times L_y=80\times 45$ and the relaxation time chosen is $\tau =1$. To produce the steady state with a drop of liquid the system has been initialized with a nonhomogeneous condition with a square spot of liquid on the lower wall.

Figure 4

(Color online) Left: typical density profile from which we extract the contact angle. The contact angle in the numeric is extracted by estimating the angle made by the surface where the density profile is equal to the average value $({\rho }_l-{\rho }_g)∕2$ being ${\rho }_l=2.65$ and ${\rho }_g=0.07$ for the present case with ${\mathcal{G}}_b=-6.0$. The computational parameters are the same used as for Fig. 3. To produce the three phases equilibrium the system has been initialized to a nonhomogeneous condition with a square spot of liquid in contact with the lower wall. Right: a vertical snapshot of the left panel for $x=45$.

Figure 5

Velocity profiles for pressure driven Stokes flow. Profiles are obtained integrating numerically the Stokes equation for unit pressure drop and kinematic viscosity (34) in a homogeneous (in the stream-wise direction) channel with a height $L=50$ grid points. The profiles are then normalized with respect to the center channel velocity of the Poiseuille flow $\left(u_c^p\right)$ for the same geometry. From top to bottom we show different values of ${\rho }_w$ corresponding to different contact angles: ${\rho }_w=0.13$ $(\theta =166°)$, ${\rho }_w=0.27$ $(\theta =145°)$, ${\rho }_w=0.49$ $(\theta =110°)$.

Figure 6

Slip length for the profiles in the Stokes approximation. In the top panel we show an extrapolation of the profiles shown in the previous figure. The slip length is obtained extrapolating to zero the parabolic profile obtained as a best fit in the bulk region of the channel. In the bottom panel we show the slip length estimated indirectly from the mass flow rate gain with respect to the Poiseuille profile for different contact angles (엯) and compare to the phenomenological result (35), solid line. The good agreement is obtained with a prefactor in (35) fixed to be 1.2.

Figure 7

Results obtained from the mechanical stability conditions (26), ◻, and thermodynamic coexistence curve (37), 엯, for liquid and gas densities. The results of numerical simulation are also reported (×). The numerical simulations are carried out in a fully periodic $2d$ setup with grid mesh $L_x\times L_y=90\times 90$ and a relaxation time $\tau =1$. The system is then initiated to have a flat interface along the $x$ direction. The top branch refers to the liquid density, ${\rho }_l$, while the bottom branch to the gas density, ${\rho }_g$.

Figure 8

Surface tension as a function of the interaction parameter ${\mathcal{G}}_b$. The surface tension calculated using both expressions in Eq. (43), the case with ${\partial }_y\rho$ (엯) and the case with ${\partial }_y\psi$ (◻). To show the importance of the normalization factor (42) we also show the surface tension calculated by simply replacing ${\partial }_y\psi \to {\partial }_y\rho$ in (28) (filled squares).

Figure 9

(Color online) Schematic configurations for the simulations of laminar flows with boundary conditions of free shear and no slip. The alternating pattern of longitudinal (top) and transverse (bottom) strips of free shear is produced by varying the boundary density from complete wetting (${\rho }_w={\rho }_l$, contact angle 0°, no slip) and perfect dewetting in a strip $H$ (${\rho }_w={\rho }_g$, contact angle 180°, free shear). With these setups the “slip” probability is $\xi =H∕L$ being $L=L_x\left(L_y\right)$ for transverse (longitudinal) strips. All the details of numerical simulations are given in Fig. 10.

Figure 10

Effective slip length as a function of the free shear percentage $\left(\xi \right)$ in a multiphase approach with alternating strips of wetting and nonwetting properties. The numerical simulations are carried out in 3D setups with $L_x\times L_y\times L_z=64\times 64\times 200$ with periodic boundary conditions in the stream-wise $\left(x\right)$ and span-wise directions $\left(y\right)$. The interaction parameter is ${\mathcal{G}}_b=-6.0$ and the relaxation time is $\tau =1$. The corresponding liquid and gas densities $({\rho }_l=2.65,{\rho }_g=0.07)$ are used to mimic no slip and free shear, respectively. Then, the steady state effective slip lengths for transverse (◻) and longitudinal (엯) strips normalized to the pattern dimension ($L_x$ for transverse and $L_y$ for longitudinal) are compared with exact analytical estimates given for stokes flow with alternating non slip and free shear (lines). Bottom panel: The velocity profile in the inlet of the channel and in the middle of a free shear strip from the integration of the incompressible lattice Boltzmann equation with mixed boundary conditions (straight lines) are compared with the results of the mesoscopic multiphase approach (◻ in the inlet and 엯 in the middle of free shear strip). The incompressible lattice Boltzmann equation has been used with the same relaxation time and an average density equal to the one of the mesoscopic approach. Both simulations have been forced with a constant pressure gradient so as to reproduce a center channel velocity equal to $u_c=0.04$ in LB units.

Figure 11

Contact angle for a two parameter model. Top panel: for a given value of ${\rho }_w$ we show the contact angle for different values of ${\mathcal{G}}_w$: ${\rho }_w=1.3$ (엯), ${\rho }_w=0.6$ (◻). The analytical estimate using the mechanical stability is compared with lattice Boltzmann simulations. The details of the numerical simulations are the same for Fig. 3 with the introduction of a wall-fluid force as in (44). Bottom panel: for the same contact angle (130°) we report different profiles obtained for different choices of the parameters $({\rho }_w,{\mathcal{G}}_w)$. In particular, the value of ${\rho }_w$ chosen is ${\rho }_w=1.3$ (엯) ${\rho }_w=0.85$ (◻) ${\rho }_w=0.3$ (x) and the corresponding value of ${\mathcal{G}}_w$ has been chosen to reproduce the same contact angle of 130° estimated from (47).

Figure 12

The schematic figure of the contact angle setup.