Lattice Boltzmann model for combustion and detonation

Abstract

In this paper we present a lattice Boltzmann model for combustion and detonation. In this model the fluid behavior is described by a finite-difference lattice Boltzmann model by Gan et al.[Physica A, 2008, 387: 1721]. The chemical reaction is described by the Lee–Tarver model [Phys. Fluids, 1980, 23: 2362]. The reaction heat is naturally coupled with the flow behavior. Due to the separation of time scales in the chemical and thermodynamic processes, a key technique for a successful simulation is to use the operator-splitting scheme. The new model is verified and validated by well-known benchmark tests. As a specific application of the new model, we studied the simple steady detonation phenomenon. To show the merit of LB model over the traditional ones, we focus on the reaction zone to study the non-equilibrium effects. It is interesting to find that, at the von Neumann peak, the system is nearly in its thermodynamic equilibrium. At the two sides of the von Neumann peak, the system deviates from its equilibrium in opposite directions. In the front of von Neumann peak, due to the strong compression from the reaction product behind the von Neumann peak, the system experiences a sudden deviation from thermodynamic equilibrium. Behind the von Neumann peak, the release of chemical energy results in thermal expansion of the matter within the reaction zone, which drives the system to deviate the thermodynamic equilibrium in the opposite direction. From the deviation from thermodynamic equilibrium, {∆}_m^{*}, defined in this paper, one can understand more on the macroscopic effects of the system due to the deviation from its thermodynamic equilibrium.

Graphical abstract

Keywords

lattice Boltzmann method / Lee–Tarver model / viscous detonation / deviation from equilibrium

1 Introduction

In recent two decades the Lattice Boltzmann (LB) method is becoming a powerful tool for simulating complex systems [1–8]. Unlike traditional computational fluid dynamics methods, the fundamental idea of LB method is to construct simplified kinetic models that incorporate the essential physics of microscopic or mesoscopic processes. In the continuum limit the LB results should obey the macroscopic equation such as Euler equation and Navier–Stokes equation. Because of its intrinsic kinetic nature, the LB model can be used to investigate many physical phenomenon at microscopic or mesoscopic level which are generally difficult for traditional methods. It has been successfully used to study various complex fluids such as magnetohydrodynamics [9, 10], flows of suspensions [11, 12], flows through porous media [13–16], compressible fluid dynamics [17–25], multiphase flows [26–34], etc. It should be pointed out that most of these studies are focused on isothermal and incompressible systems. Given the importance of shock wave and detonation in science and engineering, LB simulation of high speed compressible fluids were attempted from the early days of LB research. But most of the current LB models for compressible flows are still subject to the constraint of low Mach number. Recently, high speed compressible LB model is constructed from the following few aspects [35]. (i) Consider appropriate additional viscosity for the situation where the system is far from equilibrium, for example, for systems with shocks [20–22, 31, 36–40]. (ii) Construct Multiple-Relaxation-Time (MRT) LB model [41–45]. (iii) Establish entropic LB model [46–49]. (iv) Introduce the FIX-UP [46, 50] and/or flux-limiter schemes [29, 51, 52]. In this work our LB model for high speed compressible flows belongs to the first category.

Detonation [53] is the supersonic combustion that usually propagates through shock wave. It is different from deflagration which is the subsonic combustion usually propagating via thermal conductivity. The detonation phenomena are highly concerned issues in science and engineering. They are widely used in the acceleration of various projectiles, mining technologies, depositing of coating to a surface or cleaning of equipment, etc. However, unintentional detonation phenomenon, when deflagration is desired, is a problem in some devices. The detonation phenomenon was first recognized by Berthelot et al. [54, 55] and Mallard et al. [56] during experiments of flame propagation. The Chapmann–Jouguet (CJ) theory [57, 58] at alternation of the 20th century considers that the detonation front can be treated as a strong discontinuity plane with chemical reaction and that the chemical reaction immediately completes within the infinite thin region. The Zeldovich–von Neumann–Doering (ZND) model [59–61] presented in 1940s treats detonation front as a shock wave in which no chemical reaction occurs. In this model the chemical reaction is triggered by shock wave and proceeds at a finite rate to completion thereafter.

Numerical simulations for combustion and detonation have been significantly improved in recent 30 years along with the development in both computational methods and available computer facilities [62–64]. There are still two main challenges in simulating such processes. The first is how to describe wave front. Traditional methods include the Eulerian cell scheme [65–68], Arbitrary–Lagrangian–Eulerian(ALE) algorithm [69, 70], level set method [71, 72], Volume of Fluid(VOF) method [73], front tracking method [74–76], etc. Though the traditional methods can treat contact discontinuity surface, they are not capable of describing the detonation wave faithfully. The second is how to simulate the chemical reaction and energy releasing processes which are coupled with the flow behavior. Several chemical reaction models [77], such as the Arrhenius kinetics, Cochran’s rate function, forest fire burn, two-step model, Lee–Tarver model, etc, have been introduced. During the energy releasing procedure, the state of system is generally far from equilibrium. However, the traditional description based on computational fluid dynamics considers only the cases where the systems can be described by the continuum Euler or Navier–Stokes equations. As a special discretization of the basic equation of non-equilibrium statistical physics, the Boltzmann equation, the LB model is intrinsically suitable for handling the non-equilibrium effects of the problems mentioned above.

To the best of our knowledge, The earliest work of LB modeling combustion phenomena was finished by Succi et al.[78]. They simulated the methane-air laminar flame under the assumptions of fast chemistry and cold flames with weak heat release. In 2000, Filippova et al. [79, 80] proposed a hybrid scheme where the flow field is solved by LB method, while the transport equation, energy equation and components equation are solved by the finite difference method. As an incomplete LB scheme, the hybrid scheme does not apply fully the advantage of LB method. Later, Yamamoto et al. [81] employed a pure LB scheme for combustion simulation. The scheme uses double distribution functions. One is employed to describe the density and velocity fields, the other describes temperature field. Lee et al.[82] proposed a quasi-incompressible model for the laminar jet diffusion flame. All those LB models mentioned above are for nearly incompressible systems which are not appropriate for detonation phenomena. In addition to incompressibility, they assumed that the chemical reaction does not affect the flow fields, which is a lethal flaw for simulating detonation and most combustion problems.

In this paper, a novel model for combustion and detonation is proposed. It couples the Finite-Difference (FD) LB model by Gan et al. [31] for fluid with the Lee–Tarver [83] model for chemical reaction. The Lee–Tarver reaction rate function is composed of two terms, the first is generally referred to as the hot spot formation term and the second is generally referred to as the reaction growth term. The former is used to investigate various spot formation processes and their subsequent growth. The latter describes the growth of the reaction. The Lee–Tarver reaction model is one of the most physically justifiable models which have produced satisfying simulations of the experimental results. The rest part of the paper is organized as below. In Section 2 we introduce the FDLB model for combustion and detonation. In Section 3 we validate the new model by simulating various detonation problems. In Section 4, we define explicitly some quantities to measure the deviation from thermodynamic equilibrium from various aspects and use the new LB model to study the non-equilibrium effects in the steady detonation procedure. Conclusions and discussions are given in Section 5.

2 FDLB model for combustion and detonation

In this section we construct a hybrid LB model for simulating the combustion and detonation phenomena. Such a model should be capable of describing the high speed compressible flow influenced by the chemical reaction. For the first step, it is assumed that there are only two species present, the reactant and the product, and the reactant is converted to the product by a one step irreversible chemical reaction governed chemical kinetics equation. The flow field is described by an FDLB model.

2.1 FDLB model for flow behavior

To model the flow behavior in the combustion and detonation process, we employ the lattice Bhatnagar–Gross–Krook (BGK) model improved by Gan et al. [31] for the density, momentum and energy. The model is composed of three components: the discrete velocity model, modified Lax–Wendroff finite difference scheme and additional viscosity term.

The discrete velocity model has the following 33 discrete velocities:

\begin{array}{}& \{\begin{array}{cc}{\nu }_0=0\hfill & (1\mathrm{a})\hfill \\ {\nu }_{ki}=v_k[\cos \left({\displaystyle \frac{\mathrm{i}\mathrm{\pi }}{4}}\right),\sin \left({\displaystyle \frac{\mathrm{i}\mathrm{\pi }}{4}}\right)],i=1,2,\mathrm{\cdots },8\hfill & (1\mathrm{b})\hfill \end{array}\end{array}

where subscript k and i indicate the k-th group, the i-th direction of the particle velocities, respectively. For the compressible fluid, the evolution of distribution function f_{ki} with BGK approximation reads

\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial f_{ki}}{\partial t}}+{\nu }_{ki}\cdot {\displaystyle \frac{\partial f_{ki}}{\partial r}}=-{\displaystyle \frac{1}{\tau }}(f_{ki}-f_{ki}^{eq})\end{array}

where r and \tau are spatial coordinate and the relaxation time. The equilibrium distribution function f_{ki}^{eq} is calculated by

\begin{array}{}\text{(3)}& \begin{array}{cc}{f}_{ki}^{eq}=\hfill & \rho F_k[(1-{\displaystyle \frac{u^2}{2T}}+{\displaystyle \frac{u^4}{8T^2}})+{\displaystyle \frac{v_{kiϵ}{u}_{ϵ}}{T}}(1-{\displaystyle \frac{u^2}{2T}})\hfill \\ \hfill & +{\displaystyle \frac{v_{kiϵ}{v}_{ki\pi }{u}_{ϵ}{u}_{\pi }}{2T^2}}\left(1-{\displaystyle \frac{u^2}{2T}}\right)\hfill \\ \hfill & +{\displaystyle \frac{v_{kiϵ}{v}_{ki\pi }{v}_{ki\vartheta }{u}_{ϵ}{u}_{\pi }{u}_{\vartheta }}{6T^3}}\hfill \\ \hfill & +{\displaystyle \frac{v_{kiϵ}{v}_{ki\pi }{v}_{ki\vartheta }{v}_{ki\xi }{u}_{ϵ}{u}_{\pi }{u}_{\vartheta }{u}_{\xi }}{24T^4}}]\hfill \end{array}\end{array}

where \rho and T are density and temperature respectively. Weight factors F_k and F_0 are calculated by

\begin{array}{}\text{(3)}& \begin{array}{cc}{F}_k=\hfill & {\displaystyle \frac{1}{v_k^2(v_k^2-v_{k+1}^2)(v_k^2-v_{k+2}^2)(v_k^2-v_{k+3}^2)}}\hfill \\ \hfill & \cdot [48T^4-6(v_{k+1}^2+v_{k+2}^2+v_{k+3}^2)T^3\hfill \\ \hfill & +(v_{k+1}^2{v}_{k+2}^2+v_{k+2}^2{v}_{k+3}^2+v_{k+3}^2{v}_{k+1}^2)T^2\hfill \\ \hfill & -{\displaystyle \frac{v_{k+1}^2{v}_{k+2}^2{v}_{k+3}^2}{4}}T],k=1,2,3,4\hfill \end{array}\end{array}

\begin{array}{}\text{(3)}& F_0=1-8(F_1+F_2+F_3+F_4)\end{array}

where

\begin{array}{}\text{(3)}& \{k+l\}=\left\{\begin{array}{ll}k+l,\hfill & \mathrm{i}\mathrm{f}k+l\le 4\hfill \\ k+l-4,\hfill & \text{if }k+l>4\hfill \end{array}\right.\end{array}

We choose v_0=0 and four nonzero v_k(k=1,2,3,4).

The following moments of the equilibrium distribution function are necessary for the current FDLB model:

\begin{array}{}\text{(4)}& {\displaystyle \sum _{ki}}{f}_{ki}^{eq}=\rho \end{array}

\begin{array}{}\text{(5)}& {\displaystyle \sum _{ki}}{v}_{ki\alpha }{f}_{ki}^{eq}=\rho u_{\alpha }\end{array}

\begin{array}{}\text{(6)}& \begin{array}{c}{\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}{v}_k^2{f}_{ki}^{eq}=e_{therm}+{\displaystyle \frac{1}{2}}\rho u^2=\rho T+{\displaystyle \frac{1}{2}}\rho u^2=P+{\displaystyle \frac{1}{2}}\rho u^2\hfill \end{array}\end{array}

\begin{array}{}\text{(7)}& {\displaystyle \sum _{ki}}{v}_{ki\alpha }{v}_{ki\beta }{f}_{ki}^{eq}=e_{therm}{\delta }_{\alpha \beta }+\rho u_{\alpha }{u}_{\beta }\end{array}

\begin{array}{}\text{(8)}& \begin{array}{c}{\displaystyle \sum _{ki}}{v}_{ki\alpha }{v}_{ki\beta }{v}_{ki\gamma }{f}_{ki}^{eq}=e_{therm}(u_{\gamma }{\delta }_{\alpha \beta }+u_{\alpha }{\delta }_{\beta \gamma }+u_{\beta }{\delta }_{\gamma \alpha })+\rho u_{\alpha }{u}_{\beta }{u}_{\gamma }\hfill \end{array}\end{array}

\begin{array}{}\text{(9)}& {\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}{v}_k^2{v}_{ki\alpha }{f}_{ki}^{eq}=2e_{therm}{u}_{\alpha }+{\displaystyle \frac{1}{2}}\rho u^2{u}_{\alpha }\end{array}

\begin{array}{}\text{(10)}& \begin{array}{c}{\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}{v}_k^2{v}_{ki\alpha }{v}_{ki\beta }{f}_{ki}^{eq}=\left(2T+{\displaystyle \frac{1}{2}}{u}^2\right)e_{therm}{\delta }_{\alpha \beta }+\left(3e_{therm}+{\displaystyle \frac{1}{2}}\rho u^2\right)u_{\alpha }{u}_{\beta }\hfill \end{array}\end{array}

where u, e_{therm}, P are the hydrodynamic velocity, internal energy and pressure, respectively. We can obtain \rho, u, T and P from Eqs. (4)–(6).

Adding the dispersion term and additional viscosity to the right hand side of Eq. (2), then performing the Lax–Wendroff finite difference scheme and second-order space-centered difference to the left hand side and the right hand side of Eq. (2), respectively, we eventually obtain the following LB equation:

\begin{array}{}\text{(11)}& \begin{array}{c}{f}_{kiI}^{new}=f_{kiI}-{\displaystyle \frac{c_{ki\alpha }}{2}}(f_{kiI+1}-f_{kiI-1})-{\displaystyle \frac{\mathrm{\Delta }t}{\tau }}(f_{kiI}-f_{kiI}^{eq})\hfill \\ +{\displaystyle \frac{c_{ki\alpha }^2}{2}}(f_{kiI+1}-2f_{kiI}+f_{kiI-1})\hfill \\ +{\displaystyle \frac{c_{ki\alpha }(1-c_{ki\alpha }^2)}{12}}(f_{kiI+2}-2f_{kiI+1}+2f_{kiI-1}-f_{kiI-2})\hfill \\ +{\displaystyle \frac{{\theta }_{\alpha I}|{\kappa }_{\alpha }|(1-|{\kappa }_{\alpha }|)}{2}}(f_{kiI+1}-2f_{kiI}+f_{kiI-1})\hfill \end{array}\end{array}

where

\begin{array}{}\text{(11)}& \begin{array}{cc}{c}_{ki\alpha }=v_{ki\alpha }\mathrm{\Delta }t/\mathrm{\Delta }{r}_{\alpha },\hfill & {\kappa }_{\alpha }=u_{\alpha }\mathrm{\Delta }t/\mathrm{\Delta }{r}_{\alpha }\hfill \end{array}\end{array}

\begin{array}{}\text{(11)}& {\theta }_{\alpha I}=\eta \left|{\displaystyle \frac{P_{\alpha I+1}-2P_{\alpha I}+P_{\alpha I-1}}{P_{\alpha I+1}+2P_{\alpha I}+P_{\alpha I-1}}}\right|\end{array}

The third suffixes, I-2, I-1, I, I+1, I+2, index the meshs node in x- or y-direction. As a monotonicity preserving scheme, the Lax–Wendroff scheme improves the numerical stability. The fourth line of Eq. (11) is the additional viscosity which makes the scheme suitable for both the low and high speed fluid flows.

The FDLB scheme gives consistent results with the Naiver–Stokes equations in the continuum limit. The traditional numerical simulation for combustion and detonation is to solve the governing equations which couple the hydrodynamic equations with the chemical kinetic equation. In most traditional studies, the effects of transportation due to viscosity and heat conduction are not explicitly taken into account, so what is used are the simple Euler equations. Currently, the combustion and detonation behaviors with viscosity and heat conduction still need to be carefully investigated. The LB model contains more fundamental information than the traditional description, especially that related to non-equilibrium effects. This characteristic will be further discussed and illustrated below.

2.2 Lee–Tarver model for combustion

Selecting appropriate chemical reaction kinetics is an important step for describing the combustion and detonation phenomena under consideration. The reaction model decides the initiation, the development and the propagation of the combustion/detonation in the simulation. The chemical reaction process is very complex. It includes a variety of reaction mechanisms. So far, most of the chemical reaction kinetics are phenomenological models. The Lee–Tarver model [83] has produced many satisfying simulations of the experimental results. It is widely used in combustion and detonation studies. The Lee–Tarver chemical reaction rate law reads

\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{d}\lambda }{\mathrm{d}t}}=a{(1-\lambda )}^x{\eta }^r+b{(1-\lambda )}^x{\lambda }^y{P}^z\end{array}

\begin{array}{}\text{(13)}& \eta ={\displaystyle \frac{V_0}{V}}-1\end{array}

Here \lambda is a parameter for the chemical reaction process or the fraction of reactant that has reacted. \eta is the relative compression. V_0 and V are specific volume ahead of the shock front and behind it, respectively. a, b, x, y, z, and r are constant parameters, where b, x, and y describe the dependence on the variables like burning area, etc. The constant b is also dependent on the pressure of incident shock wave. z describes the dependence on the local pressure. This model is composed of two terms: The first term in Eq. (12) is generally referred to as the hot spot formation term, which mainly describes the formation and subsequent growth of the hot spots. The second term in Eq. (12) is generally referred to as the growth term. It mainly describes the growth of the reaction.

As the first step, in this work we consider only the simplest case where x=y=1, z=r=0. Considering the thermal initiation, the Lee–Tarver model and initiation condition are written in the following form

\begin{array}{}\text{(14)}& \frac{\text{d}\lambda }{\text{d}t}=\left\{\begin{array}{l}a(1-\lambda )+b(1-\lambda )\lambda ,T\ge T_{th}\text{ and 0 }\le \lambda \le 1\hfill \\ 0,\text{ else}\hfill \end{array}\right.\end{array}

where T_{th} is the temperature threshold for chemical reaction.

Since the chemical reaction is much faster than the process of fluid flow, the time scale of chemical reaction and macroscopic behavior differ by several orders, the term of chemical reaction can not be treated using the same time step as the terms describing fluid flow. In this paper, we apply the operator-splitting scheme to Eq. (14). In numerical simulations the evolution can be described by two steps:

Step 1, calculate the convection contribution

\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial \lambda }{\partial t}}+u\nabla \lambda =0\end{array}

Step 2, calculate the contribution of chemical reaction

\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial \lambda }{\partial t}}=a(1-\lambda )+b(1-\lambda )\lambda \end{array}

Equation (15) can be solved by the upwind scheme

the suffixes I-1, I, I+1 index the mesh node of in x- or y-direction. Equation (16) can be solved analytically. It gives

\begin{array}{}\text{(18)}& {\lambda }_I^{n+1}={\displaystyle \frac{{\mathrm{e}}^{(a+b)\mathrm{\Delta }t}+\frac{a({\lambda }_I^n-1)}{a+b{\lambda }_I^n}}{{\mathrm{e}}^{(a+b)\mathrm{\Delta }t}+\frac{b(1-{\lambda }_I^n)}{a+b{\lambda }_I^n}}}\end{array}

2.3 Coupling of the reaction and flow behavior

In previous LB studies [78–82], it was assumed that the chemical reaction does not affect the flow fields, which is a lethal flaw for simulating detonation and most of the combustion problems. In a real combustion process, the heat of reaction is coupled with the internal energy. The increment of internal energy results in the variation of pressure which naturally influences the flow behavior. Accordingly, in our LB simulation the total internal energy increasing rate \dot{e} contains two parts: the chemical reaction part {\dot{e}}_{chem} and the original thermodynamic part {\dot{e}}_{therm}. It can be written as below:

\begin{array}{}\text{(19)}& \dot{e}={\dot{e}}_{therm}+{\dot{e}}_{chem}\end{array}

\begin{array}{}\text{(20)}& {\dot{e}}_{chem}=\dot{\lambda }\rho Q\end{array}

where Q is the reaction heat per unit mass of reactant.

In our simulations, we first obtain e_{therm} from Eq. (6). Then calculate the total internal energy e via Eqs. (19) and (20). Next, update the pressure P and temperature T by replacing e_{therm} with the total internal energy e. Finally, calculate equilibrium distribution function f_{ki}^{eq} by using the updated temperature T. In this way, the chemical reaction couples naturally with the flow behavior.

3 Verification and validation

In this section, we study several typical one- and two-dimensional detonation problems using the present model: the piston problem including effects of viscosity and heat conduction which is a typical test in studies on shock initiation of explosive; collision between detonation and shock waves; regular and Mach reflection of plane detonation wave; Richtmyer–Meshkov instability caused by detonation wave.

3.1 Piston problem with viscosity and heat conduction

Consider a detonation wave in rigid tube, followed by a piston controlled by external forces so that its velocity may be specified. Three cases are considered: (i) u_p>u_{cj}, where u_p is the velocity of piston, u_{cj} is velocity of CJ point. This case has a very simple solution: a uniform state between the front and the piston, since the flow is subsonic, a rarefaction generated at the piston will overtake the front and eventually produce the uniform steady solution corresponding to the new piston velocity; (ii) u_p=u_{cj}. We still obtain a uniform state, the CJ state. As the flow is sonic, the rarefaction generated at the piston can not overtake the detonation wave front; (iii) . The front can move no slower than velocity of detonation D_{cj}, so we still have the CJ state at the front, but the rarefaction wave reduces the velocity behind the front gradually to that of the piston. Therefore, the solution consists of a rarefaction wave followed by a uniform state.

In order to demonstrate the validity of the new model, we numerically simulate the case (ii) of piston problem, see Fig. 1. The initial macroscopic quantities are set as below:

\begin{array}{}\text{(20)}& \{\begin{array}{ll}{(\rho ,u,\nu ,T,\lambda )}_L=(1\mathrm{.}34531,0.809765,0,2\mathrm{.}65175,1)\hfill & (21\text{a})\hfill \\ {(\rho ,u,\nu ,T,\lambda )}_R=(1,0,0,1,0)\hfill & (21\text{b})\hfill \end{array}\end{array}

The periodic boundary condition is adopted on the top and bottom. The left boundary is wall, and at the right boundary we impose the right states of the initial traveling wave solution. Other parameters are \mathrm{\Delta }x=\mathrm{\Delta }y=0.001, \mathrm{\Delta }t={10}^{-5}, \tau =8\times {10}^{-6}, a=1.0, b={10}^3, Q=1.0, \gamma =2.0, T_{th}=1.1 and N_x\times N_y=1000\times 3, where ix and iy are the indexes of lattice node in the x- and y-directions, N_x and N_y are the numbers of lattice node in the x- and y-direction.

Figure 1 shows the pressure profiles of case (ii) (u_p=u_{cj}) of piston problem including effect of viscosity and heat conduction at times 0, 0.05, 0.1, 0.15, 0.2, 0.25, and 0.3. The essential features of detonation, such as the von Neumann peak, the chemical reaction zone, etc can be observed obviously.

The Hugoniot relations of detonation wave read

\begin{array}{}\text{(22)}& {\rho }_0(D-u_0)={\rho }_1(D-u_1)\end{array}

\begin{array}{}\text{(23)}& P_1-P_0={\rho }_0(D-u_0)(u_1-u_0)\end{array}

\begin{array}{}\text{(24)}& e_1-e_0={\displaystyle \frac{1}{2}}(P_1+P_0)(1/{\rho }_0-1/{\rho }_1)+\lambda Q\end{array}

where the suffixes 0 and 1 index physical quantities before and after the detonation wave, respectively, D is the velocity of detonation wave. Removing \lambda Q, Eqs. (22)–(24) reduce to the Hugoniot relations of shock wave. From Fig. 1, we obtain the velocity of detonation wave D=3.1. Other physical quantities after detonation wave are (\rho ,u,P)=(1.34587,0.81226,3.54092) from the current model. To numerically validate the current model, we apply the above LB results to Eqs. (22)–(24). The values of the Left Hand Side (LHS) are 3.1, 2.54092, 1.63095 and values of the Right Hand Side (RHS) are 3.079, 2.51801, 1.58348. Compared with the values of LHS, the relative differences of the two sides are about 0.7\%, 0.9\%, and 3\%, respectively. Considering that the viscosity and heat conductivity are not taken into account in the Hugoniot relations (22)–(24), the LB results are satisfying. Simulation results fully indicate that the LB model works for simulating the shock initiation of explosives.

3.2 Collision between detonation and shock waves

The setup is as follows: a detonation (on the left side) and a shock wave (on the right side) move face to face in the opposite directions. In order to investigate the collision between detonation and shock waves, we assume the shock strength can not ignite the combustible gas, that is, the temperature after the shock is lower than temperature threshold of explosion. The physical and chemical processes after collision can be divided into several stages: First, the collision between detonation and shock is instantaneous. Then, on the left hand side of the contact discontinuity resulting from the difference of entropies, the transmitted shock interacts with combustion production. On the right hand side of the contact discontinuity, the transmitted detonation interacts with the combustible gas compressed by shock and consequently chemical reaction occurs. The zone between transmitted shock and transmitted detonation is the domain of influence. To simulate the interaction of detonation wave and shock wave, the initial strength of the shock wave must be weaker than that of the detonation wave.

Density distributions on the x–t plane for the procedure of detonation/shock collision with different shock strengths are shown in Fig. 2(a)–(d). The initial macroscopic quantities on the left boundary are {\rho }_L=1.41995, u_L=1.52086, v_L=0, T_L=6.11682, {\lambda }_L=1.0. The parameters are Q=3.5, T_{th}=3.8. Other parameters are chosen as the same as Fig. 1. Thus, the Mach number M_D of the detonation wave is about 3.5. The Mach numbers for shock waves on the right side are M_s=1.2, 1.5, 2.0, and 2.5, respectively. Incident shock, incident detonation, transmitted shock, transmitted detonation and contact surface are obtained successfully.

The physical quantities, such as density \rho, velocity u, pressure P, and temperature T before and after collision for detonation and shock with M_s=1.5 are given in Fig. 3(a)–(d). We can find that the density and temperature change after collision while the pressure and velocity remain constants. Before collision, the gas velocity after shock (the direction is to the left) is smaller than the gas velocity after detonation wave (the direction is to the right). After collision, the direction of gas velocity is to the right. Under the current parameters the density on the right of the contact surface is greater than that on the left.

From Fig. 2(b) and Fig. 3, we obtain that the velocities of incident detonation wave and incident shock are 5.0 and -2.1, respectively. The velocities of transmitted detonation and transmitted shock waves are 4.5 and -2.7 under current parameters. Now, we numerically validate our model via two sets of Hugoniot relations. One set is for the transmitted detonation wave. The other set is for the transmitted shock wave. We first apply values of the physical quantities before and after the transmitted detonation wave to Eqs. (22)–(24). The values of the LHS are 8.41906, 11.6014, and 4\mathrm{.}89218. The values of the RHS are 8\mathrm{.}48053, 11.6567, and 4\mathrm{.}92229. Compared with values of the LHS, the relative differences of the two sides are 0.7\%, 0.5\%, and 0.6\%, respectively. Likewise, we apply values of the physical quantities before and after the transmitted shock wave to the Hugoniot relations of shock wave, i.e., Eqs. (22)–(24) with \lambda Q=0. The values of the LHS are -5\mathrm{.}99507,5.58746,1\mathrm{.}75327 and values of the RHS are -5\mathrm{.}97939, 5\mathrm{.}53874, 1\mathrm{.}75295. Compared with values of the LHS, the relative differences of the two sides are 0.3\%, 0.9\%, and 0.02\%, respectively. It is clear that our numerical results agree well with Hugoniot relations of detonation/shock waves.

3.3 Regular reflection and Mach reflection

We consider the regular and Mach reflection of detonation wave. Figure 4 illustrates a rectangular computational domain where combustible mixture is uniformly congested, two points A and B are measured 2L. Ignition arises at A and B, then two symmetric detonation waves labeled by “a” generate. Two waves collide at the symmetry plane and obtain the waves labeled by “b”. After collision, detonation waves continue to propagate, at the same time, the regular reflection is observed. As the interaction time goes on, the regular reflection changes gradually to the Mach reflection. The initial state is set as below:

\begin{array}{}\text{(25)}& (\rho ,u,\nu ,T,\lambda ){|}_{x,y,0}=\left\{\begin{array}{l}(1.35826,0.816497,0,2.59709,1),\text{ at points }A\text{ and }B\hfill \\ (1,0,0,1,0),\text{ else}\hfill \end{array}\begin{array}{cc}& \end{array}\right.\end{array}

The resolution is set to \mathrm{\Delta }x=\mathrm{\Delta }y=0\mathrm{.}001, the time step is set to \mathrm{\Delta }t={10}^{-5}, the relaxation time is \tau ={10}^{-5}, the reaction heat is Q=1.0 and the lattice size N_x\times N_y=400\times 400. The parameters of chemical reaction model are a=1.0, b={10}^3, respectively.

We show the transition process from the regular reflection to Mach reflection in Fig. 5. For simplicity, only one half iso-pressure contour lines at different times are shown. Figure 5(a) is for the process of ignition arising and detonation waves generating. Figure 5(b) shows the regular reflection. When the incident angle {\alpha }_0 exceeds the critical value {\alpha }_c [84], the regular reflection is replaced by the Mach reflection, see Fig. 5(c). The critical value satisfies following equation:

\begin{array}{}\text{(26)}& \begin{array}{c}{M}^2{(1-{\mu }^2)}^2(t_0-t_2)+M\{{(1-{\mu }^2)}^2-{(t_0-t_2)}^2-{({\mu }^2+t_0{t}_2)}^2\}-(t_0-t_2)=0\hfill \end{array}\end{array}

where M is the Mach number, t_2=\tan {\alpha }_c, {\mu }^2=\frac{\gamma -1}{\gamma +1}, t_2 is the unique real root of Eq. (26). The analytical value of critical angle in this case is {\alpha }_c={27.668}^{\circ }. The critical angle from LB result is {\alpha }_c={26.4323}^{\circ }. Compared with the analytical solution, the relative difference is about 4\%. Figure 5(d) is for the Mach reflection of detonation wave. Triple point A appears near the interaction point of detonation wave. In Fig. 5(d) the front AB shows the incident detonation wave, the front AC shows the Mach rod, and AD denotes the reflected shock wave.

3.4 Richtmyer–Meshkov instability by detonation wave

Two main types of Richtmyer–Meshkov (RM) instability problems caused by detonation wave are discussed. The first occurs when a detonation wave travels from a light medium to a heavy one. The second occurs when the detonation wave travels from a heavy medium to a light one.

Case (I) Detonation wave travels from light to heavy media.An incident detonation is a strong shock wave propagating into a reactant from the left, followed by a thin zone of reaction which supports the shock. The reactant is heated by the shock via compression, so that the ignition arises, then detonation hits an interface with sinusoidal perturbation. The initial macroscopic quantities are set as below:

\begin{array}{}\text{(26)}& \left\{\begin{array}{ll}{(\rho ,u,\nu ,T,\lambda )}_L=(1,25581,0.34570,0,1.26346,0)\hfill & (27\text{a})\hfill \\ {(\rho ,u,\nu ,T,\lambda )}_M=(1,0,0,1,0)\hfill & (27\text{b})\hfill \\ {(\rho ,u,\nu ,T,\lambda )}_R=(5.04,0,0,0.198413,0)\hfill & (27\text{c})\hfill \end{array}\right.\end{array}

The resolution is set to \mathrm{\Delta }x=\mathrm{\Delta }y=0.001, the time step is \mathrm{\Delta }t={10}^{-5}, the relaxation time is \tau ={10}^{-5} and lattice size is N_x\times N_y=600\times 100, \gamma =2.0, a=1.0, and b={10}^3 in the whole computational domain. The initial sinusoidal perturbation at the interface is x=0.15\times N_x\times \mathrm{\Delta }x+0.008\times \cos (40\pi y). At the open ends, i.e., at the left and right boundaries, we impose the in flow and outflow boundary conditions. The upper and lower boundaries are periodic. Figure 6 shows the snapshots of density field, (a)–(d) correspond to t=0, 0.06, 0.2, and 1.0, respectively. Mushroom structures for detonation wave travel from heavy medium to light one and vesicular structures for detonation travel from light medium to heavy one are successfully obtained, which is similar to the case with shock wave [37–39].

Case (II) Detonation wave travels from heavy to light media. Consider a detonation with Mach number 2.5 impinges on a sinusoidal perturbation x=0.15\times N_x\times \mathrm{\Delta }x+0.008\times \cos (40\pi y). The initial condition are set to

\begin{array}{}\text{(26)}& \left\{\begin{array}{ll}{(\rho ,u,\nu ,T,\lambda )}_L=(2.27273,1.9799,0,3.52,0)\hfill & (28\text{a})\hfill \\ {(\rho ,u,\nu ,T,\lambda )}_M=(1,0,0,1,0)\hfill & (28\text{b})\hfill \\ {(\rho ,u,\nu ,T,\lambda )}_R=(0.33333,0,0,3,0)\hfill & (28\text{c})\hfill \end{array}\begin{array}{rr}\hfill & \hfill \begin{array}{ccc}& & \end{array}\end{array}\right.\end{array}

In the simulation, we adopt the same boundary condition as in case (I). The computational domain is 0.6\times 0.1, the common parameters are \mathrm{\Delta }x=\mathrm{\Delta }y=0.001, \mathrm{\Delta }t={10}^{-5}, \tau ={10}^{-5}, \gamma =2.0, a=1.0, and b={10}^3. Figure 7 shows snapshots of density field, (a)–(d) correspond to the time t=0, 0.02, 0.05, and 0.1 respectively. In this case, interface reversal is successfully observed.

4 LB study on detonation phenomena

4.1 Phase diagram of viscous detonation

Presently, most detonation numerical studies reported are based on solving governing equations which couples Euler equation and chemical kinetics equation. This method neglects the effect of transport processes. The LB scheme, which naturally includes the effects of viscosity and heat conduction, describes more faithfully real physical procedure than traditional method.

Consider the one-dimensional steady flow of detonation [53]. Wood has shown the \lambda–V/V_0 phase plane for detonation [85], a sketch is shown in Fig. 8(a). The vertical lines \lambda =0 and \lambda =1 indicate Hugoniot curves for shock and detonation, respectively. The solid line is Rayleigh curve, which intersects the vertical line \lambda =1 at points S and W, corresponding to the strong detonation and weak detonation, respectively. The dashed lines are integral curves for detonation depending on different values of detonation velocity D. In this plane, the state point in the ZND model jumps discontinuously from initial point O to von Neumann point N, then moves up the Rayleigh curve to point S. The analogous solution of the viscous problems is an integral curve leaving point O or O^{\prime } and terminating at point S, without passing through von Neumann point N.

The \lambda–V/V_0 phase diagram for the viscous detonation obtained by our LB model is showed in Fig. 8(b), see the solid circles, where the line is for the integral curve solved by Mathematica 8.0. The parameters used here are the same as those in Fig. 1. The deviation of the integral curve from LB result originates from the different values of viscosity used in the integration of procedure by Mathematica 8.0. For Mathematica 8.0, the viscosity is P_{cj}\tau, while it equals P\tau and is variable in our LB simulating. Figure 9(a) and (b) show the influences of viscosity and heat conductivity. As the relaxation time \tau increases, the curve for V/V_0 versus \lambda become flatter, and the von Neumann peaks becomes less evident, the viscosity smoothes wave front of pressure.

4.2 Thermodynamic equilibrium versus chemical reaction

Among the seven moment relations (4)–(10) required by the present LB model, only for the first three, i.e., the definitions of density, momentum and energy [see Eqs. (4)–(6)], the equilibrium distribution function f_{ki}^{eq}can be replaced by the distribution function f_{ki}. If we replace f_{ki}^{eq} by f_{ki} in the left hand side of any one of Eqs. (7)–(10), the value of left hand side will have a deviation from that of the right hand side. This deviation may work as a measure for the deviation of system from its thermodynamic equilibrium [35]. We introduce

\begin{array}{}\text{(29)}& {∆}_m=I_m(f_{ki})-I_m(f_{ki}^{eq}), m=4,5,6,7\end{array}

where the subscript m indicate the m-th moment, I_m(f_{ki}) and I_m(f_{ki}^{eq}) are definitions of moments calculate by f_{ki} and f_{ki}^{eq}, respectively. Specifically,

\begin{array}{}\text{(30)}& \begin{array}{cc}{\mathrm{\Delta }}_{4,\alpha \beta }\hfill & =I_{4,\alpha \beta }(f_{ki})-I_{4,\alpha \beta }(f_{ki}^{eq})\hfill \\ \hfill & ={\displaystyle \sum _{ki}}{v}_{ki\alpha }{v}_{ki\beta }{f}_{ki}-{\displaystyle \sum _{ki}}{v}_{ki\alpha }{v}_{ki\beta }{f}_{ki}^{eq}\hfill \end{array}\end{array}

where {\mathrm{\Delta }}_{4,\alpha \beta }, I_{4,\alpha \beta }(f_{ki}) and I_{4,\alpha \beta }(f_{ki}^{eq}) are the components xx, xy, yx, and yy, respectively.

\begin{array}{}\text{(31)}& \begin{array}{cc}{\mathrm{\Delta }}_{5,\alpha \beta \gamma }\hfill & =I_{5,\alpha \beta \gamma }(f_{ki})-I_{5,\alpha \beta \gamma }(f_{ki}^{eq})\hfill \\ \hfill & ={\displaystyle \sum _{ki}}{v}_{ki\alpha }{v}_{ki\beta }{v}_{ki\gamma }{f}_{ki}-{\displaystyle \sum _{ki}}{v}_{ki\alpha }{v}_{ki\beta }{v}_{ki\gamma }{f}_{ki}^{eq}\hfill \end{array}\end{array}

where {\mathrm{\Delta }}_{5,\alpha \beta \gamma }, I_{5,\alpha \beta \gamma }(f_{ki}) and I_{5,\alpha \beta \gamma }(f_{ki}^{eq}) are components xxx, xxy, xyx, xyy, yxx, yxy, yyx, and yyy, respectively.

\begin{array}{}\text{(32)}& \begin{array}{cc}{\mathrm{\Delta }}_{6,\beta }=\hfill & I_{6,\beta }(f_{ki})-I_{6,\beta }(f_{ki}^{eq})\hfill \\ \hfill & ={\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}{v}_k^2{v}_{ki\alpha }{f}_{ki}-{\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}{v}_k^2{v}_{ki\alpha }{f}_{ki}^{eq}\hfill \end{array}\end{array}

where {\mathrm{\Delta }}_{6,\beta }, I_{6,\beta }(f_{ki}) and I_{6,\beta }(f_{ki}^{eq}) are the x- and y-components, respectively.

\begin{array}{}\text{(33)}& \begin{array}{cc}{\mathrm{\Delta }}_{7,\alpha \beta }\hfill & =I_{7,\alpha \beta }(f_{ki})-I_{7,\alpha \beta }(f_{ki}^{eq})\hfill \\ \hfill & ={\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}{v}_k^2{v}_{ki\alpha }{v}_{ki\beta }{f}_{ki}-{\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}{v}_k^2{v}_{ki\alpha }{v}_{ki\beta }{f}_{ki}^{eq}\hfill \end{array}\end{array}

where {\mathrm{\Delta }}_{7,\alpha \beta }, I_{7,\alpha \beta }(f_{ki}) and I_{7,\alpha \beta }(f_{ki}^{eq}) are the components xx, xy, yx, and yy, respectively. In summary, among the four moments I_4 to I_7, I_4 and I_7 are second-order tensors, I_5 is a third-order tensor and I_6 is a first-order tensor or vector.

It is clear that {∆}_m and I_m contains the information of the macroscopic flow velocity u. Furthermore, we replace {\nu }_{ki} by {\nu }_{ki}-u in any one of Eqs. (30)–(33), named {∆}_m^{*} and I_m^{*}. {∆}_m^{*} and I_m^{*} are only the manifestation of the thermo-fluctuations of molecules relative to the macroscopic flow velocity u. Specifically,

\begin{array}{}\text{(34)}& \begin{array}{cc}{\mathrm{\Delta }}_{4,\alpha \beta }^{*}\hfill & =I_{4,\alpha \beta }^{*}(f_{ki})-I_{4,\alpha \beta }^{*}(f_{ki}^{eq})\hfill \\ \hfill & ={\displaystyle \sum _{ki}}(f_{ki}-f_{ki}^{eq})(v_{ki\alpha }-u_{\alpha })(v_{ki\beta }-u_{\beta })\hfill \end{array}\end{array}

\begin{array}{}\text{(35)}& \begin{array}{cc}{\mathrm{\Delta }}_{5,\alpha \beta \gamma }^{*}=\hfill & I_{5,\alpha \beta \gamma }^{*}(f_{ki})-I_{5,\alpha \beta \gamma }^{*}(f_{ki}^{eq})={\displaystyle \sum _{ki}}(f_{ki}-f_{ki}^{eq})\hfill \\ \hfill & \cdot (v_{ki\alpha }-u_{\alpha })(v_{ki\beta }-u_{\beta })(v_{ki\gamma }-u_{\gamma })\hfill \end{array}\end{array}

\begin{array}{}\text{(36)}& \begin{array}{cc}{\mathrm{\Delta }}_{6,\beta }^{*}=\hfill & I_{6,\beta }^{*}(f_{ki})-I_{6,\beta }^{*}(f_{ki}^{eq})={\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}(f_{ki}-f_{ki}^{eq})\hfill \\ \hfill & \cdot {(v_{ki\alpha }-u_{\alpha })}^2(v_{ki\beta }-u_{\beta })\hfill \end{array}\end{array}

\begin{array}{}\text{(37)}& \begin{array}{cc}{\mathrm{\Delta }}_{7,\alpha \beta }^{*}=\hfill & I_{7,\alpha \beta }^{*}(f_{ki})-I_{7,\alpha \beta }^{*}(f_{ki}^{eq})={\displaystyle \sum _{ki}}{\displaystyle \frac{1}{2}}(f_{ki}-f_{ki}^{eq})\hfill \\ \hfill & \cdot {(v_{ki\alpha }-u_{\alpha })}^2(v_{ki\beta }-u_{\beta })(v_{ki\gamma }-u_{\gamma })\hfill \end{array}\end{array}

where, {∆}_m^{*} and I_m^{*} have the same property components as {∆}_m and I_m respectively.

Now, we use the newly introduced concepts and theory to study a simple case of detonation. For the case of CJ detonation shown in Fig. 1, Figs. 10–13 show the fourth moments and its deviations to the seventh moments and its deviations respectively. Figures 10(a)–13(a) show the profiles of physical quantities: the density \rho, pressure P, temperature T, x-component of velocity u and the fraction of product \lambda. Figures 10(b)–13(b) are for I_4 to I_7 respectively. Figures 10(c)–13(c) are for {∆}_4 to {∆}_7 respectively. Figures 10(d)–13(d) are for I_4^{*} to I_7^{*}respectively. Figures 10(e)–13(e) are for {∆}_4^{*} to {∆}_7^{*}respectively. In Figs. 10(b)–13(b) and 10(d)–13(d), the symbols are for moments calculated from f_{ki} and the solid lines are for moments calculated from f_{ki}^{eq}. The specific correspondences between the components and the symbols/lines are referred to the legends. To show clearly the system state from various sides, in Figs. 10(a)–13(a), the profiles of \rho, P, T, u and \lambda are repeatedly shown two times in the two columns to guide the eyes. The vertical dashed line in each plot indicates the von Neumann peak. Here the time t=0.3. Figure 14 shows the profiles of {∆}_4^{*} to {∆}_7^{*} with relaxation factor \tau ={10}^{-4} at time t=0.3.

From Fig. 10(a) we can find that, due to the compression effects of the shock wave, the density, pressure, temperature and flow velocity increase sharply and reach the von Neumann peak. When the temperature arrives at the threshold value for explosion, chemical reaction occurs. The fraction of product \lambda increases from 0 to 1 within the reaction zone. Since the propagation of detonation wave is faster than the pure shock wave, compared with the case of pure shock wave action, the matter behind the front of detonation wave expand quickly within the reaction time. Therefore, the density, pressure and flow velocity decrease quickly within the reaction zone. As for the variation of temperature, the situation is a little more complex. The volume expansion results in a decrease of temperature, while the chemical reaction releases heat to increase the temperature. Finally, the variation of temperature is dependent on the competition of the two mechanisms. More specifically, it is dependent on the equation of state and the reaction model. It should be pointed out that in most of the traditional studies which are based on the Euler equations, the width of the shock front is assumed to be zero. So, the density, pressure, temperature and flow velocity reach their von Neumann peak values instantaneously. Recent studies have shown that this assumption is not strictly true. Our simulation results clearly show the increasing processes of these quantities. For the case shown in Fig. 1 or Fig. 10(a), we clearly observe that the chemical reaction occurs before the temperature arrives at its von Neumann peak. This result can also be understood from Fig. 8. If the threshold temperature increases, the starting point of the chemical reaction will be closer to the von Neumann peak. If the threshold temperature is higher than the value at the von Neumann peak, no chemical reaction occurs. What we observe will be a pure shocking process. The density, pressure, temperature and flow velocity will keep constant after the front of the shock wave.

From Fig. 10 it is clear that the moments calculated from f_{ki}behavior qualitatively are the same as those from f_{ki}^{eq}. The I_4 and I_4^{*} have components showing a peak at the same position as the von Neumann peak and have components being nearly zero. The xx and yy components of I_4^{*} have the same amplitude. The {∆}_4 and {∆}_4^{*}, have components deviating significantly from zero during the reaction procedure. It is very interesting to find that, at the von Neumann peak, the system is much closer to its thermodynamic equilibrium. In the front of and behind the von Neumann peak, xx and yy components of system deviates from its equilibrium in opposite directions with the same deviation amplitude. It should also be pointed out that the deviation from thermodynamic equilibrium starts from the beginning of the shocking procedure, instead of the beginning of the chemical reaction. With the finish of the chemical reaction, the system goes back to its thermodynamic equilibrium gradually. The features of I_7, I_7^{*}, {∆}_7 and {∆}_7^{*} shown in Fig. 13 are similar. Figure 11 shows that the I_5 and I_5^{*} have components showing a peak at the same position near the von Neumann peak and have components being nearly zero. At the von Neumann peak, the system is near its thermodynamic equilibrium. Since {∆}_6and {∆}_6^{*} shown in Fig. 12 has only two components, the amplitude of its y component keeps zero, only the x component shows the deviation from thermodynamic equilibrium. Figure 14 shows that, when \tau increases to {10}^{-4}, at the von Neumann peak, the system is closer to its thermodynamic equilibrium, while in front of the peak, the amplitude of deviating from equilibrium becomes much larger.

All the non-equilibrium effects in Figs. 10 to 13 can be consistently interpreted as below. Among the four physical fields for the density, momentum, pressure and temperature, the temperature gradient is the most fundamental driving force triggering the non-equilibrium effects. The gradient of any other triggers the non-equilibrium effects via triggering macroscopic transportation which leads to temperature gradient. The temperature gradient first initiates the variance of the internal energy in the same degree of the freedom as that of the temperature gradient. Then, part of internal energy variance is transferred to other degrees of freedoms via collisions of molecules. Then, the internal energy in this degree of freedom further varies according to the temperature gradient, and so on. Only when the temperature gradient vanishes, the system can arrive at its thermodynamic equilibrium.

5 Conclusions and discussion

The combustion and detonation phenomena are widely used in the acceleration of various projectiles, mining technologies, depositing of coating to a surface or cleaning of equipment, etc. However, unintentional detonation phenomenon, when deflagration is desired, is a problem in some devices. The combustion and detonation phenomena have become highly concerned issues in science and engineering. In this work, a LB model for combustion and detonation has been presented. This model is composed of descriptions of two processes, the fluid flow and the chemical reaction. The former is described by the FDLB model [31], which gives the same results as the Navier–Stokes equations in the hydrodynamic limit. The chemical reaction is described by the Lee–Tarver model [83]. In this scheme, the heat of reaction is coupled with the internal energy. The increment of internal energy results in the variation of pressure which naturally influences the flow behavior. Since the time scale in the chemical reaction process is much smaller than that in the thermodynamical process, an operator-splitting scheme is necessary to obtain a successful simulation. In order to indicate the validity of the new model, several problems are studied: (i) the piston problem with viscosity and heat conduction which is a typical test in studies on shock initiation of explosive, (ii) collision between detonation and shock waves, (iii) regular and Mach reflection of plane detonation wave, (iv) the Richtmyer–Meshkov instability, (v) the phase diagram for viscous detonation. Simulation results fully indicate that the present model works for fundamental processes and problems of shock initiation of explosives, can capture the essential features of combustion and detonation, such as the von Neumann peak, chemical reaction zone, and etc.

In contrast with Yamatomo’s model [81], the scheme describes the density, momentum and energy using only one distribution function. Compared with the model by Filippova et al. [79, 80] and Lee et al. [82], the present model is completely thermal and compressible. Compared with all the previous LB models for combustion [78–82], the present model realizes nature coupling of the chemical reaction and flow behavior. It works not only for the combustion systems with low-Mach number but also for the detonation phenomena with high-Mach number. It can be further used to investigate the fundamental physics in various phenomena, particularly those non-equilibrium effects, related to combustion and detonation. For example, problems on the triple point trajectory, the Mach reflection and regular reflection of detonation waves, problems on shock-to-detonation and deflagration-to-detonation transitions, and problems on the Rayleigh–Taylor, Richtmyer–Meshkov and Kelvin–Helmholtz instabilities induced by detonation wave, etc. The same ideas, including that for the operator-splitting, can be used to construct appropriate LB models and simulate various combustion and detonation phenomena via choosing corresponding compressible LB and reaction models.

As a specific application of the new model, we have studied the simple steady detonation phenomenon. To show the merit of LB model over the traditional ones, we focus on the reaction zone to study the non-equilibrium effects. It is found that, at the von Neumann peak, the system is near its thermodynamic equilibrium. In the front of and behind the von Neumann peak, the system deviates from its equilibrium in opposite directions. The deviation from thermodynamic equilibrium starts from the beginning of the shocking procedure, instead of the beginning of the chemical reaction. With the finish of the chemical reaction, the system goes back to its thermodynamic equilibrium gradually. Even though most of the combustion and detonation phenomena show three-dimensional effects, if we observe from a region being small enough so that the wave front can be roughly regarded as a plane, the one-dimensional conclusions still work. From the deviation from the thermodynamic equilibrium, {∆}_m^{*}, defined in this paper, we can understand more on the macroscopic effects of the system deviating from its thermodynamic equilibrium. If we assume the normal direction of the wave front is n and t is an arbitrarily chosen tangential direction, then we have the following observations. The nn and the tt components of {∆}_4^{*} behave qualitatively in the opposite directions with the same amplitudes. The nn (tt) component shows a positive (negative) peak and a negative (positive) peak in the front of and behind the von Neumann peak, respectively. The other two components of {∆}_4^{*} keep zero during the whole procedure. The features of {∆}_7^{*} are similar. As for {∆}_6^{*} which has only two components, the amplitude of its t component keeps zero, only the n component shows the deviation from thermodynamic equilibrium. Among the eight components of {∆}_5^{*}, the nnn component has the largest amplitude, then the ones for ttn, tnt and ntt. The others keep zero. With the same method, the complex three-dimensional effects can be carefully studied and they are out of the scope of the present paper. With increasing the viscosity, the amplitude of deviation from thermodynamic equilibrium becomes larger at the two sides of the von Neumann peak, while at the von Neumann peak the system is closer to its thermodynamic equilibrium.

References

[1] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, New York: Oxford University Press, 2001

[2] R. Benzi, S. Succi, and M. Vergassola, Phys. Rep., 1992, 222(3): 145

[3] A. Xu, G. Gonnella, and A. Lamura, Phys. Rev. E, 2003, 67(5): 056105

[4] A. Xu, G. Gonnella, and A. Lamura, Phys. Rev. E, 2006, 74(1): 011505

[5] A. Xu, G. Gonnella, and A. Lamura, Physica A, 2004, 331: 10

[6] A. Xu, G. Gonnella, and A. Lamura, Physica A, 2004, 334: 750

[7] A. Xu, G. Gonnella, and A. Lamura, Physica A, 2006, 362: 42

[8] C. Aidun and J. Clausen, Annu. Rev. Fluid Mech., 2010, 42(1): 439

[9] S. Chen, H. Chen, D. Martnez, and W. Matthaeus, Phys. Rev. Lett., 1991, 67(27): 3776

[10] G. Vahala, B. Keating, M. Soe, J. Yepezand, and L. Vahala, Commun. Comput. Phys., 2008, 4: 624

[11] A. J. C. Ladd, J. Fluid Mech., 1994, 271: 285

[12] A. J. C. Ladd, J. Fluid Mech., 1994, 271: 311

[13] S. Succi, E. Foti, and F. Higuera, Europhys. Lett., 1989, 10(5): 433

[14] Y. Xu, H. Li, S. Guo, and G. Huang, Commun. Theor. Phys., 2004, 41: 949

[15] Y. Xu, Y. Liu, and G. Huang, Chin. Phys. Lett., 2004, 21: 2454

[16] Y. Xu, Y. Liu, X. Yang, and F. Wu, Commun. Theor. Phys., 2008, 49: 1319

[17] M. Watari and M. Tsutahara, Phys. Rev. E, 2003, 67(3): 036306

[18] A. Xu, Europhys. Lett., 2005, 69(2): 214

[19] A. Xu, Phys. Rev. E, 2005, 71(6): 066706

[20] Y. Gan, A. Xu, G. Zhang, and Y. Li, Commun. Theor. Phys., 2008, 50(2): 201

[21] Y. Gan, A. Xu, G. Zhang, and Y. Li, Commun. Theor. Phys., 2011, 56(3): 490

[22] Y. Gan, A. Xu, G. Zhang, and Y. Li, Phys. Rev. E, 2011, 83(5): 056704

[23] Q. Li, Y. L. He, Y. Wang, and W. Q. Tao, Phys. Rev. E, 2007, 76(5): 056705

[24] Q. Li, Y. L. He, Y. Wang, and G. H. Tang, Phys. Lett. A, 2009, 373(25): 2101

[25] Y. Wang, Y. L. He, T. Zhao, G. H. Tang, and W. Q. Tao, Int. J. Mod. Phys. C, 2007, 18(12): 1961

[26] M. R. Swift, W. R. Osborn, and J. M. Yeomans, Phys. Rev. Lett., 1995, 75(5): 830

[27] X. He, S. Chen, and R. Zhang, J. Comput. Phys., 1999, 152(2): 642

[28] M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama, and F. Toschi, Phys. Rev. E, 2007, 75(2): 026702

[29] V. Sofonea, A. Lamura, G. Gonnella, and A. Cristea, Phys. Rev. E, 2004, 70(4): 046702

[30] A. Cristea, G. Gonnella, A. Lamura, and V. Sofonea, Commun. Comput. Phys., 2010, 7: 350

[31] Y. Gan, A. Xu, G. Zhang, and Y. Li, Physica A, 2008, 387(8–9): 1721

[32] Y. Gan, A. Xu, G. Zhang, Y. Li, and H. Li, Phys. Rev. E, 2011, 84(4): 046715

[33] Y. Gan, A. Xu, G. Zhang, P. Zhang, and Y. Li, Europhys. Lett., 2012, 97(4): 44002

[34] Y. Gan, A. Xu, G. Zhang, and Y. Li, Front. Phys., 2012, 7(4): 481

[35] A. Xu, G. Zhang, Y. Gan, F. Chen, and X. Yu, Front. Phys., 2012, 7(5): 582

[36] X. F. Pan, A. Xu, G. Zhang, and S. Jiang, Int. J. Mod. Phys. C, 2007, 18(11): 1747

[37] F. Chen, A. Xu, G. Zhang, Y. Gan, T. Cheng, and Y. Li, Commun. Theor. Phys., 2009, 52: 681

[38] F. Chen, A. Xu, G. Zhang, and Y. Li, Commun. Theor. Phys., 2010, 54(6): 1121

[39] F. Chen, A. Xu, G. Zhang, and Y. Li, Commun. Theor. Phys., 2011, 55(2): 325

[40] A. N. Gorban and D. Packwood, Phys. Rev. E, 2012, 86(2): 025701(R)

[41] F. J. Higuera, S. Succi, and R. Benzi, Europhys. Lett., 1989, 9(4): 345

[42] P. Lallemand and L. S. Luo, Phys. Rev. E, 2000, 61(6): 6546

[43] P. Lallemand and L. S. Luo, Phys. Rev. E, 2003, 68(3): 036706

[44] F. Chen, A. Xu, G. Zhang, and Y. Li, Phys. Lett. A, 2011, 375(21): 2129

[45] F. Chen, A. Xu, G. Zhang, Y. Li, and S. Succi, Europhys. Lett., 2010, 90(5): 54003

[46] F. Tosi, S. Ubertini, S. Succi, H. Chen, and I. V. Karlin, Math. Comput. Simul., 2006, 72(2–6): 227

[47] S. Ansumali and I. V. Karlin, J. Stat. Phys., 2002, 107(1/2): 291

[48] S. S. Chikatamarla and I. V. Karlin, Phys. Rev. Lett., 2006, 97(19): 190601

[49] S. S. Chikatamarla and I. V. Karlin, Phys. Rev. E, 2009, 79(4): 046701

[50] Y. Li, R. Shock, R. Zhang, and H. Chen, J. Fluid Mech., 2004, 519: 273

[51] Y. Gan, A. Xu, G. Zhang, and Y. Li, Commun. Theor. Phys., 2011, 56(3): 490

[52] F. Chen, A. Xu, G. Zhang, and Y. Li, Commun. Theor. Phys., 2011, 56(2): 333

[53] W. Fickett and W. C. Davis, Detonation,Theory and Experiment, New York: Dover Publications, Inc. Mineola, 1979

[54] M. Berthelot, P. Vielle, and C. R. Hebd, Sceances Acad. Sci., 1881, 93: 18

[55] M. Berthelot, P. Vielle, and C. R. Hebd, Sceances Acad. Sci., 1882, 94: 149

[56] E. Mallard, H. Le Chatelier, and C. R. Hebd, Sceances Acad. Sci., 1881, 93: 145

[57] D. L. Chapmann, Philos. Mag., 1899, 47: 90

[58] E. Jouguet, J. Math. Pures Appl., 1905, 1: 347

[59] Ya. B. Zeldovich and S. A. Kompaneets, Zh. Eksp. Teor. Fiz., 1940, 10: 542

[60] J. Von Neumann, Theory of Detonation Waves, New York: Macmillan, 1942

[61] W. Doering, Ann. Phys., 1943, 43: 421

[62] C. L. Mader, Numerical Modeling of Explosives and Propellants, New York: CRC Press, 2008

[63] C. Wang, X. Zhang, C. W. Shu, and J. Ning, J. Comput. Phys., 2012, 231(2): 653

[64] S. Tan, C. Wang, C. W. Shu, and J. Ning, J. Comput. Phys., 2012, 231(6): 2510

[65] S. Karni, J. Comput. Phys., 1994, 112(1): 31

[66] A. Marquina and P. Mulet, J. Comput. Phys., 2003, 185(1): 120

[67] J. J. Quirk and S. Karni, J. Fluid Mech., 1996, 318: 129

[68] K. M. Shyue, J. Comput. Phys., 1998, 142(1): 208

[69] R. Loubre, P. H. Maire, M. Shashkov, J. Breil, and S. Galera, J. Comput. Phys., 2010, 229: 4724

[70] S. Galera, P. H. Maire, and J. Breil, J. Comput. Phys., 2010, 229(16): 5755

[71] S. Ssher and R. P. Fedkiw, J. Comput. Phys., 2001, 169: 463

[72] M. Sussman, P. Smereka, and S. Osher, J. Comput. Phys., 1994, 114(1): 146

[73] R. Scardovelli and S. Zaleski, Annu. Rev. Fluid Mech., 1999, 31(1): 567

[74] G. Tryggvason, B. Bunner, A. Esmaeeli, and D. Juric, J. Comput. Phys., 2001, 169: 708

[75] J. Glimm, J. W. Grove, X. L. Li, and D. C. Tan, SIAM J. Sci. Comput., 2000, 21(6): 2240

[76] D. K. Mao, J. Comput. Phys., 2007, 226(2): 1550

[77] J. Sun and J. Zhu, Theory of Detonation Physics, Beijing: National Defense Industry Press, 1995 (in Chinese)

[78] S. Succi, G. Bella, and F. Papetti, J. Sci. Comput., 1997, 12(4): 395

[79] O. Filippova and D. Hanel, J. Comput. Phys., 2000, 158(2): 139

[80] O. Filippova and D. Hanel, Comput. Phys. Commun., 2000, 129(1–3): 267

[81] K. Yamamoto, X. He, and G. D. Doolen, J. Stat. Phys., 2002, 107(1/2): 367

[82] T. Lee, C. L. Lin, and L. D. Chen, J. Comput. Phys., 2006, 215(1): 133

[83] E. L. Lee and C. M. Tarver, Phys. Fluids, 1980, 23(12): 2362

[84] R. Courant and K. O. Friedrichs, Supersonic Flow and Shock Waves, New York: Interscienc Publishers Inc., 1948

[85] W. W. Wood, Phys. Fluids, 1963, 6(8): 1081

Rights & permissions

Higher Education Press and Springer-Verlag Berlin Heidelberg 2013