Spectra of spatial fluctuations in lattice gas cellular automata fluids

Spectra of spatial fluctuations in lattice gas cellular automata fluids

PHYSiCA ELSEVIER Physica A 224 (1996) 180-187 Spectra of spatial fluctuations in lattice gas cellular automata fluids Shankar P. Das School of Physi...

366KB Sizes 2 Downloads 68 Views

PHYSiCA ELSEVIER

Physica A 224 (1996) 180-187

Spectra of spatial fluctuations in lattice gas cellular automata fluids Shankar P. Das School of Physical Sciences, Jawaharlal Nehru University, New Delhi 110067, India

Abstract Spectra or decay rates of spatial fluctuations for different Lattice Gas Cellular Automata (LGCA) Fluids are discussed using the Boltzmann equation. We demonstrate how the spectra provide information on the time and length scales of hydrodynamic behavior, isotropy and existence of spurious conservation laws in different LGCA models. The range of validity of generalized hydrodynamics with wave-number-dependent sound speed and transport coefficients is indicated. These properties give a quantitative explanation of several existing, but unexplained simulation results on the so called "negative bulk viscosities", and on the strong dispersion in the simulation results on sound damping. The use of the spectra to analyse the spatial instabilities occurring in certain models of LGCA is also indicated.

1. The eigenvalue problem The eigenvalue problem studied in this paper is based on the Boltzmann equation, where correlations between time-dependent fluctuations are neglected. In formulating the problem we follow essentially the presentation of Ref. [ 1-3] including the discussion of thermal models [4] with multispeed lattice gases. The average occupation f ( r , ci, t) of a state (r, ci) satisfies the nonlinear lattice Boltzmann equation,

f ( r + ci, ci, t + 1) = f ( r , ci, t) + l i ( f ( t ) ),

(1)

where i = ( 1 , 2 . . . . . b) labels the velocity channels of a b-bit model, and I i ( f ) represents the nonlinear collision term, accounting for gains and losses in the f ' s . The stationary solution of this equation, f ( r , cj, oo) = f~y, with I i ( f °) = 0, is the equilibrium distribution, and we restrict ourselves to basic equilibria where the total momentum vanishes and the fluid is macroscopically at rest. The equation for the linear excitation, 8 f i( r, ci, t) = f ( r, ci, t) - f ~ i , is obtained by linearizing the collision operator li around 0378-4371/96/$15.00 (~) 1996 Elsevier Science B.V. All rights reserved SSDI 0378-4371 (95)00309-6

S.P. Das/Physica A 224 (1996) 180-187

181

the stationary solution if/, i.e. l i ( f ° + 6 f ) = /2ijtfj + O ( 8 f 2 ) , where summation convention has been used. A spatial fluctuation has the general form 6 f ( r , c , t ) = Oa ( k, c) exp [ i k . r + za ( k ) t]. The fluctuation at wave number k denoted by Oa (k, c), satisfies the basic eigenvalue equation, [e zA(k)+ik'c - 1 - / 2 ] l¢~b,~(k, c ) = 0

(2)

in matrix notation. Here ~Oa( k, c ) is a b-dimensional vector with components ~ba( k, c j) ( j = 1,2 . . . . . b) and/2, x and exp(ik, c ) are b x b matrices. The last two are diagonal, i.e. Kij = Ki~ij = f~i ( 1 -- f~i )C~ij and [exp ik.c ]ij = ~ij exp( ik.cj ). The linearized collision operator can be expressed in terms of the basic transition probabilities that define the collision rules of the lattice gas [5]. If the model satisfies the detailed balance condition, then the symmetry /2ijKj = Ki/2ij (no summation) holds. The eigenmodes of (1) are defined through K e n ( k , c, t)

=

[e-ik'c( 1 + /2)it K~b,~(k, c ) = exp[z~(k)t] Ken(k, c ).

(3)

The b x b matrix g2ij can be calculated analytically or numerically from its definition in terms of the collision rules for the model [5]. In an athermal model the collision matrix /2 and its eigenvalues depend only on the density. In a thermal model they depend both on density and temperature. Once 12 is known, the eigenvalues of the matrix (2) can be calculated numerically as a function of the wave number k and the thermodynamic parameters. The eigenvalues, za(k) = Reza(k) + i l m z a ( k ) , are in general complex, and Im z a ( k ) is defined modulo 2~-. We first consider a two-dimensional lattice gas, 7-bit FHP-UI lattice gas, named after Frisch, Hasslacher and Pomeau [5,6]. The resulting dispersion relations are shown in Fig. 1, where (a) and (b) show the real and imaginary parts of all seven eigenvalues in the direction of highest symmetry, i.e. k parallel to the y-axis. The real parts of z + ( k ) coincide and are shown as the lower branch in Fig. la. The real parts of hydrodynamic and kinetic eigenvalues are well separated for k < kl, where kl is the wave number at which for the first time a hydrodynamic and kinetic eigenvalue are equal. Wave numbers on the order of kl characterize the kinetic regime. In Fig. 1, where the density p = 1.4, one has kl ~- 2.1. In the limit of small k (actually k£0 << 1 where £0 is the mean free path) the eigenvalues are given by the hydrodynamic dispersion relations, z+ ( k ) = :Ficsk - F k 2, z ± ( k ) = -~,k 2,

(4)

where Cs is the speed of sound, ~, is the shear viscosity and F the sound damping constant, given in athermal models by F = ½(~, + ~'), with ~" the bulk viscosity. Because of periodicity the same behavior pertains around every site Q of the reciprocal lattice, with Ikl in (4) replaced by Ik - Q I . Fig. lb also shows a propagation gap ( I m z ± ( k ) = +Tr) in the interval 3.2 < Ikl _ 4.1. Similar propagation gaps of sound modes have been found from the revised Enskog theory for dense hard sphere systems and have been observed in the dynamic structure factor S(k, to) in neutron scattering experiments

S.P Das/Physica A 224 (1996) 180-187

182

3.0

I

I

i

(a) 7

2.0 4 N ¢)

6

n~ '

5

1.0

0.0

0.0

2.0

4.0

6.0

kY 4.0

i

i

I

(b) 2.0

t'q

0.0

E I

-2.0

-4.0 0.0

2.0

4.0

6.0

kY Fig. 1. Real and imaginary part of the spectrum z~(k), for the 7-bit FHP-III model at density p = 1.4, with k = (0,k) in (a), (b) respectively. The labels refer to the "soft" shear (_1_) and sound (:E) modes, the staggered modes (s), and the "hard" kinetic (4, 5, 6, 7) modes.

on liquid Argon and liquid metals [7]. With the increase of k, the dispersion relations of classical hydrodynamics with k-independent transport coefficients break down and we enter the regime of generalized hydrodynamics, where the dispersion relations can be represented approximately by (4) with a slowly varying k-dependent speed of sound cs(k) and k-dependent transport coefficients F(k) and v(k). Inspection of Fig. lb shows that the speed of sound cs(k) is constant up to k "~ 1. The other important aspect for fluid dynamics simulations with lattice gases is the length scale at which the transport coefficients become k dependent. For this we consider

183

S.P Das/Physica A 224 (1996) 180-187

1.5

1

I

i

v(k) ..... F(k) l

1.0

0.5

0.0

0.1

0.2

0.3

0.4

0.5

k Fig. 2. Shear viscosityv(k) and sound dampingconstant F(k) versus k, for the 6-bit FHP-I model at p = 1.8 with k II .~. a lattice gas with six velocity states ci per node, and referred to as the athermal 6-bit FHP-I model [6]. For the FHP-I model the shear viscosity ~,(k) = - z ± ( k ) / k 2 (solid line) and the sound damping constant F ( k ) =_ - R e z ± ( k ) / k 2 (dashed line) as a function of k with k parallel to the y-direction are shown in Fig. 2. In this model as k ~ 0 the sound damping constant F ( k ) ~ F = ½9, since the bulk viscosity vanishes in a single speed model. Inspection of Fig. 2 shows that the transport coefficients are approximately constant for kh ~ 0.1. Consequently, if one wants to measure the constant transport coefficients in computer simulations from the decay of a sinusoidal wave, its wavelength A = 2zr/k should satisfy ?t > Ah ~ 60 lattice units. The analysis presented in this paper is based on the linear Boltzmann equation, which accounts for short-ranged spatial correlations (i.e. k-dependent transport coefficients) and all memory effects are neglected. Therefore the present theory does not yield any frequency-dependent transport coefficient. Next, we give a brief description of a few useful applications of the study of the spectra for the LGCA models.

Ikl _<

2. Negative bulk viscosity The generalized hydrodynamic effects on transport coefficients, shown in Fig. 2, also provide a quantitative explanation for some older puzzling simulation results in the literature [6,8], described as a "negative bulk viscosity". The simulation data for the FHP-I model [6] are shown in Fig. 3 as dots, referring to the shear viscosity, and squares, referring to the sound damping constant. The computer experiments were carried out at wavelengths between 30 and 80 lattice units, corresponding to k ~ 0.2 and k '~ 0.08 respectively, and the results were averaged over this k-interval. In Fig. 3 the simulation

S.P Das/Physica A 224 (1996) 180-187

184 1.6

i

• v(k) (simulationT] o 2r(k)(simulation) I

I

- v(k) (theory) |

~[.

- 2F(k)(theor),) ]

!.2

0.8

~\ ' U

0

\o

0.4

i

0.0

i

i

i

o

/

./

-/

i

1.0

i

2.0

i

3.0

P Fig. 3. k-dependent shear viscosity v and sound damping constant F versus density p, for the 6-bit FHP-I model, with k = 0.2. The simulation data were taken from Ref. [9]. The graphs show that the k-dependent bulk viscosity ((k) = 2F(k) - v(k) is negative for this k-direction. data are compared with the k-dependent shear viscosity u ( k ) (solid line) and the kdependent sound damping constant 2 F ( k ) (dashed line) for k = 0.2. In addition to the (positive) wave-number-dependent damping coefficients v ( k ) and F ( k ) one may also introduce a wave-number-dependent bulk viscosity, ( ( k ) = 2 F ( k ) - v ( k ) , based on the analogous formula, ( = 2 F - v, valid in the long wavelength limit, where ( > 0 in the FHP-III model and ( = 0 in the FHP-I model. The observation from Fig. 3, that 2Fsim < Psim [6,8], has been misinterpreted as a "negative bulk viscosity". However, the good agreement between theory and simulations, especially near the large-k-end of the averaging interval, convincingly shows that the experiments were done outside the range of classical hydrodynamics. It is important to note that the "negative bulk viscosity" effect is not a finite size effect. It also occurs on an infinite lattice at the same k-vector.

3. T h e r m a l effects We consider the simplest thermal lattice gas to be a 9-bit model, defined on a square lattice, with four speed 1, four speed ~ particles and one rest particle. It is discussed in Refs. [9,10]. In thermal models energy is conserved. This leads to an additional hydrodynamic mode, the heat mode, labeled a = T, with a dispersion relation for k~0,

zr( k ) = - D r k 2,

(5)

where Dr is thermal diffusivity. In Fig. 4 we show the k-dependent heat diffusivity computed from the corresponding hydrodynamic eigenvalue ZT through the relation D r ( k ) - - z r ( k ) / k 2. The corresponding k-dependent viscosity is shown in Fig. 5.

S.P. Das/Physica A 224 (1996) 180-187 1.2

.

,

.

,

.

-

185

,

k)/(~,0)

(t,i) (.).8 [-

0.4

0.0

,

0.0

i

0.5

,

i

,

1.0

,

,

2.0

1.5

k Fig. 4. Heat diffusivity Dr(k) versus k, for the 9-bit model, with k II (1,0) (solid line) and k (dashed line). 0.4

,

,

II

(1, 1)

,

kll

(I,0)

1

0.3

'-~.~ 0.2 >

0.1

0.0 0.0

~

~

0.5

1.0

'

'

1.5

2.0

k Fig. 5. Viscosity ~,(k) versus k, for the 9-bit model with k II (1,0) (solid line) and k II (1, 1) (dashed line).

Generalized hydrodynamic effects are particularly strong in Fig. 4, where D r ( k ) is constant only for A = 2¢r/k > 1~h ~'~ 50 lattice units. For ,~ < Ah the heat diffusivity D r ( k ) rapidly falls off by an order of magnitude. For the shear viscosity the effects of generalized hydrodynamics are less pronounced and become only noticeable for A < 25. This demonstrates the fact that a 9-bits model is not suitable to study Fourier's heat law with constant transport coefficients. The relation between the heat flux and the temperature gradient is strongly non-local, extending over a distance Ah ~-- 50, as compared to Ah '~ 3 in the FHP-III model. In simulations of the 9-bit square lattice gas one would need a system size a factor ( 5 0 / 3 ) 2 larger than in FHP-III to have an

186

s.P Das/Physica A 224 (1996) 180-187

equivalent range of k-values where standard hydrodynamic behavior holds.

4. Isotropy The isotropy of the spectrum of the FHP-III lattice gas model is due to the hexagonal symmetry of the triangular lattice underlying these models. The real and imaginary parts of za(k) can be computed [1] for various directions of k. For small k the eigenvalues are approximately independent of the direction of the wave vector. At larger k-values the isotropy breaks down for all lattices. With decreasing density the isotropy breaks down at smaller and smaller k-values. On square lattices the spectrum is anisotropic even on the longest wave lengths, as can be seen from the study of the spectra of eigenvalues [ 1 ]. The shear viscosity is computed from the decay rate of the transverse shear mode. The fourth rank viscosity tensor appearing in the macrodynamic equations on a square lattice is not isotropic, and contains in principle three different types of scalar viscosities, u, O and ( ( ( --- 0 in the 9-bit square lattice gas), as discussed in Ref. [ 10]. From that reference one readily derives the angular dependence of the damping constant of the shear mode as Ik[ ~ 0. The result is p ( k ) = u±~±e( k ) = u cos 2 ~b + O sin 2 ~b,

(6)

where cos(½~b) = ~:x is the x-component of the unit vector k = k/Ik[. In Fig. 5 the solid line and the dashed line correspond to ~b = 0 and ~b = 7r/4 respectively.

5. Spurious soft modes Another conspicuous feature in the spectra za (k) is the occurrence of a soft mode, at k = (0, 27r/x/3) in Fig. 2a. These soft modes are the well-known staggered momentum densities [ 11 ] appearing in lattice gases as a consequence of extra conserved quantities due to the discrete spatial symmetry in these systems. The staggered diffusivity is not isotropic, and can be decomposed into a longitudinal (DII) and a transverse ( D ± ) diffusivity. For example with k = (0, ky) (see Fig. la), one has

( z o ( k y ) = -Oil

2 ky-

v/~j

.

(7)

The staggered diffusivities have been calculated in Refs. [ 11,12] for 6- and 7-bit athermal triangular lattice gases and in Ref. [ 10] for the thermal 9-bit lattice gas.

6. Spatial instability A recent and very interesting application of the spectra deals with the dynamics of phase separation in lattice gas models for thermodynamically unstable systems (spinodal

S.P Das/Physica A 224 (1996) 180-187

187

decomposition), as occurring in the long-range Van der Waals-type model o f Appert and Zaleski [ 13,14] or in the biased lattice gas o f Ref. [15]. Calculations o f spectra show that the real part o f one or more eigenvalues z a ( k ) becomes positive for Ikl < k0, implying that the corresponding modes with a wavelength larger than 2 ~ / k o are unstable. The onset o f coarsening in spinodal decomposition and the typical wavelength at which this occurs are well described by the eigenvalue spectrum o f the linear theory. In fact the Boltzmann equation for the above lattice gas provides a microscopic model that leads to the phenomenological C a h n - H i l l i a r d equation for spinodal decomposition [ 16]. In concluding this paper we emphasize that our discussion has been restricted to the more generic features o f the eigenvalue spectra o f single particle fluctuations, using an analysis based on mean field theory. The main importance o f analyzing the spectral properties is that they provide clear and crucial information about the basic criteria that lattice gases have to satisfy in order to qualify as models for non-equilibrium fluids.

Acknowledgements I sincerely acknowledge the organizers for the excellent facilities provided at the conference on "Dynamics o f Complex Systems" in Calcutta, in August 1995. The "Stichting voor Fundamenteel Onderzoek der Materie ( F O M ) " , which is sponsored by the "Nederlandse Organistie voor Wetenschappelijk Onderzoek ( N W O ) " is acknowledged for financial support during a visit in July 1992 and in December 1994 to the University o f Utrecht.

References [1] S.E Das, H. Bussemaker and M.H. Ernst, Phys. Rev. E 48 (1993) 245. [2] R. Brito, M.H. Ernst and T.R. Kirkpatrick, in: Discrete models of Fluid Dynamics, Series on Advances in Mathematics for Applied Sciences, A.S. Alves, ed. (World Scientific, Singapore, 1991 ) p. 198. 13] L.S. Luo, H. Chen, S. Chen, G.D. Doolen and Y.C. Lee, Phys. Rev. A 45 (1991) 7097. [4] M.H. Ernst and S.E Das, J. Stat. Phys. 66 (1992) 465. [5] U. Frisch, D. d'Humi~res, B. Hasslacher E LaUemand, Y. Pomeau and J.-E Rivet, Complex Systems l (1987) 649. [6] D. d'Humi~res and E Lallemand, Complex Systems 1 (1987) 599. [7] I.M. de Schepper and E.G.D. Cohen, J. Stat. Phys. 27 (1982) 223. [8] J.P. Rivet and U. Frisch, C.R. Acad. Sc. Paris 302 (1986) 267. [9] D. d'Humi~res, E Lallemand and U. Frisch, Europhys. Lett. 2 (1986) 291. 110] S.E Das and M.H. Ernst, Physica A 187 (1992) 191. [ ll] G. Zanetti, Phys. Rev. A 40 (1989) 1539. [12] R. Brito, M.H. Ernst and T. Kirkpatrick, J. Stat. Phys. 62 (1991) 283. [ 13] C. Appert and S. Zaleski, Phys. Rev. Lett. 64 (1990) I. [ 14] M. Gerits, M.H. Ernst and D. Frenkel, Phys. Rev. E 48 (1993) 988. [ 15] H.J. Bussemaker, Ph. D Thesis, University of Utrecht, The Netherlands (1995). [ 16l K. Binder, in: Phase Transformations in Materials, E Haasen, ed., Materials Science and Technology, Vol. 5 (VCH, Weinheim).