# Two dimensional MRT LB model for compressible and incompressible flows

###### Abstract

In the paper we extend the Multiple-Relaxation-Time (MRT) Lattice Boltzmann (LB) model proposed in [Europhys. Lett. 90, 54003 (2010)] so that it is suitable also for incompressible flows. To decrease the artificial oscillations, the convection term is discretized by the flux limiter scheme with splitting technique. New model is validated by some well-known benchmark tests, including Riemann problem and Couette flow, and satisfying agreements are obtained between the simulation results and analytical ones. In order to show the merit of LB model over traditional methods, the non-equilibrium characteristics of system are solved. The simulation results are consistent with physical analysis.

###### pacs:

47.11.-j, 51.10.+y, 05.20.DdKeywords: lattice Boltzmann method; multiple-relaxation-time; flux limiter technique; Prandtl numbers effect; non-equilibrium characteristic

## I Introduction

In recent years, the Lattice Boltzmann (LB) method has emerged as an attractive computational approach for complex physical systemSS (); BSV (); XGL1 (). The lattice Bhatnagar-Gross-Krook (BGK) model, based on a single-relaxation-time approximation, is the simplest and the most popular form. However, this simplicity also leads to some deficiencies, such as the numerical stability problem, and fixed Prandtl number. To overcome these deficiencies of BGK model, the Multiple-Relaxation-Time (MRT) lattice Boltzmann methodHSB (); 13 () has been developed, and successfully used in simulating various fluid flow problemsLallemand1 (); Lallemand2 (); PRE036701 (); JCP539 (); CF957 (); CMA1392 (); JCP453 (); JCP7774 (); PRE016705 (); JPA055501 (). Most of the existing MRT models work only for isothermal system. To simulate system with temperature field, many attempts have been madepre036706 (); pre026705 (); Mezrhab2010 ().

Besides the models mentioned above, we proposed a MRT Finite Difference lattice Boltzmann model for compressible flows with arbitrary specific heat ratio and Prandtl number in previous workchenepl (). In the model, the kinetic moment space and the equilibria of nonconserved moments are constructed according to the seven-moment relations associated with the local equilibrium distribution function. Numerical experiments showed that compressible flows with strong shocks can be well simulated by this model.

In the paper we extend the MRT LB model so that it is suitable also for incompressible flows. In order to efficiently decrease the unphysical oscillations, the flux limiter schemeSofonea1 (); Sofonea3 (); chenctp () with splitting technique is incorporated into the new model. When the system deviates more from equilibrium, the LB simulation can give more physical informationxufop (); yan (); JCP (), such as the non-equilibrium characteristics of system. Here, in the new MRT LB model, the non-equilibrium characteristics of system are solved through a dynamic procedure where a shock wave propagates from a heavy medium to a light one.

The rest of the paper is organized as follows. Section II presents the extended MRT LB model. Section III describes the finite difference schemes. section IV is for the validation and verification of the new LB model. Non-equilibrium characteristics are shown and analyzed in section V. Section VI makes the conclusion for the present paper.

## Ii Model description

According to the main strategy of MRT LB method, the MRT LB equation can be described as:

(1) |

where and are the particle distribution function in the velocity space and the kinetic moment space respectively, is the discrete particle velocity, , ,, is the number of discrete velocities, the subscript indicates or . The matrix is the diagonal relaxation matrix. is the transformation matrix between the velocity space and the kinetic moment space. , is an element of the transformation matrix. is the equilibrium value of distribution function in the kinetic moment space.

In the previous work, we constructed a two-dimensional MRT LB model based on the model by Kataoka and TsutaharaKataoka () (see Fig. 1):

where cyc indicates the cyclic permutation. Transformation matrix and the equilibrium distribution function in the moment space are chosen according to the seven-moment relations (see Appendix for details). At the continuous limit, the above formulation recovers the following Navier-Stokes (NS) equations:

(2a) | |||||

(2b) | |||||

(2c) | |||||

(2d) | |||||

where , is twice of the total energy, and is a constant related to the specific-heat-ratio . | |||||

In order to maintain the isotropy constraint of viscous stress tensor and heat conductivity, some of the relaxation parameters should be equal to one another, namely , . The above NS equations reduce to |

(3a) | |||

(3b) | |||

(3c) | |||

where the viscosity , the bulk viscosity , , . | |||

However, the viscous coefficient in the energy equation (3c) is not consistent with that in the momentum equation (3b). By modifying the collision operators of the moments related to energy flux: |

(4a) | |||

(4b) | |||

we get the following energy equation: |

(5) |

where the thermal conductivity .

## Iii Finite Difference Scheme

In the original LB modelchenepl (), the time evolution is based on the usual first-order forward Euler scheme, while space discretization is performed through a Lax-Wendroff scheme. In this work, the flux limiter scheme with splitting technique corresponding to the MRT model is adopted. The proposed flux limiter scheme can efficiently decrease the unphysical oscillations around the interfaces.

Figure 2 shows the characteristic lines in the flux limiter scheme and corresponding projections in and directions. and are corresponding projections of node in the and directions. Let be the value of distribution function at time in the node along the direction , we rewrite the evolution of in node at time step as follows,

(6) |

where

(7) |

() and () are and components of the outgoing (incoming) flux in node along the direction ,

(8a) | |||

(8b) | |||

(8c) | |||

(8d) | |||

The flux limiter is expressed as |

(9) |

where the smoothness functions are

(10a) | |||

(10b) |

The Lax-Wendroff scheme is recovered for the flux limiter , and the first order upwind scheme is recovered when .

## Iv validation and verification

### iv.1 Performance on discontinuity

In order to check the performance of flux limiter scheme on discontinuity, we construct the following problem

(11) |

is the length of computational domain. In the direction, is set, where the macroscopic quantities adopt the initial values. In the direction, the periodic boundary condition is adopted. The physical quantities on the two sides satisfy the Hugoniot relations. Fig. 3 shows the simulation results of density, pressure, component of velocity, and temperature at time using different space discretization schemes. The parameters are , , , , and other collision parameters are . The simulations with Lax-Wendroff scheme have strong unphysical oscillations in the shocked region. The second order upwind scheme results in unphysical ‘overshoot’ phenomena at the shock front. The simulation results with flux limiter scheme are much more accurate, and this scheme has the ability to decrease the unphysical oscillations at the discontinuity.

### iv.2 Lax shock tube problem

The initial condition of the problem is:

(12) |

The profiles of density, pressure, component of velocity, and temperature at are shown in Fig. 4, where the exact solutions are presented with solid lines for comparison. The parameters are , , , , , and other collision parameters are . Obviously, the simulation results agree well with the exact solutions.

The above simulations show that compressible flows, especially those with discontinuity and shock waves, can be well simulated by the present model.

### iv.3 Couette flow

Here we conduct a series of numerical simulations of Couette flow. In the simulation, the left wall is fixed and the right wall moves at speed , . The initial state of the fluid is , , , . The simulation results are compared with the analytical solution:

(13) |

where and are temperatures of the left and right walls (, ), is the width of the channel. Periodic boundary conditions are applied to the bottom and top boundaries, and the left and right walls adopt the nonequilibrium extrapolation method. Fig. 5 shows the comparison of LB results with analytical solutions for thermal Couette Flows. (a) corresponds to , and (b) corresponds to . It is clearly shown that the simulation results of new model are in agreement with the analytical solutions, and the Prandtl number effects are successfully captured. New model is suitable for incompressible flows.

## V Non-equilibrium characteristic

To show the merit of LB method over traditional ones, in this section we study the non-equilibrium characteristics using the new model. Among the moment relations required by each LB model, only for the first three (density, momentum and energy), the equilibrium distribution function can be replaced by the distribution function . If we replace by in the left hand of other moment relations, the value of left side will have a difference from that of the right side. This difference represents the deviation of system from its thermodynamic equilibriumxufop (); yan (); JCP (). In this MRT LB model, the kinetic moment space and the corresponding equilibria of nonconserved moments are constructed according to the seven-moment relations. So, the deviation from equilibrium in this model can be defined as . contains the information of macroscopic flow velocity . Furthermore, we replace by in the transformation matrix M, named (see Appendix for details). is only the manifestation of molecular thermalmotion and does not contain the information of macroscopic flow.

Now, we study the following dynamic procedure. An incident shock wave with Mach number travels from a heavy medium and hits a light one, where the two different fluids are separated by an unperturbed interface. The initial macroscopic quantities are as follows:

where the subscripts , , indicate the shock wave region, the heavy medium region, and the light medium region. In our simulations, the computational domain is , and divided into mesh-cells. The initial position of shock wave is , the unperturbed interface lies at the position . Inflow boundary is applied at the left side, outflow boundary is applied at the right side, and periodic boundary conditions are applied at the top and bottom boundaries. in the whole domain. The density, pressure, component of velocity and temperature profiles and () on the center line at time are shown in Fig. 6. The parameters are , , and other collision parameters are .

In the figures, the system shows three different interfaces, rarefaction wave, material interface, and shock wave. Physical quantities change significantly at the three interfaces, and vertical lines indicate the positions of interfaces. The system starts to deviate from equilibrium once the physical quantities starts to change. When the physical quantities arrives at its steady-state required by the Hugoniot relations, the system goes back to its equilibrium state. The peak values of deviations at shock wave interface are larger than the others. This is because the shock dynamic procedure is faster than the other two processes, and the system has less time to relax to its thermodynamic equilibrium.

At the interfaces, , and have small amplitudes. contains two parts, and components of internal translational kinetic energy. This indicates that the two parts deviate from equilibrium in opposite directions with the same amplitude. shows an opposite deviation for the rarefaction wave interface and the shock interface. The physical reason is as below. The temperature gradient first initiates variance of the internal kinetic energy in the direction of temperature gradient. (Here, the temperature shows gradient in the direction.) Then, part of internal kinetic energy variance is transferred to other degrees of freedoms via collisions of molecules. The internal kinetic energy in the temperature gradient direction further varies, and so on. The shock wave increases density, pressure and temperature, while the rarefaction wave decreases those quantities. So, shows a negative deviation for the rarefaction wave interface, while shows a positive deviation for the shock interface. The values of at material interface and shock wave interface have the same order, and are much larger than that at rarefaction wave. This is because the sizes of temperature variation near the material interface and shock wave differ little, and larger than that near the rarefaction wave. When the temperature gradient vanishes, the system attains its thermodynamic equilibrium.

## Vi Conclusions

In the paper a MRT LB model for compressible flows is extended so that it is suitable also for incompressible flows. In order to efficiently decrease the unphysical oscillations, space discretization adopts flux limiter scheme with splitting technique. It is validated and verified via same well-known benchmark tests, including Riemann problem and Couette flow, and satisfying agreements are obtained between the new model results and analytical ones. In order to show the merit of LB model over traditional methods, we studied the behaviors of system deviating from its equilibrium through a dynamic procedure where shock wave propagates from a heavy material to a light one. The simulation results are consistent with the physical analysis.

## Acknowledgements

The authors would like to sincerely thank S. Succi and C. Lin for many instructive discussions. We acknowledge support of National Natural Science Foundation of China[under Grant Nos. 11075021 and 11047020]. AX and GZ acknowledge support of the Science Foundation of CAEP [under Grant Nos. 2012B0101014 and 2011A0201002] and the Foundation of State Key Laboratory of Explosion Science and Technology[under Grant No.KFJJ14-1M].

## Appendix A Transformation matrix and equilibria of the nonconserved moments

In the model by Kataoka and Tsutahara, the local equilibrium distribution function satisfies the following relations:

(14a) | |||

(14b) | |||

(14c) | |||

(14d) | |||

(14e) | |||

(14f) | |||

(14g) |

where a parameter is introduced, in order to describe the extra-degrees of freedom corresponding to molecular rotation and/or vibration, where for , ,, and for , , .

The transformation matrix in the MRT model is composed as below: ,

(15a) | |||

(15b) | |||

(15c) | |||

(15d) | |||

(15e) | |||

(15f) | |||

(15g) | |||

(15h) | |||

(15i) | |||

(15j) | |||

(15k) | |||

(15l) | |||

(15m) | |||

(15n) | |||

(15o) | |||

(15p) |

where .

Replacing by in the transformation matrix M, matrix is expressed as follows: ,

(16a) | |||

(16b) | |||

(16c) | |||

(16d) | |||

(16e) | |||

(16f) | |||

(16g) | |||

(16h) | |||

(16i) | |||

(16j) | |||

(16k) | |||

(16l) | |||

(16m) | |||

(16n) | |||

(16o) | |||

(16p) |

where .

The equilibria of nonconserved moments are as follows:

(17a) | |||

(17b) | |||

(17c) | |||

(17d) | |||

(17e) | |||

(17f) | |||

(17g) | |||

(17h) | |||

(17i) | |||

(17j) | |||

(17k) | |||

(17l) |

## References

- (1) S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, New York (2001).
- (2) R. Benzi, S. Succi, and M. Vergassola, Phys. Rep. 222 (1992) 145.
- (3) Aiguo Xu, G. Gonnella, and A. Lamura, Phys. Rev. E 67 (2003) 056105; Phys. Rev. E 74 (2006) 011505; Physica A 331 (2004) 10; Physica A 344 (2004) 750; Physica A 362 (2006) 42.
- (4) F. J. Higuera, S. Succi, and R. Benzi, Europhys. Lett. 9 (1989) 345. F. J. Higuera, J. Jimenez, Europhys. Lett. 9 (1989) 662.
- (5) D. d’Humières, Generalized lattice-Boltzmann equations, in Rarefied Gas Dynamics: Theory and Simulations, edited by B. D. Shizgal and D. P Weaver, Progress in Astronautics and Aeronautics, Vol. 159(AIAA Press, Washington, DC, 1992), pp450-458.
- (6) P. Lallemand, L. S. Luo, Phys. Rev. E 61 (2000) 6546.
- (7) D. d’Humières, M. Bouzidi, and P. Lallemand, Phys. Rev. E 63 (2001) 066702.
- (8) M. E. McCracken, J. Abraham, Phys. Rev. E 71 (2005) 036701.
- (9) K. N. Premnath, J. Abraham, J. Comp. Phys. 224 (2007) 539.
- (10) H. D. Yu, L. S. Luo, and S. S. Girimaji, Computers Fluids 35 (2006) 957.
- (11) P. Asinari, Comput. Math. Appl. 55 (2008) 1392.
- (12) I. Rasin, S. Succi, and W. Miller, J. Comput. Phys. 206 (2006) 453.
- (13) H. Yoshida, M. Nagaoka, J. Comput. Phys. 229 (2010) 7774.
- (14) Z. H. Chai, T. S. Zhao, Phys. Rev. E 86 (2012) 016705.
- (15) J. J. Huang, H. B. Huang, C. Shu, Y. T. Chew and S. L. Wang, J. Phys. A: Math. Theor. 46 (2013) 055501.
- (16) P. Lallemand, L. S. Luo, Phys. Rev. E, 68 (2003) 036706.
- (17) L. Zheng, B.C. Shi, Z.L. Guo, Phys. Rev. E, 78 (2008) 026705.
- (18) A. Mezrhab, M. A. Moussaouia, M. Jami, H. Naji, and M. Bouzidi, Phys. Lett. A, 374 (2010) 3499.
- (19) F. Chen, Aiguo Xu, G. C. Zhang, Y. J. Li, and S. Succi, Europhys. Lett. 90 (2010) 54003.
- (20) A. Cristea, V. Sofonea, Proceedings of the Romanian Academy, Series A 4 (2003) 59; Central European J. Phys. 2 (2004) 382.
- (21) V. Sofonea, A. Lamura, G. Gonnella, A. Cristea, Phys. Rev. E 70 (2004) 046702.
- (22) F. Chen, Aiguo Xu, G. C. Zhang, Y. J. Li, Commun. Theor. Phys. 56 (2011) 333.
- (23) Aiguo Xu. G. C. Zhang, Y. B. Gan, F. Chen, X. Yu, Front. Phys. 7 (2012) 582.
- (24) B. Yan, Aiguo Xu. G. C. Zhang, Y. J. Ying, H. Li, Front. Phys. 8 (2013) 94.
- (25) C. Lin, Aiguo Xu, G. Zhang, Y. Li, and S. Succi, e-print arXiv:1302.7104v1.
- (26) T. Kataoka, M. Tsutahara, Phys. Rev. E 69 (2004) 035701(R).