Next: Simulation Results Up: Main Previous: Introduction

II. BRIEF REVIEW OF THE LATTICE BOLTZMANN MODEL AND THREAD BREAKUP IN BULK

In the LB model that we employ [27, 29], the fluid within volume element is described in terms of the particle velocity distribution function at each point in space, is the number of particles per unit volume at node x, time t with velocity, e_a, where the subscript (a = 1,...,b) indicates the velocity direction and superscript i labels the fluid component. Time t evolves in discrete time steps, and the fluid particles can collide with each other as they move under applied forces. For this study, we use the D3Q19 (three dimensional lattice with b=19) lattice [29] where each e_a corresponds to the velocity of particles that stream to nearest neighbor sites (1 a 6) and next nearest sites (7 a 18) on a cubic lattice, while e₁₉ = 0 corresponds to a rest particle.  The units of e_a are the lattice constant divided by the time step.

Macroscopic quantities such as the number density, nⁱ(x, t), and the fluid velocity, uⁱ, of each fluid component, i, are obtained by taking suitable moment sums of .  Note that while the velocity distribution function is defined only over a discrete set of velocities, the actual macroscopic velocity field of the fluid is continuous. The time evolution of the particle velocity distribution function satisfies the following LB equation:

where the collision operator  describes the rate of change of the particle distribution due to collisions. This quantity is simplified by use of the single relaxation time approximation

where   is the equilibrium distribution at (x, t) and τ_i controls the rate of approach to equilibrium. The equilibrium distribution can be represented in the following forms for particles of each type:

with t_a = 1/18 for 1 a 6 and t_a = 1/36 for 7 a 18. The number density nⁱ and equilibrium velocity u of the fluid mixture are then defined by the weighted averages

where the sum over a ranges from 1 to 19. Similarly, d_o can be related to the temperature T by the following average:

This leads to the relation T = (1 − d_o)/2 (units are chosen such that Boltzmann's constant equals 1). The continuum limit [27] of these LB equations leads to a velocity field that is a solution of the Navier-Stokes equation where v is the kinematic viscosity,

with for i=i′. G is a "coupling constant" that controls the strength of thermodynamic interaction between the fluids. (G is the analog of the usual "exchange or van der Waals interaction" in the usual lattice models of fluid mixtures.) It has been shown that this interaction leads to phase separation and associated critical properties (phase boundaries, correlation length, surface tension) of this model fluid mixture have recently been reported [28]. Phase separation occurs in the model upon cooling for a fixed value of G or for a critical value G_c at fixed T [28]. Comparison to measurement is facilitated by expressing G in terms of a reduced variable τ_G that ranges between 0 and 1 in the two-phase region where the fluids are immiscible, a larger value of τ_G implying a larger quench depth and interfacial tension and a smaller interfacial width between the coexisting phases [28]. In real fluids, τ_G should be taken as directly comparable to the reduced temperature T-T _c /T_c where T_c is the critical temperature for phase separation. The LB model is notably a mean field model and does not account for fluctuation effects that can renormalize critical properties near T_c [28]. In the simulations below, we take τ_G=0.329 which for a T_c near room temperature (300 K) corresponds to a quench depth of 99 K, a moderately deep quench. All distances are reported in lattice spacings or as dimensionless ratios of scales. At this quench depth, the relative volume fractions in the two coexisting phases for our model are determined to equal 0.998 and 0.002 [28].

A fluid-surface interaction is incorporated by modifying Eq. (6) in the region surrounding the fluid [28, 29]. While n^i′ (x + e_a,t) normally corresponds to a particle number density, it is assigned a value 1 in Eq. (6) when x + e_a resides in the solid where the value of  , is then set to allow the solid to attract (wet) or repulse (dewet) the fluid component i [29]. A no-slip boundary condition is maintained at the wall by utilizing a second order 'bounce back' boundary condition. Here, Eq. (1a) is replaced by ,  where a' is defined such that e_a' = −e_a.

The simulations were initialized in the following way. Denoting the two fluid components as A and B, a thread of radius R with composition c_A = 0.98 and c_B = 0.02 was surrounded with a fluid having composition c_A = 0.02 and c_B =0.98. We then allowed the system to equilibrate. As a result, the thread radius shifted slightly. It should be appreciated that the interface between the thread liquid and the surrounding fluid is diffuse because of the moderate quench depth into the two-phase region (see Ref. [28]), creating some uncertainty about the interface position. To make its location specific, we defined the thread "boundary" as the location where the local fluid volume fraction is equal to 50%. This criterion normally corresponds to the inflection point of the composition profile governing the liquid-liquid interface [28]. We then defined the ratio of the tube radius R_tube to the thread radius after equilibration R_thread as a dimensionless measure of confinement, Λ ≡ R_tube / R_thread.

In the discussion below, we express time in terms of the rate of growth of the capillary instability in bulk, as described by the theory of Tomotika [18]. According to this theory, the amplitude α(t) of a perturbation of the thread boundary (infinitely long thread) grows exponentially at "early" times (see Fig. 1):

The growth rate q depends on the interfacial tension σ, the shear viscosity η_m of the "matrix" fluid, and the viscosity of the thread η_thr. The "dimensionless capillary wave growth rate factor" Φ(λ, p) is a function of the matrix-thread viscosity ratio (p = η_thr / η_m) and wavelength λ of the perturbation and this function is tabulated by Tomotika [18] (see below). The growth rate q is maximal for the wavelength, λ_max = 2πR_thread / δ, in bulk fluids (i.e., threads not subjected to confinement) where the "dimensionless wave number" δ also depends on p [18, 30]. This scale grows to predominate at long times and is thus "selected" in the late stages of thread breakup. Simulation times are expressed in the reduced time, t_red = q_∞t, where q_∞(R_tube = ∞) is the bulk capillary instability growth at the sclae, λ_max. Confinement alters the rate of growth of the capillary instability q [or equivalently Φ(λ, p)] and the wavelength λ_max of the disintegrating thread and theoretical and LB simulation results relating to these finite size effects are discussed below.

Our treatment of finite size effects on the rate of thread breakup is restricted to the case in which the fluid and matrix viscosities are equal (p=1). In this case, the Tomitika theory predicts that δ=0.56 and Φ(λ_max,p=1)=0.0714 [18, 19]. We note for comparison that δ=0.69 for the case where the viscosity of both the fluid thread and matrix are neglected ("inviscid" fluids) and that δ approaches 0 for a viscous fluid in an inviscid medium (i.e., the instability wavelength relative to the initial thread radius diverges) [16].

Our restriction to simulations of viscosity-matched fluids is made because of the large parameter space that must be investigated. Moreover, previous work on thread breakup has indicated that this assumption leads to results that are "typical" for the case where the viscosity mismatch is not large (p ≈ 1) [24, 25] and evidence for this insensitivity is described in Sec. III A.

Although no qualitative changes in the nature of the thread breakup should occur for modest values of viscosity mismatch, the growth rate q is affected by p and its variation is often of practical interest. Some insight into the magnitude of this effect can be obtained from tabulated values of Φ(λ_max,p) in the bulk limit [18]. While no concise analytic expression exists for Φ(λ_max,p), even in bulk fluids, we note that Φ(λ_max,p) is reasonably well described by the Padé approximant [31]

with a=0.963, b=456.5, c=806.2, and d=12 199. This relation agrees with the exact results of Tomotika [18] to within a maximum deviation of 0.03 for p in the broad range between 10⁻⁵ to 10⁵ where Φ(λ_max,p) varies monotonically over a corresponding range between 1 and 0. From this expression we see that thread breakup generally occurs more slowly when the thread is less viscous than the matrix fluid. Simulation of thread breakup with a "large" viscosity mismatch ( p 0.1 or p 5 requires an alternate formulation of the present LB method [32]. We defer this more general investigation to the future.

Next: Simulation Results Up: Main Previous: Introduction