Numerical Example

In the previous sections several methods for the incorporation of full energy conservation in discrete Boltzmann methods were described. At this point it would be useful to illustrate these ideas with a numerical example and verify that, for example, the total energy is conserved. Here the method of rescaling the velocity and temperature (or internal energy) in a thermal equilibrium distribution in order to construct the modified BGK collision term will be utilized.

Consider the case of a low density system of particles governed by a 6-12 Lennard-Jones (LJ) interaction potential. In this low density regime we make use of the molecular chaos approximation and, for simplicity, approximate the two point correlation function, to lowest order, by $g(r)= exp[-\beta V(r)]$ where $\beta= 1/KT$ [14]. It is assumed that the density varies slowly over the range the LJ interaction. It is then easy to numerically evaluate $\phi, \vec{Z},$ and $\stackrel {\leftrightarrow}{B}$. Note $\stackrel {\leftrightarrow}{B}$ is diagonal for this case since V(r)is a central potential. Once $\phi, \vec{Z},$and $\stackrel {\leftrightarrow}{B}$ are evaluated the velocity and temperature are rescaled according to Eqs. 21 and 23. The time derivative of $\phi$ was numerically determined to first order from information stored from the previous time step. Although a crude approximation, it serves to illustrate the method. The details of the thermal lattice Boltzmann model are given in the appendix.

The fluid system is initialized with a small random perturbation in the density and then allowed to relax to a equilibrium configuration. To demonstrate energy conservation, the global change of potential energy $\Delta \Phi$and change in kinetic energy, $\Delta KE$ as a function of time is shown in Figure 1. Note the sum $\Delta \Phi + \Delta KE=0$ to the numerical accuracy of the computer. Similar results concerning conservation of mass and momentum were obtained. For this simulation it should be noted that a significant fraction of the exchange of potential energy and kinetic energy took place over 1 to 2 mean free collision times. This result is similar to that predicted in a study of LJ fluids [26].

Figure 1: The global change in potential energy and kinetic energy for a fluid system initially subject to a small density perturbation. The system size studied was 121³ in lattice units, $\Delta x$. Energy is in units of $m(\Delta x/\Delta t)^2$ where m=1 and the time step $\Delta t=1$ for simplicity in this simulation. The solid line corresponds to the change in kinetic energy and the dashed line corresponds to the potential energy change at each time step. The sum of $\Delta \Phi$ and $\Delta KE$ (dotted line) is clearly equal to zero demonstrating the algorithm conserves energy.

It is interesting to see how introducing total energy conservation effects the time decay of a density fluctuation autocorrelation function, $G(t)=<\delta \rho(t)^2>$. Figure 2 shows the time decay of G(t)for different Lennard-Jones interaction strengths. Here, results are presented for the case of a quasi-two dimensional system (1024x1024x2) to help reduce finite size effects. Clearly, the constraint of conserving both kinetic and potential energy slows the decay of G(t). This effect diminishes over long times (and for weaker interaction strength) as the system approaches thermal equilibrium. The simulations carried out were consistent with $G(t) \sim t^{-1}$ at long times. Such scaling behavior is associated with the decay of the thermal diffusivity mode in the study of hydrodynamic fluctuations [27] which predicts $G(t)\sim t^{-D/2}$ where D is the spatial dimension.

Figure 2: The density autocorrelation function, G(t), versus time, t, for a fluid system initially subject to a small density perturbation. The solid curves correspond to the fully energy conserving equation for different Lennard-Jones interaction strengths. From top to bottom the strength of interaction is 0.04, 0.03, 0.02, 0.004, and 0.0 in units of KT. The straight line has a slope of -1 to help guide the eye. The dashed line corresponds to the Shan Chen model with interaction strength equal to .04.

Figure 2 also includes G(t) for the case where the constraint of energy conservation is removed while still including a forcing due to the interaction potential. This case is similar to the Shan-Chen model. Comparing to the fully energy conserving case with the same LJ interaction strength, there is a clear difference in the early rate of decay of G(t) indicating that the short time behavior of this system strongly depends on whether or not one imposes energy conservation. This may be an important consideration when modeling systems that are dynamically driven and subject to frequent perturbations over time.