A stability property for a mono-dimensional three velocities scheme with relative velocity

1. Introduction

Studying stability of lattice Boltzmann schemes is a non-trivial problem. Classically for this purpose, the scheme is linearized around a constant state and a Fourier analysis is performed. We refer to the work of Lallemand and Luo [9] for the D2Q9 scheme applied to hydrodynamics. Note also the work of Ginzburg et al. [8] extending the Fourier analysis to a wide variety of different two and three dimensional lattice Boltzmann schemes. Instabilities and their interpretattion in terms of bulk viscosity has been proposed by Dellar [2]. But no mathematical analysis has been performed.

A new way of improving stability is proposed by Geier [7], who proposed a new generalized lattice Boltzmann scheme with the approach of relative velocities and utilized it for hydrodynamics applications [4]. The tentative of analysis of this method for a two-dimensional scalar linear equation has been also proposed [5].

Even if it is a difficult task, it is well known that Fourier analysis is not the best method for analyzing nonlinear hyperbolic equations. Total variation diminishing schemes, developed for suppressing oscillations in higher order CFD algorithms, provide an alternative nonlinear stability analysis tool for analysing the schemes for nonlinear wave propagation. The convergence of such schemes is well established [10]. The underlying stability notion concerns the maximum principle. This notion can be extended to nonlinear cases and a first attempt has been proposed in [1] for lattice Boltzmann schemes for the D1Q2 scheme for scalar nonlinear hyperbolic equations.

In this contribution, we propose to investigate the stability in the maximum sense, of a linear mono-dimensional lattice Boltzmann scheme with three velocities. More precisely, we look to a positivity constraint for a particle distribution function in the context of a relative velocities.

In Section 2, we describe the scheme and the underlying advection model. More precisely, the local relaxation step is written as a linear operator on the particle distribution functions. If all the coefficients of the underlying matrix are nonnegative, the non-negativity of the distribution is maintained during this step. Because the transport step is just a change of locus, the non-negativity is maintained for the whole time step of the scheme. The question is then to find appropriate conditions to handle this property. In Section 3, a necessary and sufficient condition is derived on the parameters in order to ensure that the scheme has the stability property. In Section 4, we completely describe the classical case where the relative velocity is reduced to zero. In Section 5, the general case is presented. With an analytical study for necessary conditions and a numerical one for a complete description of the stability zones. In Section 6, numerical experiments show the correlation of the positivity constraint for a particle distribution and the presence of oscillations for discontinuous profiles.

2. Description of the framework

2.1. Description of the scheme

In this contribution, we investigate a mono-dimensional three velocities linear lattice Boltzmann scheme with relative velocity [4]. Denoting $Δ x$ the spatial step, $Δ t$ the time step, and $λ = Δ x / Δ t$ the lattive velocity, this scheme can be described in a generalized D. d’Humière’s framework [3]:

[(1)]
the 3 velocities $c_{1} = - 1$ , $c_{2} = 0$ , and $c_{3} = 1$ ;
the 3 associated distributions $f_{1}$ , $f_{2}$ , and $f_{3}$ ;
the 3 moments
$ρ$ , $q (u)$ , and $ε (u)$ given by

$ρ=\smashoperator[]∑1≤j≤3fj,q(u)=λ\smashoperator[]∑1≤j≤3(cj−u)fj,ε(u)=3λ2\smashoperator[]∑1≤j≤3(cj−u)2fj−2λ2\smashoperator[]∑1≤j≤3fj,$

where $u$ is a given scalar representing the relative velocity;
the equilibrium value of the 3 moments

$ρ\eq=ρ,q\eq(u)=λ(V−u)ρ,ε\eq(u)=λ2(3u2−6uV+α)ρ,$

where $V$ and $α$ are given scalars (without loose of generality, we assume that $V > 0$ );
the 2 relaxation parameters $s$ and $s^{'}$ such that the relaxation phase reads

$q⋆(u)=(1−s)q(u)+sq\eq(u),ε⋆(u)=(1−s′)ε(u)+s′ε\eq(u).$

In this formalism, the moments are defined as polynomial functions of the discrete velocities and the discrete distribution functions. Indeed, introducting $P_{1} = 1$ , $P_{2} = λ X$ , and $P_{3} = λ^{2} (3 X^{2} - 2)$ , the three moments read

ρ = \sum 1 \leq j \leq 3 P_{1} (c_{j} - u) f_{j}, q (u) = \sum 1 \leq j \leq 3 P_{2} (c_{j} - u) f_{j}, ε (u) = \sum 1 \leq j \leq 3 P_{3} (c_{j} - u) f_{j} .

The equilibrium values are chosen such that the equilibrium distributions do not depend on the relative velocity $u$ . Indeed, we have:

f\eqj=16ρ(2+3cjV+(3c2j−2)α),1≤j≤3.

Note that this scheme can be used (see e.g. [6]) to simulate a scalar transport equation with constant velocity $λ V$ given by

\partial_{t} ρ + λ V \partial_{x} ρ = 0.

Denoting $x_{k} = k Δ x$ , $k \in Z$ , the discret spatial points and $t^{n} = n Δ t$ , $n \in N$ , the discret time points, one time step of the scheme reads

f_{j} (t^{n} + Δ t, x_{k} + c_{j} Δ t) = f_{j} (t^{n + 1}, x_{k + j}) = f_{j}^{⋆} (t^{n}, x_{k}), 1 \leq j \leq 3.

2.2. Matrix notation for the relaxation step

In this contribution, we are concerned with the non-negativity of the particle distribution functions. This property can be viewed as refering to a weak maximum principle for linear schemes. Indeed, it is always possible, by adding a constant, to assume that all the particle distribution functions are initialy nonnegative. If the scheme ensures that this property of non-negativity remains as time marches, each particle distribution function is then bounded, as their total sum is conserved. Moreover, as the transport step consists simply in exchanging the position of the particle distribution functions, we focus on the relaxation step.

We use a matrix notation for the relaxation step as it can be read as a multiplication by a matrix. As this step is local in space, we omit the dependancy on time and space. We define the vector of the distribution functions

\vectF

\vectF=(f1f2f3)T.

One relaxation step then reads

\vectF⋆=\matR\vectF,

where the matrix $\matR$ is defined by

\matR=\matiM\matiT(\matid+\matS(\matT\matE\matiT−\matid))\matT\matM,

with

	$\matM=⎛⎜⎝111−λ0λλ2−2λ2λ2⎞⎟⎠,\matT=⎛⎜⎝100−λu103λ2u2−6λu1⎞⎟⎠,$
	$\matS=⎛⎜⎝0000s000s′⎞⎟⎠,\matE=⎛⎜⎝100Vλ00αλ200⎞⎟⎠,\matid=⎛⎜⎝100010001⎞⎟⎠.$

The coefficients of the matrix $\matM$ are obtained by the relations $\matMk,j=Pk(cj)$ , $1 \leq k, j \leq 3$ and those of the matrix $\matT$ by the change of basis formula: $\matTk,l$ is the coefficient of the $l^{th}$ -element $P_{l} (c_{j})$ in the definition of $P_{k} (c_{j} - u)$ according to

Pk(cj−u)=∑1≤l≤3\matTk,lPl(cj),1≤k,j≤3.

The matrix $\matM$ is then the change of basis that transforms the vector $\vectF$ into the vector $\vectMz=(ρ,q(0),ε(0))T$ :

\vectMz=\matM\vectF,\vectM=\matT\matM\vectF.

The matrix $\matT$ can then be viewed as the change of basis matrix from the classical moments without relative velocity toward the moments with relative velocity.

2.3. Remark on the choice of the moments

Note that the last moment $ε (u)$ that is chosen in this contribution is not the energy but a moment that is orthogonal to the two first ones, $ρ$ and $q (u)$ . In this section, we show that all the results of the contribution would be identical by choosing the last moment as the energy: the relaxation matrix $\matR$ would still be the same.

We consider two schemes with two different choices of polynomials: the moments of the first scheme are defined by $(P_{1}, P_{2}, P_{3})$ while the moments of the second scheme by $({ˆ P}_{1}, {ˆ P}_{2}, {ˆ P}_{3})$ . The first moment is the same in both schemes to be able to simulate the same transport equation. We then have ${ˆ P}_{1} = P_{1}$ . We define $\matC$ the change of basis matrix associated to the tranformation $\matM$ into $\matMh$ :

\matMh=\matC\matM.

The first line of $\matC$ is then $(1, 0, 0)$ .

We assume that the equilibrium values of the distribution functions are the same, that is $\matEh=\matC\matE$ , and that the relaxation parameters are the same, that is $\matSh=\matC\matS$ . Then, we have $\matRh=\matR$ for all $(s, s^{'})$ iff

ˆP2∈\spn(P1,P2),ˆP3∈\spn(P1,P3),inR[X]/X(X−1)(X+1).

Proof.

First, we immediately obtain the following relations by identifying the coefficients of $\matMh$ and of $\matM$ :

ˆPk(cj)=∑1≤l≤3\matCk,lPl(cj),1≤j≤3.

We deduce that

\matTh=\matC\matT\matiC.

Moreover, as the second and the third columns of $\matE$ are zero (the equilibrium values depend only on the first moment $ρ$ ), and as the first line of $\matC$ is then $(1, 0, 0)$ , we have

\matE\matC=\matE.

We have

	$\matRh$	$=\matiMh\matiTh(\matid+\matSh(\matTh\matEh\matiTh−\matid))\matTh\matMh$
		$=\matiM\matiT\matiC(\matid+\matS(\matC\matT\matE\matC\matiT\matiC−\matid))\matC\matT\matM$
		$=\matiM\matiT(\matid+\matiC\matS\matC(\matT\matE\matiT−\matid))\matT\matM.$

Then

\matRh−\matR

As the matrices $\matM$ , $\matT$ , and $\matC$ are invertible, the condition $\matRh=\matR$ is equivalent to $(\matS\matC−\matC\matS)\matT(\matE−\matid)=0$ . Denoting

\matC=⎛⎜⎝100c21c22c23c31c32c33⎞⎟⎠,

a straightforward calculation yields

(\matS\matC−\matC\matS)\matT(\matE−\matid)=(s−s′)⎛⎜⎝000c23λ2(α−6V)c236λu−c23c32λV−c320⎞⎟⎠.

Then the property $\matRh=\matR$ for all values of $s$ and $s^{'}$ is equivalent to $c_{23} = c_{32} = 0$ , that ends the proof. ∎

2.4. Positivity of the $\matR$ matrix

The velocity $V$ being fixed, we propose to give a full description of the sets

ΩV,u={(s,s′,α)∈R4 such that% \matR is a non-negative matrix},u∈R.

Indeed, the non-negativity of the matrix $\matR$ imposes that all the distributions $f_{j}$ , $1 \leq j \leq 3$ , remain non-negative if they are so at the initial time. These sets are first described by a set of nine inequalities that can be joined into just one. Numerical illustrations are then given to visualize it in the characteristic cases including SRT, MRT, and relative velocity scheme.

3. Positivity of the iterative matrix

The nine inequalities obtained from the matrix $\matR$ can be combined neatly into one formula.

The inequalities are

(1)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} R_{0, 0} & = V s u - \frac{1}{2} V s - V s^{'} u + \frac{1}{6} α s^{'} + s u - \frac{1}{2} s - s^{'} u - \frac{}{1} 6 s^{'} + 1 & ⩾ 0, R_{0, 1} & = V s u - \frac{1}{2} V s - V s^{'} u + \frac{1}{6} α s^{'} + \frac{1}{3} s^{'} & ⩾ 0, R_{0, 2} & = V s u - \frac{1}{2} V s - V s^{'} u + \frac{1}{6} α s^{'} - s u + \frac{1}{2} s + s^{'} u - \frac{1}{6} s^{'} & ⩾ 0, R_{1, 0} & = - 2 V s u + 2 V s^{'} u - \frac{1}{3} α s^{'} - 2 s u + 2 s^{'} u + \frac{1}{3} s^{'} & ⩾ 0, R_{1, 1} & = - 2 V s u + 2 V s^{'} u - \frac{1}{3} α s^{'} - \frac{1}{3} 2 s^{'} + 1 & ⩾ 0, R_{1, 2} & = - 2 V s u + 2 V s^{'} u - \frac{1}{3} α s^{'} + 2 s u - 2 s^{'} u + \frac{1}{3} s^{'} & ⩾ 0, R_{2, 0} & = V s u + \frac{1}{2} V s - V s^{'} u + \frac{1}{6} α s^{'} + s u + \frac{1}{2} s - s^{'} u - \frac{1}{6} s^{'} & ⩾ 0, R_{2, 1} & = V s u + \frac{1}{2} V s - V s^{'} u + \frac{1}{6} α s^{'} + \frac{1}{3} s^{'} & ⩾ 0, R_{2, 2} & = V s u + \frac{1}{2} V s - V s^{'} u + \frac{1}{6} α s^{'} - s u - \frac{1}{2} s + s^{'} u - \frac{1}{6} s^{'} + 1 & ⩾ 0. \end{matrix}

We prove now that the previous nine inequalities can be written in a much more lucid way.

We introduce the reduced parameters $¯ ¯ ¯ u$ and $γ$ according to

(2)

¯ ¯ ¯ u = 2 u (s - s^{'}), γ = \frac{s^{'}}{6} (1 - α) - u (s - s^{'}) V .

Then the nine previous inequalities $R_{i, j} \geq 0$ displayed in (1) are equivalent to

(3)

max (s^{'} - 1, | ¯ ¯ ¯ u |) \leq 2 γ \leq min (2 - s - | ¯ ¯ ¯ u - s V |, s - | ¯ ¯ ¯ u + s V |, s^{'} - | s V |) .

Proof.

Consider first the two inequalities associated with $R_{0, 0}$ and $R_{2, 2}$ :

{\begin{matrix} V s u - V s^{'} u + \frac{1}{6} α s^{'} - \frac{1}{2} s - \frac{1}{6} s^{'} + 1 & ⩾ \frac{1}{2} V s - s u + s^{'} u, V s u - V s^{'} u + \frac{1}{6} α s^{'} - \frac{1}{2} s - \frac{1}{6} s^{'} + 1 & ⩾ - \frac{1}{2} V s + s u - s^{'} u . \end{matrix}

They can be synthetized in the following form

\frac{1}{2} ∣ ∣ ¯ ¯ ¯ u - V s ∣ ∣ \leq 1 - \frac{1}{2} s - (\frac{1}{} 6 s^{'} (1 - α) + u (s^{'} - s) V) = 1 - \frac{1}{2} s - γ

and we can write this relation as

(4)

2 γ \leq 2 - s - | ¯ ¯ ¯ u - s V | .

Write now the inequalities (1) associated with $R_{0, 1}$ and $R_{2, 1}$ :

{\begin{matrix} V s u - V s^{'} u + \frac{1}{6} α s^{'} + \frac{1}{3} s^{'} & ⩾ \frac{1}{2} V s, V s u - V s^{'} u + \frac{1}{6} α s^{'} + \frac{1}{3} s^{'} & ⩾ - \frac{1}{2} V s . \end{matrix}

In other words, $| V s / 2 | \leq - γ + s^{'} / 2$ . Then

(5)

2 γ \leq s^{'} - | s V | .

We now focus on the inequalities (1) associated with $R_{0, 2}$ and $R_{2, 0}$ :

{\begin{matrix} V s u - V s^{'} u + \frac{1}{6} α s^{'} + \frac{1}{2} s - \frac{1}{6} s^{'} & ⩾ \frac{1}{2} V s + s u - s^{'} u, V s u - V s^{'} u + \frac{1}{6} α s^{'} + \frac{1}{2} s - \frac{1}{6} s^{'} & ⩾ - \frac{1}{2} V s - s u + s^{'} u . \end{matrix}

We have $| V s / 2 + s u - s^{'} u | \leq V s u - V s^{'} u + α s^{'} / 6 + s / 2 - s^{'} / 6 = s / 2 - γ$ . In consequence,

(6)

2 γ \leq s - | ¯ ¯ ¯ u + s V | .

Considering the inequalities (1) with $R_{1, 0}$ and $R_{1, 2}$ , we have

{\begin{matrix} - V s u + V s^{'} u - \frac{1}{6} α s^{'} + \frac{1}{6} s^{'} & ⩾ s u - s^{'} u, - V s u + V s^{'} u - \frac{1}{6} α s^{'} + \frac{1}{6} s^{'} & ⩾ - s u + s^{'} u, \end{matrix}

and $| ¯ ¯ ¯ u / 2 | \leq - V s u + V s^{'} u - α s^{'} / 6 + s^{'} / 6 = γ$ . In consequence,

(7)

| ¯ ¯ ¯ u | \leq 2 γ .

The last inequality $R_{1, 1} \geq 0$ can be written as $- V s u + V s^{'} u - α s^{'} / 6 - s^{'} / 3 + 1 / 2 ⩾ 0$ and this inequality is equivalent to $γ - s^{'} / 2 + 1 / 2 ⩾ 0$ . In other terms,

(8)

s^{'} - 1 \leq 2 γ .

The inequalities (4), (5), and (6) establish a triple majoration of $2 γ$ whereas the inequalities (7) and (8) show a double minoration of the same quantity. The proof is completed. ∎

4. The particular case $u = 0$

In this section, we suppose that the relative velocity $u$ is reduced to zero. Then the necessary and sufficient conditions (3) for the stability can be written as

(9)

max (s^{'} - 1, 0) \leq \frac{s^{'}}{3} (1 - α) \leq min (2 - s - | s V |, s - | s V |, s^{'} - | s V |) .

To fix the ideas, we suppose that the advection velocity $V$ is non-negative:

(10)

V \geq 0.

The case $V \leq 0$ follows directly. When $u = 0$ , the reduced stability conditions

(11)

max (s^{'} - 1, 0) \leq min (2 - s - | s V |, s - | s V |, s^{'} - | s V |) .

are equivalent to the following conditions for the relaxation parameters

(12)

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ \begin{matrix} 0 \leq s, s^{'} \leq 2, s^{'} \geq s V, s \leq 2 / (1 + V), s^{'} \leq min (3 - (1 + V) s, 1 + (1 - V) s) \end{matrix}

joined with a natural Courant type condition for explicit schemes on the advection velocity

(13)

V \leq 1.

Of course, the conditions (9) have still to be imposed for the equilibrium parameter $α$ when the pair $s, s^{'}$ is given. In particular,

(14)

α \leq 1

and

(15)

s^{'} \leq \frac{3}{α + 2} .

Figure 1. Necessary and sufficient stability regions described by the inequalities (11) for a null relative velocity $u$ . Illustration proposed for $V = 2 / 3$ .

Proof.

We first observe that $0 \leq max (s^{'} - 1, 0) \leq | s V | \leq min (2 - s, s)$ . Then $0 \leq s \leq 2$ . Secondly, we have $| s V | \leq s^{'}$ and because both $s$ and $V$ are positive, we have $0 \leq s V \leq s^{'}$ . We have also $s V \leq s$ and (13) is established.

Moreover, $s \geq 0$ and $| s V | \leq 2 - s$ implies $s \leq 2 / (1 + V)$ . Due to the positivity of the parameter $s$ , we deduce from (9) a new set of inequalities: $s^{'} - 1 \leq max (s^{'} - 1, 0) \leq min (2 - s - s V, s - s V)$ . Then $s^{'} \leq min (3 - (1 + V) s, 1 + (1 - V) s$ . Moreover, $s^{'} \leq 3 - (1 + V) s \leq 2$ because $V \geq 0$ .

Conversely, if the relations (12) and (13) are satisfied, we have $s^{'} - 1 \leq 2 - s - s V$ , $s^{'} - 1 \leq s - s V$ and $s^{'} - 1 \leq s^{'} - s V$ . Then $s^{'} - 1 \leq min (2 - s - s V, s - s V, s^{'} - s V$ . Moreover, $0 \leq 2 - s - s V$ , $0 \leq s - s V$ because $V \leq 1$ and $0 \leq s^{'} - s V$ . Thus $0 \leq min (2 - s - s V, s - s V, s^{'} - s V$ .

Finally the inequalities (11) are established and the proposition is proven. ∎

5. The general case

We analyse the general case of non-zero $u$ in this section, with suitable illustrations of the stability region for various ranges of the parameters.

5.1. Necessary conditions for stability

In this subsection, we prove the following proposition.

We suppose that the necessary and sufficient stability conditions (3) are satisfied:

max (s^{'} - 1, | ¯ ¯ ¯ u |) \leq 2 γ \leq min (2 - s - | ¯ ¯ ¯ u - s V |, s - | ¯ ¯ ¯ u + s V |, s^{'} - | s V |)

with the notations (2):

¯ ¯ ¯ u = 2 u (s - s^{'}), γ = \frac{s^{'}}{6} (1 - α) - u (s - s^{'}) V .

We suppose also that the advection velicity $= V$ is positive. Then if the scheme is stable, the point $(s, s^{'})$ satisfies the following inequalities

(16)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} 0 \leq s V \leq s^{'} \leq 2 0 \leq s V \leq 1 0 \leq s \leq 2 s^{'} \leq min (2 - s V, s + 1, 3 - s) s \leq 2 / (1 + V) . \end{matrix}

With these necessary stability conditions, the parameter $¯ ¯ ¯ u$ has been eliminated. We have also, necessarily

(17)

V \leq 1

and

(18)

| ¯ ¯ ¯ u | \leq \frac{1}{2} .

Proof.

We start from the inequalities (3). Then we have $0 \leq | ¯ ¯ ¯ u | \leq max (s^{'} - 1, | ¯ ¯ ¯ u |) \leq min (2 - s - | ¯ ¯ ¯ u - s V |, s - | ¯ ¯ ¯ u + s V |, s^{'} - | s V |) \leq s^{'} - | s V |$ and $s^{'} \geq | s V | \geq s V$ . With the same type of inequalities, $| ¯ ¯ ¯ u | + | ¯ ¯ ¯ u - s V | \leq 2 - s$ and $| ¯ ¯ ¯ u | + | ¯ ¯ ¯ u + s V | \leq s$ .

If we remark that $| ¯ ¯ ¯ u | \leq | ¯ ¯ ¯ u - s V | / 2 + | ¯ ¯ ¯ u + s V | / 2$ , we deduce that $| ¯ ¯ ¯ u | \leq (2 - s + s) / 2 = 1$ and the inequality (18) is proven.

We have the triangular inequality $| s V | \leq | ¯ ¯ ¯ u - s V | + | ¯ ¯ ¯ u |$ . Then from the genaral stability conditions (3), we deduce $| ¯ ¯ ¯ u - s V | + | ¯ ¯ ¯ u | \leq 2 - s$ and $| s V | \leq 2 - s$ . In a similar way, $| s V | \leq | ¯ ¯ ¯ u + s V | + | ¯ ¯ ¯ u |$ , $| ¯ ¯ ¯ u + s V | + | ¯ ¯ ¯ u | \leq s$ from (3) and finally $| s V | \leq 2 - s$ . We put together the two inequalities and we have $0 \leq | s V | \leq min (s, 2 - s)$ .

In consequence, we have $0 \leq s \leq 2$ and the third point is proven. Since we made the choice of $V \geq 0$ then $s \geq 0$ . The first inequality of the two first points are established. From $s V \leq s$ we have $V \leq 1$ and the relation (17) is true. Moreover, $s V \leq 2 - s$ and the last inequality of (16) is true.

Consider now the inequalities $s^{'} - 1 \leq 2 γ \leq min (2 - s - | ¯ ¯ ¯ u - s V |, s - | ¯ ¯ ¯ u + s V |)$ . We deduce $s^{'} - 1 + | ¯ ¯ ¯ u - s V | \leq 2 - s$ and $s^{'} - 1 + | ¯ ¯ ¯ u + s V | \leq s$ . Due to the positvity of the absolute values, we have also $s^{'} \leq 3 - s$ and $s^{'} \leq s + 1$ . A part of the fifth inequality of (16) is proven.

Finally, due to the triangular inequality, $s V = | s V | \leq | ¯ ¯ ¯ u - s V | / 2 + | ¯ ¯ ¯ u + s V | / 2$ . We add this inequality with the two following ones: $s^{'} - 1 + | ¯ ¯ ¯ u - s V | \leq 2 - s$ and $s^{'} - 1 + | ¯ ¯ ¯ u + s V | \leq s$ . Then $s^{'} - 1 + s V \leq (2 - s + s) / 2 = 1$ and $s^{'} \leq 2 - s V$ . The fifth inequality of (16) is completely established. Because $s V \geq 0$ , we have also $s^{'} \leq 2$ and the first inequality of (16) is also established. ∎

We illustrate the zones of necessary stability in the Figure 2 for five particular velocities: $V \in {0, 1 / 4, 1 / 3, 1 / 2, 2 / 3, 1}$ .

Figure 2. Necessary stability regions given by inequalities (16) for $V = 0$ (top-left), $V = 1 / 4$ (top-right), $V = 1 / 3$ (middle-left), $V = 1 / 2$ (middle-right), $V = 2 / 3$ (bottom-left), $V = 1$ (bottom-right).

5.2. Numerical study of necessary and sufficient conditions for stability

We now illustrate the necessary and sufficient stability regions for $V = 0$ , $V = 1 / 4$ , $V = 1 / 3$ , $V = 1 / 2$ , $V = 2 / 3$ , and $V = 1$ for various ranges of $u$ in the Figure 3.

Figure 3. Nunerical study of necessary and sufficient stability regions

For each value of the velocity $V$ , the necessary and sufficient stability regions given by the relations (3) are displayed for several values of the relative velocity: $u \in {- 2 V, - V, 0, V / 2, V, 2 V}$ . In dotted line, the necessary stability region of the Proposition 5.1 is added for remind. The particular value $u = 0$ is enhanced by filling the region in gray.

Some analysis can be drawn from these figures:

the stability region changes with the relative velocity;
the maximal value of the first relaxation parameter $s$ is obtained (not only) for $u = 0$ ;
the stability region is not clearly more interessant (larger or including greater value of $s$ ) for $u = V$ ;
the point $(s, s^{'}) = (1, 1)$ is always in the stability region.

To conclude this section, this notion of stability allows a large set of values for the relaxation parameters. If the scheme is used to simulate the hyperbolic advection equation without diffusive term, we can try to minimize the numerical diffusion while maintaining this stability property. This task is complicated and out of the scope of the paper as the second-order term reads as a non-linear formula that links all the parameters $s$ , $α$ , and $V$ .

6. Numerical illustrations

In this section, we illustrate the stability property with numerical simulations involving the D1Q3 model with and without relative velocity used to simulate the linear advection equation. The stability is demonstrated for the parameters chosen according to the analysis presented in the previous sections. Oscillations are seen whenever the parameters go beyond the stability limits presented, as highlighted in the results.

The parameters chosen for the simulations are the following:

	$V$	$u$	$s$	$s^{'}$	$α$
left (stable)	0.25	0.0	1.6	1.3	0.3076923076923076
	0.25	0.25	1.6	1.3	-0.17548076923076938
right (unstable)	0.25	0.0	1.9	1.4	0.14285714285714302
	0.25	0.25	1.9	1.4	-0.10491071428571441

For the left figures, the parameters are chosen in order to satisfy the stability property. We observe numerically also a maximum principle, even if it is not formally the stability notion that we investigate. For the right figures, the parameters are chosen in order not to break the stability property.

Figure 4. Smooth profile with a continuous derivative. The parameters for the D1Q3 are tuned in order to have (left) or not (right) the non-negativity property

Figure 5. Continuous profile. The parameters for the D1Q3 are tuned in order to have (left) or not (right) the non-negativity property

Figure 6. Discontinuous profile. The parameters for the D1Q3 are tuned in order to have (left) or not (right) the non-negativity property

If the profile is smooth, we have observed that no numerical oscillations occur even if the non-negativity property of the matrix is not satisfied. If the profile is just continuous, small negative values of the macroscopic quantity are observed when our non-negativity property is not satisfied. Last but not least, classical oscillations are visible for discontinuous profiles if our non-negativity property is not satisfied. These oscillations are eliminated when the non-negativity property of the matrix is realized.

7. Conclusion

In this contribution, we have investigated a stability property for a classical mono-dimensional linear three velocities lattice Boltzmann scheme with relative velocity. This property ensures that the particle distribution functions remain nonnegative if they are at initial time. We then give a necessary and sufficient condition to describe the stability region. The case without relative velocity is completely describe and simplest necessary conditions are given for the general case. We finally propose some numerical simulations that illustrate the stability property: even if the stability notion that we investigate is not exactly a cosntraint of convexity, a numerical maximum principle is observed if the parameters are in the stability region whereas numerical oscillations appear (in particular for non smooth profils) if the parameters are not there.

Moreover, relative velocities modify the stability array in a nontrivial manner. For instance, intuition might have suggested that the stability region for the relative velocity equal to the advection velocity contains all the others but it is realy not the case. For a given advection velocity, relative velocities cannot be used to increase the value of the first relaxation parameter, the one involved in the numerical diffusion.

The non-negativity of the relaxation matrix could be extended to nonliear schemes. The theoretical study will be much more technical and has not been performed. Nevertheless, numerical experiments for the Burgers equation show that the behavior of the D1Q3 scheme, and in particular the possibility or not of oscillations, is analogous to the linear case.

References

[1] F. Caetano, F. Dubois, and B. Graille. A result of convergence for a mono-dimensional two-velocities lattice Boltzmann scheme. arXiv:1905.12393 [math], May 2019. arXiv: 1905.12393.
[2] P. J. Dellar. Nonhydrodynamic modes and a priori construction of shallow water lattice Boltzmann equations. Physical Review E, 65(3):036309, February 2002.
[3] D. d’Humières. Generalized Lattice-Boltzmann equations. In Rarefied Gas Dynamics: Theory and simulation, volume 159, pages 450–458. AIAA Progress in astronomics and aeronautics, 1992.
[4] F. Dubois, T. Février, and B. Graille. Lattice Boltzmann schemes with relative velocities. Communications in Computational Physics, 17(4):1088–1112, 2015.
[5] F. Dubois, T. Février, and B. Graille. On the stability of a relative velocity lattice Boltzmann scheme for compressible Navier–Stokes equations. Comptes Rendus Mécanique, 343(10-11):599–610, 2015.
[6] F. Dubois and P. Lallemand. Towards higher order lattice Boltzmann schemes. Journal of Statistical Mechanics: Theory and Experiment, 2009:P06006, 2009.
[7] M. Geier, A. Greiner, and J. C. Korvink. Cascaded digital lattice Boltzmann automata for high reynolds number flow. Physical Review E, 73:066705, 2006.
[8] I. Ginzburg.
Truncation Errors, Exact And Heuristic Stability Analysis Of Two-Relaxation-Times Lattice Boltzmann Schemes For Anisotropic Advection-Diffusion Equation.
Communications in Computational Physics, 11(5):1439–1502, May 2012.
[9] P. Lallemand and L.-S. Luo. Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability. Physical Review E, 61:6546–6562, 2000.
[10] S. Osher. Convergence of generalized muscl schemes. SIAM Journal on Numerical Analysis, 22(5):947–961, 1985.

A stability property for a mono-dimensional three velocities scheme with relative velocity

Authors

Quasi-linear analysis of dispersion relation preservation for nonlinear schemes

Stabilized integrating factor Runge-Kutta method and unconditional preservation of maximum bound principle

Energy Stability of Explicit Runge-Kutta Methods for Non-autonomous or Nonlinear Problems

Maximum bound principle preserving integrating factor Runge-Kutta methods for semilinear parabolic equations

Energy plus maximum bound preserving Runge-Kutta methods for the Allen-Cahn equation

Type-Preserving Matrices and Security of Block Ciphers

An upwind DG scheme preserving the maximum principle for the convective Cahn-Hilliard model

1. Introduction