NUMERICAL SCHEME BASED ON FRACTIONAL
STEP METHOD IN DISCRETE KINETIC THEORY

Kokou Anani Agosseme, Amah Séna d'Almeida

► To cite this version:

Kokou Anani Agosseme, Amah Séna d'Almeida. NUMERICAL SCHEME BASED ON FRACTIONAL STEP METHOD IN DISCRETE KINETIC THEORY. International Journal of research and Reviews in Applied Sciences, 2018, 37 (1), pp.22-32. hal-05214377

HAL Id: hal-05214377

https://hal.science/hal-05214377v1

Submitted on 19 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# NUMERICAL SCHEME BASED ON FRACTIONAL STEP METHOD IN DISCRETE KINETIC THEORY

Kokou Anani Agosseme¹ & Amah Séna d'Almeida²

¹,² Université de Lomé - Département de mathématiques

¹ka.agosseme@gmail.com, ²dal_me@yahoo.fr

## ABSTRACT

In this paper we use a numerical scheme based on the classical fractional step method to solve an initial-boundary value problem resulting of the modelling of a flow in plane micro-channel by the eight velocity spatial Broadwell model. We prove the convergence of the scheme and compare the numerical solution to an analytical one. A good agreement is observed.

**Keywords:** Discrete Kinetic Theory, Fractional Step Method, Rarefied Gas, Micro-channel.

## 1. INTRODUCTION

Discrete velocity models of gas are simplified models of the Boltzmann equation but the kinetic equations are nonlinear and difficult to solve exactly. However some exact analytical solutions exist to the initial-boundary value problem [1, 2, 4, 6, 8] and to the pure boundary value problem [5, 7] and a lot of work has been done in the proof of existence of the solution to the two kinds of problem in the one-dimensional spatial case. The situation is different in the two-dimensional spatial case. For the pure boundary value problem we can mention the existence result of Cercignani and al [3]. In the unsteady case, R. Temam used in 1969 [9] a fractional step method to prove the existence of solution for the Carleman model with null Dirichlet conditions. Later in 1991 using the same method, Toscani and Walus prove the same result for the four velocity Broadwell model with specular reflection. The aim of this paper is to use a numerical scheme based on the fractional step method to solve an initial-boundary value problem for the eight velocity Broadwell model with the more realistic boundary conditions of diffuse reflection in order to investigate flows in micro-channels. The paper is organised as follows : in section 2 we state the physical problem, the scheme is then proposed and its numerical convergence is proved in the section 3, we then present some numerical results in section 4 and end with comparison with analytical results in section 5.

## 2. PROBLEM FORMULATION

Two identical semi-infinite tanks considered as two-dimensional areas said upstream and downstream regions containing the same gas are connected by a micro-channel. The physical space is related to the orthonormal reference $R = (O, x', y', z')$. The tanks respectively occupy the half-spaces $x' < 0$ and $x' > L$, ($L > 0$). The channel walls are two infinite impermeable plates located at $y = -h/2$ and $y' = h/2$, ($h > 0$). The Figure 1 shows the channel geometry.

<figure><img src="image_2.png" /><figcaption>Figure 1 Channel geometry</figcaption></figure>

The velocities of the eight velocity Broadwell model in the basis $(x', y', z')$ are $u_1 = c(-1, 1, 1)$, $u_2 = c(1, 1, 1)$, $u_3 = c(-1, -1, 1)$, $u_4 = c(1, -1, 1)$, $u_5 = c(-1, 1, -1)$, $u_6 = c(1, 1, -1)$, $u_7 = c(-1, -1, -1)$, $u_8 = c(1, -1, -1)$, where c is an
arbitrary constant. We assume that the flow is two dimensional depending only upon the spatial variables $(x', y')$ and the time $t'$. We denote by $N_i(t', x', y')$ the number density of particles of velocity $\mathbf{u}_i$ in point $M(x', y')$ and at time $t'$. The kinetic equations for this model with binary collisions are the equations (1.1)-(1.8) [5]:

$$ \frac{\partial N_1}{\partial t'} - c \frac{\partial N_1}{\partial x'} + c \frac{\partial N_1}{\partial y'} = c s \sqrt{2} (N_2 N_3 - N_1 N_4 + N_2 N_5 - N_1 N_6 + N_3 N_5 - N_1 N_7) + \frac{cs\sqrt{3}}{2} (N_2 N_7 + N_3 N_6 + N_4 N_5 - 3N_1 N_8), \quad (1.1) $$

$$ \frac{\partial N_2}{\partial t'} + c \frac{\partial N_2}{\partial x'} + c \frac{\partial N_2}{\partial y'} = c s \sqrt{2} (N_1 N_4 - N_2 N_3 + N_1 N_6 - N_2 N_5 + N_4 N_6 - N_2 N_8) + \frac{cs\sqrt{3}}{2} (N_1 N_8 + N_3 N_6 + N_4 N_5 - 3N_2 N_7), \quad (1.2) $$

$$ \frac{\partial N_3}{\partial t'} - c \frac{\partial N_3}{\partial x'} - c \frac{\partial N_3}{\partial y'} = c s \sqrt{2} (N_1 N_4 - N_2 N_3 + N_1 N_7 - N_3 N_5 + N_4 N_7 - N_3 N_8) + \frac{cs\sqrt{3}}{2} (N_2 N_7 + N_1 N_8 + N_4 N_5 - 3N_3 N_6), \quad (1.3) $$

$$ \frac{\partial N_4}{\partial t'} + c \frac{\partial N_4}{\partial x'} - c \frac{\partial N_4}{\partial y'} = c s \sqrt{2} (N_2 N_3 - N_1 N_4 + N_2 N_8 - N_4 N_6 + N_3 N_5 - N_4 N_7) + \frac{cs\sqrt{3}}{2} (N_2 N_7 + N_3 N_6 + N_1 N_8 - 3N_4 N_5), \quad (1.4) $$

$$ \frac{\partial N_5}{\partial t'} - c \frac{\partial N_5}{\partial x'} + c \frac{\partial N_5}{\partial y'} = c s \sqrt{2} (N_1 N_6 - N_2 N_5 + N_1 N_7 - N_3 N_5 + N_6 N_7 - N_5 N_8) + \frac{cs\sqrt{3}}{2} (N_2 N_7 + N_3 N_6 + N_1 N_8 - 3N_4 N_5), \quad (1.5) $$

$$ \frac{\partial N_6}{\partial t'} + c \frac{\partial N_6}{\partial x'} + c \frac{\partial N_6}{\partial y'} = c s \sqrt{2} (N_2 N_5 - N_1 N_6 + N_2 N_8 - N_4 N_6 + N_5 N_8 - N_6 N_7) + \frac{cs\sqrt{3}}{2} (N_2 N_7 + N_1 N_8 + N_4 N_5 - 3N_3 N_6), \quad (1.6) $$

$$ \frac{\partial N_7}{\partial t'} - c \frac{\partial N_7}{\partial x'} - c \frac{\partial N_7}{\partial y'} = c s \sqrt{2} (N_5 N_8 - N_6 N_7 + N_3 N_5 - N_1 N_7 + N_3 N_8 - N_4 N_7) + \frac{cs\sqrt{3}}{2} (N_1 N_8 + N_3 N_6 + N_4 N_5 - 3N_2 N_7), \quad (1.7) $$

$$ \frac{\partial N_8}{\partial t'} + c \frac{\partial N_8}{\partial x'} - c \frac{\partial N_8}{\partial y'} = c s \sqrt{2} (N_4 N_6 - N_2 N_8 + N_4 N_7 - N_3 N_8 + N_6 N_7 - N_5 N_8) + \frac{cs\sqrt{3}}{2} (N_2 N_7 + N_3 N_6 + N_4 N_5 - 3N_1 N_8). \quad (1.8) $$

Due to the symmetry of the model and that of the physical problem, we assume that $N_1=N_5$, $N_2=N_6$, $N_3=N_7$, $N_4=N_8$. The macroscopic variables of the flow are the mean density $N$, the longitudinal velocity $U$ and the transversal velocity $V$ given by:
$$ 
\begin{aligned}
N &= 2(N_1 + N_2 + N_3 + N_4), \\
NU &= 2c(-N_1 + N_2 - N_3 + N_4), \\
NV &= 2c(N_1 + N_2 - N_3 - N_4).
\end{aligned} 
\qquad (2) $$

The third component of the macroscopic velocity vanishes according to the symmetry assumption. The Maxwellian densities of the model associated with the macroscopic variables N, U and V are:

$$ 
\begin{aligned}
N_{1M} &= \frac{N}{8} \left( 1 - \frac{U}{c} \right) \left( 1 + \frac{V}{c} \right), \quad N_{2M} = \frac{N}{8} \left( 1 + \frac{U}{c} \right) \left( 1 + \frac{V}{c} \right), \\
N_{3M} &= \frac{N}{8} \left( 1 - \frac{U}{c} \right) \left( 1 - \frac{V}{c} \right), \quad N_{4M} = \frac{N}{8} \left( 1 + \frac{U}{c} \right) \left( 1 - \frac{V}{c} \right).
\end{aligned} 
\qquad (3) $$

The microscopic densities of the discrete gas in Maxwellian equilibrium with a wall, denoted $N_{lw}^{\pm}$ are the Maxwellian densities associated with 1 and the longitudinal and transversal velocities of the wall respectively denoted by $U_w^{\pm}$ and $V_w^{\pm}$. Let $\lambda^{\pm}$ the respective accommodation coefficients. The boundary conditions of diffuse reflection [7, 5] are:

$$ 
\begin{aligned}
N_1(t', x', -h/2) &= \lambda^{-}(t', x') N_{1w}^{-} = \frac{\lambda^{-}(t', x')}{8} \left( 1 - \frac{U_w^{-}}{c} \right) \left( 1 + \frac{V_w^{-}}{c} \right), \\
N_2(t', x', -h/2) &= \lambda^{-}(t', x') N_{2w}^{-} = \frac{\lambda^{-}(t', x')}{8} \left( 1 + \frac{U_w^{-}}{c} \right) \left( 1 + \frac{V_w^{-}}{c} \right), \\
N_3(t', x', h/2) &= \lambda^{+}(t', x') N_{3w}^{+} = \frac{\lambda^{+}(t', x')}{8} \left( 1 - \frac{U_w^{+}}{c} \right) \left( 1 - \frac{V_w^{+}}{c} \right), \\
N_4(t', x', h/2) &= \lambda^{+}(t', x') N_{4w}^{+} = \frac{\lambda^{+}(t', x')}{8} \left( 1 + \frac{U_w^{+}}{c} \right) \left( 1 - \frac{V_w^{+}}{c} \right).
\end{aligned} 
\qquad (4) $$

The impermeability of the plates means that the normal velocity near the plates vanishes. Therefore:

$$ \mathbf{U}^{-} \cdot \mathbf{n}^{-} = 0, \quad \mathbf{U}^{+} \cdot \mathbf{n}^{+} = 0 \qquad (5) $$

where $\mathbf{n}^{-}$ and $\mathbf{n}^{+}$ denote the inward-pointing (i.e. into the gas) unit vectors normal to the plates and $\mathbf{U}^{-}$ and $\mathbf{U}^{+}$ the velocities of the discrete gas at M(x', -h/2) and M(x', h/2) respectively. At the inlet (x' = 0) and the outlet (x' = L) of the channel we assume that the microscopic densities of the molecules are known and are equal to the microscopic densities of the gas molecules of upstream (x' < 0) and downstream (x' > L) regions respectively. These are denoted $N_{l,0}$ and $N_{l,L}$ respectively. Let $N_0, U_0, V_0$ and $N_L, U_L, V_L$ be the macroscopic variables respectively in the upstream and downstream regions. We have:

$$ 
\begin{aligned}
N_1(t', L, y') &= N_{1,L} = \frac{N_L}{8} \left( 1 - \frac{U_L}{c} \right) \left( 1 + \frac{V_L}{c} \right), \quad N_2(t', 0, y') = N_{2,0} = \frac{N_0}{8} \left( 1 + \frac{U_0}{c} \right) \left( 1 + \frac{V_0}{c} \right) \\
N_3(t', L, y') &= N_{3,L} = \frac{N_L}{8} \left( 1 - \frac{U_L}{c} \right) \left( 1 - \frac{V_L}{c} \right), \quad N_4(t', 0, y') = N_{4,0} = \frac{N_0}{8} \left( 1 + \frac{U_0}{c} \right) \left( 1 - \frac{V_0}{c} \right).
\end{aligned} 
\qquad (6) $$

We assume that, initially, the gas in the channel has the macroscopic variables of the gas in the downstream region. So:

$$ N_l(0, x', y') = N_{l,L}, \quad l \in \{1,2,3,4\}. \qquad (7) $$
The problem is put in the dimensionless form using the following dimensionless quantities :

$$
\begin{align}
x &= x'/L, y &= y'/h, t &= ct'/h, \varepsilon &= h/L, Kn &= (sN_0h)^{-1}, n_l &= N_l/N_0, n_{lw}^{\pm} &= N_{lw}^{\pm}, \nonumber \\
\rho_0 &= N_0/N_0 = 1, \rho_L &= N_L/N_0, u_w^{\pm} &= U_w^{\pm}/c, v_w^{\pm} &= V_w^{\pm}/c, u_0 &= U_0/c, v_0 &= V_0/c, \tag{8} \\
u_L &= U_L/c, v_L &= V_L/c, u &= U/c, v &= V/c. \nonumber
\end{align}
$$

The dimensionless problem to solve is :

$$
\begin{align}
& \frac{\partial n_1}{\partial t} - \varepsilon \frac{\partial n_1}{\partial x} + \frac{\partial n_1}{\partial y} = \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_2 n_3 - n_1 n_4), \notag \\
& \frac{\partial n_2}{\partial t} + \varepsilon \frac{\partial n_2}{\partial x} + \frac{\partial n_2}{\partial y} = \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_1 n_4 - n_2 n_3), \notag \\
& \frac{\partial n_3}{\partial t} - \varepsilon \frac{\partial n_3}{\partial x} - \frac{\partial n_3}{\partial y} = \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_1 n_4 - n_2 n_3), \notag \\
& \frac{\partial n_4}{\partial t} + \varepsilon \frac{\partial n_4}{\partial x} - \frac{\partial n_4}{\partial y} = \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_2 n_3 - n_1 n_4) \notag \\
& n_l(0, x, y) = n_{l,L}, \quad (x, y) \in [0,1] \times [-1/2,1/2], \quad l = 1,2,3,4 \notag \\
& n_l(t, 0, y) = n_{l,0}, \quad y \in [-1/2,1/2], \quad l \in \{2,4\} \notag \\
& n_l(t, 1, y) = n_{l,L}, \quad y \in [-1/2,1/2], \quad l \in \{1,3\} \tag{9} \\
& n_1(t, x, -1/2) = \lambda^-(t, x) n_{1w}^-, \quad x \in [0,1], \notag \\
& n_2(t, x, -1/2) = \lambda^-(t, x) n_{2w}^-, \quad x \in [0,1], \notag \\
& n_3(t, x, 1/2) = \lambda^+(t, x) n_{3w}^+, \quad x \in [0,1], \notag \\
& n_4(t, x, 1/2) = \lambda^+(t, x) n_{4w}^+, \quad x \in [0,1], \notag \\
& n_1(t, x, -1/2) + n_2(t, x, -1/2) - n_3(t, x, -1/2) - n_4(t, x, -1/2) = 0, \quad x \in ]0,1[ \notag \\
& n_1(t, x, 1/2) + n_2(t, x, 1/2) - n_3(t, x, 1/2) - n_4(t, x, 1/2) = 0, \quad x \in ]0,1[. \notag
\end{align}
$$

3. NUMERICAL SCHEME

We present in this section the numerical scheme used to solve the problem. The time step is Δt and $n_l^m$ is the density $n_l$ at time $t=m\Delta t$, ($m=0,1,2,...$), $n_l^{m+1/2}$ the density in the middle time:

$$
\begin{align}
\frac{n_1^{m+1/2} - n_1^m}{\Delta t} &= \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_2^{m+1/2} n_3^{m+1/2} - n_1^{m+1/2} n_4^{m+1/2}), \quad (10.1) \\
\frac{n_2^{m+1/2} - n_2^m}{\Delta t} &= \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_1^{m+1/2} n_4^{m+1/2} - n_2^{m+1/2} n_3^{m+1/2}), \quad (10.2) \\
\frac{n_3^{m+1/2} - n_3^m}{\Delta t} &= \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_1^{m+1/2} n_4^{m+1/2} - n_2^{m+1/2} n_3^{m+1/2}), \quad (10.3) \\
\frac{n_4^{m+1/2} - n_4^m}{\Delta t} &= \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_2^{m+1/2} n_3^{m+1/2} - n_1^{m+1/2} n_4^{m+1/2}), \quad (10.4)
\end{align}
$$
$$
\begin{align}
\frac{n_1^{m+1} - n_1^{m+1/2}}{\Delta t} - \varepsilon \frac{\partial n_1^{m+1}}{\partial x} + \frac{\partial n_1^{m+1}}{\partial y} &= 0, \quad (11.1) \\
\frac{n_2^{m+1} - n_2^{m+1/2}}{\Delta t} + \varepsilon \frac{\partial n_2^{m+1}}{\partial x} + \frac{\partial n_2^{m+1}}{\partial y} &= 0, \quad (11.2), \\
\frac{n_3^{m+1} - n_3^{m+1/2}}{\Delta t} - \varepsilon \frac{\partial n_3^{m+1}}{\partial x} - \frac{\partial n_3^{m+1}}{\partial y} &= 0, \quad (11.3) \\
\frac{n_4^{m+1} - n_4^{m+1/2}}{\Delta t} + \varepsilon \frac{\partial n_4^{m+1}}{\partial x} - \frac{\partial n_4^{m+1}}{\partial y} &= 0, \quad (11.4)
\end{align}
$$

$$
\begin{align}
& n_{2w}^{-} n_{1}^{m+1} - n_{1w}^{-} n_{2}^{m+1} = 0, \quad y = -1/2, x \in [0,1], \notag \\
& n_{4w}^{+} n_{3}^{m+1} - n_{3w}^{+} n_{4}^{m+1} = 0, \quad y = 1/2, x \in [0,1], \notag \\
& \nu^{m+1} = 0, \quad y = \pm 1/2, x \in ]0,1[, \tag{12} \\
& n_{l}^{m+1} = n_{l,0}, \quad x = 0, \quad y \in [-1/2,1/2], l \in \{2,4\} \notag \\
& n_{l}^{m+1} = n_{l,L}, \quad x = 1, \quad y \in [-1/2,1/2], l \in \{1,3\}. \notag
\end{align}
$$

We recall that the system (10) represents the spatially homogeneous part of the flow problem and the system (11)
the free molecular regime part. The equations (12) are the boundary conditions taken at the time $t=(m+1)\Delta t$. We
deduce from equations (10), the relations:

$$
\begin{align}
n_1^{m+1/2} &= \frac{n_1^m + \alpha(n_1^m + n_2^m)(n_1^m + n_3^m)}{1 + \alpha(n_1^m + n_2^m + n_3^m + n_4^m)}, & n_2^{m+1/2} &= \frac{n_2^m + \alpha(n_1^m + n_2^m)(n_2^m + n_4^m)}{1 + \alpha(n_1^m + n_2^m + n_3^m + n_4^m)}, \tag{13} \\
n_3^{m+1/2} &= \frac{n_3^m + \alpha(n_1^m + n_3^m)(n_3^m + n_4^m)}{1 + \alpha(n_1^m + n_2^m + n_3^m + n_4^m)}, & n_4^{m+1/2} &= \frac{n_4^m + \alpha(n_2^m + n_4^m)(n_3^m + n_4^m)}{1 + \alpha(n_1^m + n_2^m + n_3^m + n_4^m)}, \notag
\end{align}
$$

where $\alpha = (\sqrt{2} + \sqrt{3})\Delta t / Kn$. The quantities $n_l^m$ and $n_l^{m+1/2}$ depend upon x and y. We perform a regular grid of the domain $[0,1]\times[-1/2,1/2]$ with the steps $\Delta x = 1/(J-1)$ and $\Delta y = 1/(K-1)$ where $J, K \in \mathbb{N} \setminus \{0,1\}$. Let $n_{l,j,k}^{m+1}$ be the value of $n_l^{m+1}$ at the point $(x_j, y_k) \in [0,1] \times [-1/2,1/2]$.

$$
\begin{align}
\frac{n_{1,j,k}^{m+1} - n_{1,j,k}^{m+1/2}}{\Delta t} - \varepsilon \frac{n_{1,j+1,k}^{m+1} - n_{1,j,k}^{m+1}}{\Delta x} + \frac{n_{1,j,k}^{m+1} - n_{1,j;k-1}^{m+1}}{\Delta y} &= 0, \quad j=1, \dots, J-1, k=2, \dots, K \notag \\
\frac{n_{2,j,k}^{m+1} - n_{2,j,k}^{m+1/2}}{\Delta t} + \varepsilon \frac{n_{2,j,k}^{m+1} - n_{2,j-1;k}^{m+1}}{\Delta x} + \frac{n_{2,j,k}^{m+1} - n_{2,j;k-1}^{m+1}}{\Delta y} &= 0, \quad j=2, \dots, J, k=2, \dots, K, \notag \\
\frac{n_{3,j,k}^{m+1} - n_{3,j,k}^{m+1/2}}{\Delta t} - \varepsilon \frac{n_{3,j+1,k}^{m+1} - n_{3,j;k}^{m+1}}{\Delta x} - \frac{n_{3,j,k+1}^{m+1} - n_{3,j;k}^{m+1}}{\Delta y} &= 0, \quad j=1, \dots, J-1, k=1, \dots, K-1 \notag \\
\frac{n_{4,j,k}^{m+1} - n_{4,j,k}^{m+1/2}}{\Delta t} + \varepsilon \frac{n_{4,j,k}^{m+1} - n_{4,j-1;k}^{m+1}}{\Delta x} - \frac{n_{4,j,k+1}^{m+1} - n_{4,j;k}^{m+1}}{\Delta y} &= 0, \quad j=2, \dots, J, k=1, \dots, K-1 \notag
\end{align}
$$
$$
\begin{align*}
n_{2,1,k}^{m+1} &= n_{2,0}, \quad n_{4,1,k}^{m+1} = n_{4,0}, \quad k = 1, \dots, K \\
n_{1,J,k}^{m+1} &= n_{1,L}, \quad n_{3,J,k}^{m+1} = n_{3,L}, \quad k = 1, \dots, K \\
n_{2w}^{-} n_{1,j,l}^{m+1} - n_{1w}^{-} n_{2,j,l}^{m+1} &= 0, \quad j = 1, \dots, J \\
n_{4w}^{+} n_{3,j,K}^{m+1} - n_{3w}^{+} n_{4,j,K}^{m+1} &= 0, \quad j = 1, \dots, J \\
n_{1,j,l}^{m+1} + n_{2,j,l}^{m+1} - n_{3,j,l}^{m+1} - n_{4,j,l}^{m+1} &= 0, \quad j = 2, \dots, J-1 \\
n_{1,j,K}^{m+1} + n_{2,j,K}^{m+1} - n_{3,j,K}^{m+1} - n_{4,j,K}^{m+1} &= 0, \quad j = 2, \dots, J-1.
\end{align*}
\tag{15}
$$

**3.1. Consistency**

Consider the equation (10.1) and the equation (11.1). By addition one has:

$$
\frac{n_{1,j,k}^{m+1} - n_{1,j,k}^{m}}{\Delta t} - \varepsilon \frac{n_{1,j+1,k}^{m+1} - n_{1,j,k}^{m+1}}{\Delta x} + \frac{n_{1,j,k}^{m+1} - n_{1,j,k-1}^{m+1}}{\Delta y} = \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_{2}^{m+1/2} n_{3}^{m+1/2} - n_{1}^{m+1/2} n_{4}^{m+1/2}). \quad (16)
$$

Making a Taylor serie expansion we have:

$$
\begin{align}
\frac{n_{1,j,k}^{m+1} - n_{1,j,k}^{m}}{\Delta t} &= \frac{\partial n_1}{\partial t}(t_{m+1}, x_j, y_k) + O(\Delta t), \nonumber \\
\frac{n_{1,j+1,k}^{m+1} - n_{1,j,k}^{m+1}}{\Delta x} &= \frac{\partial n_1}{\partial x}(t_{m+1}, x_j, y_k) + O(\Delta x), \tag{17} \\
\frac{n_{1,j,k}^{m+1} - n_{1,j,k-1}^{m+1}}{\Delta y} &= \frac{\partial n_1}{\partial y}(t_{m+1}, x_j, y_k) + O(\Delta y). \nonumber
\end{align}
$$

Then

$$
\left( \frac{n_{1,j,k}^{m+1} - n_{1,j,k}^{m}}{\Delta t} - \varepsilon \frac{n_{1,j+1,k}^{m+1} - n_{1,j,k}^{m+1}}{\Delta x} + \frac{n_{1,j,k}^{m+1} - n_{1,j,k-1}^{m+1}}{\Delta y} \right) - \left( \frac{\partial n_1}{\partial t}(t_{m+1}, x_j, y_k) - \varepsilon \frac{\partial n_1}{\partial x}(t_{m+1}, x_j, y_k) + \frac{\partial n_1}{\partial y}(t_{m+1}, x_j, y_k) \right) = O(\Delta t + \Delta x + \Delta y). \quad (18)
$$

The same argument holds for $l \in \{2,3,4\}$. We thus conclude that the scheme is accurate of order 1 in time and space.

**3.2. Stability**

We use Fourier analysis to study the stability of the scheme. We put:

$$
\begin{align}
n_{l,j,k}^{m} &= \tilde{n}_{l}^{m}(\xi, \eta) \exp(i(\xi j \Delta x + \eta k \Delta y)), \nonumber \\
\rho_{jk}^{m} &= 2 \sum_{l=1}^{4} n_{l,j,k}^{m} = \tilde{\rho}^{m}(\xi, \eta) \exp(i(\xi j \Delta x + \eta k \Delta y)) \tag{19}
\end{align}
$$

with $\tilde{\rho}^m(\xi, \eta) = 2(\tilde{n}_1^m(\xi, \eta) + \tilde{n}_2^m(\xi, \eta) + \tilde{n}_3^m(\xi, \eta) + \tilde{n}_4^m(\xi, \eta))$, where $(\xi, \eta)$ is an arbitrary wave vector and $i$ is the complex number such that $i^2 = -1$. The boundedness of $\tilde{n}_l^m(\xi, \eta)$, $l \in \{1, 2, 3, 4\}$ is equivalent to
that of $\tilde{\rho}^m(\xi, \eta)$. Using the conservation of mass in equations (10), one can write $\tilde{\rho}^{m+1/2}(\xi, \eta) = \tilde{\rho}^m(\xi, \eta)$. We have:

$$ 
\begin{aligned}
n_{l,j-1,k}^m &= n_{l,j,k}^m \exp(-i\xi\Delta x), \\
n_{l,j+1,k}^m &= n_{l,j,k}^m \exp(i\xi\Delta x), \\
n_{l,j,k-1}^m &= n_{l,j,k}^m \exp(-i\eta\Delta y), \\
n_{l,j,k+1}^m &= n_{l,j,k}^m \exp(i\eta\Delta y).
\end{aligned} 
\quad (20) $$

We replace these relations in the equations (14) to obtain:

$$ F_l(\xi, \eta)n_{l,j,k}^m = n_{l,j,k}^{m+1/2}, \quad l \in \{1,2,3,4\} \quad (21) $$

with

$$ 
\begin{aligned}
F_1(\xi, \eta) &= 1 + \sigma + \gamma - \gamma \exp(i\xi\Delta x) - \sigma \exp(-i\eta\Delta y), \\
F_2(\xi, \eta) &= 1 + \sigma + \gamma - \gamma \exp(-i\xi\Delta x) - \sigma \exp(-i\eta\Delta y), \\
F_3(\xi, \eta) &= 1 + \sigma + \gamma - \gamma \exp(i\xi\Delta x) - \sigma \exp(i\eta\Delta y), \\
F_4(\xi, \eta) &= 1 + \sigma + \gamma - \gamma \exp(-i\xi\Delta x) - \sigma \exp(i\eta\Delta y).
\end{aligned} 
\quad (22) $$

and $\gamma = \varepsilon\Delta t / \Delta x$ and $\sigma = \Delta t / \Delta y$. By taking the modulus, we can write:

$$ 
\begin{aligned}
|F_1(\xi, \eta)|^2 &= [1 + \gamma(1 - \cos(\xi\Delta x)) + \sigma(1 - \cos(\eta\Delta y))]^2 + [-\gamma \sin(\xi\Delta x) + \sigma \sin(\eta\Delta y)]^2, \\
|F_2(\xi, \eta)|^2 &= [1 + \gamma(1 - \cos(\xi\Delta x)) + \sigma(1 - \cos(\eta\Delta y))]^2 + [\gamma \sin(\xi\Delta x) + \sigma \sin(\eta\Delta y)]^2, \\
|F_3(\xi, \eta)|^2 &= [1 + \gamma(1 - \cos(\xi\Delta x)) + \sigma(1 - \cos(\eta\Delta y))]^2 + [-\gamma \sin(\xi\Delta x) - \sigma \sin(\eta\Delta y)]^2, \\
|F_4(\xi, \eta)|^2 &= [1 + \gamma(1 - \cos(\xi\Delta x)) + \sigma(1 - \cos(\eta\Delta y))]^2 + [\gamma \sin(\xi\Delta x) - \sigma \sin(\eta\Delta y)]^2.
\end{aligned} 
\quad (23) $$

As for any $X \in \mathbb{R}, 1 - \cos(X) \ge 0$, we have $|F_l(\xi, \eta)| > 1, \quad l \in \{1,2,3,4\}$. Thus all the amplification factors $1/F_l(\xi, \eta)$ satisfy $|1/F_l(\xi, \eta)| < 1$. Then:

$$ n_l^{m+1}(\xi, \eta) \le n_l^{m+1/2}(\xi, \eta), \quad l \in \{1,2,3,4\}. \quad (24) $$

Making the sum, we have:

$$ 
\begin{aligned}
\rho^{m+1}(\xi, \eta) &\le \rho^{m+1/2}(\xi, \eta), \quad \forall m \\
&\le \rho^m(\xi, \eta), \quad \forall m
\end{aligned} 
\quad (25) $$

Finally

$$ \rho^m(\xi, \eta) \le \rho^0(\xi, \eta), \quad \forall m. \quad (26) $$

We can then conclude to the stability of the scheme and therefore to its convergence.
4. NUMERICAL RESULTS

We compute numerically the flow in a channel with high to length ratio $\varepsilon = 0.1$ and $Kn = 0.05$. The gas in the upstream and downstream regions are at rest. The plates are also at rest. So $u_0 = v_0 = u_L = v_L = u_w^\pm = v_w^\pm = 0$. The flow is induced by density gradient such as $\rho_0 = 1 > \rho_L = 0.5$. The time step is $\Delta t=0.001$ and $J=K=21$. We give some results at the steady state.

When $Kn$ tends towards zero we have a Poiseuille flow. We have the non slip conditions at the plates (Figure 2(a)). For higher value of $Kn$ we have velocity slips whose thickness depends upon the Knudsen number. For $Kn$ tending towards $+\infty$ the velocity tends towards a constant value (Figure 2(b)). The velocity slip increases with $Kn$. However its rate of variation is not uniform. The variation is rapid for low and transitional Knudsen numbers and weak for large Knudsen numbers. The Figure 3(a) and 3(b) shows the dependence of the profiles of the macroscopic velocity vector field with Knudsen number: the profiles are parabolic for low and transitional Knudsen number and linear for higher Knudsen number.

<figure><img src="image_3.png" /><figcaption>Figure 2. Velocity profile at the steady state : (a) Longitudinal velocity as a function of y for x= 0.5, (b) Velocity slip</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3. Velocity vector U(u, v) at the steady state: (a) Kn=0.05 (b) Kn=10</figcaption></figure>

5. COMPARISON WITH ANALYTICAL SOLUTION

To validate the scheme, we compare the numerical results to those of the analytical solution for Couette flow. We consider a Couette flow between two parallel plates and we solve the boundary value problem in the steady state. The notation are the same as in the above. We suppose here that $u_w^{\pm} \neq 0$, and $v_w^{\pm} = 0$. As it is usual in the study of Couette flow [5, 8], the flow is supposed to be one dimensional and depends only on the variable y. The mathematical problem is:

$$ 
\begin{align} 
\frac{\partial n_1}{\partial y} &= \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_2 n_3 - n_1 n_4), \nonumber \\ 
\frac{\partial n_2}{\partial y} &= \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_1 n_4 - n_2 n_3), \nonumber \\ 
-\frac{\partial n_3}{\partial y} &= \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_1 n_4 - n_2 n_3), \nonumber \\ 
-\frac{\partial n_4}{\partial y} &= \frac{\sqrt{2} + \sqrt{3}}{Kn} (n_2 n_3 - n_1 n_4) \tag{27} 
\end{align} 
$$

$$ 
\begin{align*} 
n_1(-1/2) &= \lambda^{-}(1-u_w^{-})/8, \quad n_2(-1/2) = \lambda^{-}(1+u_w^{-})/8, \\ 
n_3(1/2) &= \lambda^{+}(1-u_w^{+})/8, \quad n_4(1/2) = \lambda^{+}(1+u_w^{+})/8, \\ 
n_1(-1/2) + n_2(-1/2) - n_3(-1/2) - n_4(-1/2) &= 0, \\ 
n_1(1/2) + n_2(1/2) - n_3(1/2) - n_4(1/2) &= 0. 
\end{align*} 
$$

The exact analytical solution $(n_1, n_2, n_3, n_4)$ of this equations is given by:

$$ n_1(y) = \frac{\beta k_2}{16} y + k_1, \quad n_2(y) = \frac{1}{4} - n_1(y), \quad n_3(y) = \frac{k_2}{4} + n_1(y), \quad n_4(y) = \frac{1}{4} - \frac{k_2}{4} + n_1(y), \quad (28) $$
where

$$k_1 = \frac{1-u_w^-}{8} + \frac{\beta(u_w^- - u_w^+)}{16\beta + 64}, \quad k_2 = \frac{2(u_w^- - u_w^+)}{\beta + 4}, \quad \beta = \frac{\sqrt{2} + \sqrt{3}}{Kn}. \qquad (29)$$

Then the longitudinal velocity u is:

$$u(y) = 2[-n_1(y) + n_2(y) - n_3(y) + n_4(y)] = \frac{\beta(u_w^+ - u_w^-)}{\beta + 4} y + \frac{u_w^- + u_w^+}{2}. \qquad (30)$$

We perform the numerical computations with $\varepsilon = 0.1, u_0 = v_0 = u_L = v_L = v_w^\pm = 0, Kn = 0.05$ and

$u_w^- = -u_w^+ = -0.2$. The numerical value of the mean density is 1 in the upstream and downstream region. The flow is induced by the movement of the wall. The comparison is made for $x = 0.5$ that is at the half of length of the channel where the end effects would not affect the flow which can be considered as one dimensional shear flow depending merely on y. We find a very good agreement of analytical and numerical results as shown on the figure 4.

<figure><img src="image_5.png" /><figcaption>Figure 4: Longitudinal velocity comparison (a) Velocity vector U(u, v) (b) Numerical and exact solution</figcaption></figure>

## 6. CONCLUSIONS

We solve an initial-boundary value problem with a numerical scheme based on the fractional step method. We show that the scheme converges. To test its accuracy we compare the numerical results obtained with this method to analytical ones in the case of Couette flows. A good agreement is obtained. Moreover the method can be easily used for more complex models with at least two moduli of velocity in order to study flows including thermal process.
7. REFERENCES

[1] A. V. Bobylev, On the nonlinear BE and exact solution of the carlemen model system, *Congress of Mathematicians Abstracts* **29**, Warsaw, Poland, (1983).

[2] H. Cabannes, The Discrete Boltzmann Equation (Theory and Applications), *Lecture Note given at the University of California at Berkeley during the Springer Quarter*, 1980.

[3] C. Cercignani, R. Illner, and M. Shinbrot, A boundary value problem for the two dimensional broadwell model, *Commun. Math. Phys.*, **114**, 687–698 (1988).

[4] H. Cornille and T. Platkowski, Exact solutions of a hierarchy of mixing speed models, *J. Math. Phys.*, **33**, 2587–2606 (1992).

[5] A. d'Almeida, Étude des solutions des équations de Boltzmann discrétes et applications. PhD thesis, Thèse de l'Université Pierre et Marie Curie, Paris, 24 février 1995.

[6] A. d'Almeida, Exact solutions for discrete velocity models, *Mechanics Research Communications*, **34**, 405–409 (2007).

[7] A. d'Almeida and R. Gatignol, Boundary conditions for discrete models of gases and applications to Couette flows. In Springer-Verlag Eds, editor, Computational Fluid Dynamics, 115–130. D. Leutloff, R.C. Srivastava, (1995).

[8] T. Natta, K. A. Agosseme, and A. S. d'Almeida, Exact solution of the four velocity broadwell model. Global Journal of Pure and Applied Mathematics, **13**(10), 7035–7050 (2017).

[9] R. Temam, Sur la résolution exacte et approchée d'un problème hyperbolique non linéaire de T, *Carleman. Arch. Rat. Mech. Anal.*, **35**, 351–362 (1969).