Difference Finite Scheme derived from MRT Lattice Boltzmann for Convective-Diffusive Equation

Based on the works of Bellotti et al. [7] and Chen et al. [8], we will develop a difference finite description of lattice Boltzmann equation in the multi-relaxation-times (MRT) framework.

MRT Lattice Boltzmann equation for Convective-Diffusive Equation

MRT Lattice Boltzmann equation for Convective-Diffusive Equation

To describe convective-diffusive problems, the MRT lattice Boltzmann formulation is described by:

(62)\[ f_i( x_{\alpha} + e_{i,\alpha} \delta t, t+\delta t) = f_i(x_{\alpha}, t) - (\textbf{M}^{-1}\textbf{S}\textbf{M})_{ik}\left( f_{k} - f_{k}^{eq} \right), \]

where $\textbf{M}$ is matrix of moments, $\textbf{S}$ is the diagonal matrix of relaxation times, and $\textbf{M}^{-1}$ is the inverse matrix of $\textbf{M}$. Considering the problem description in a lattice D1Q3, we have:

\[\begin{split} \textbf{M} =\begin{pmatrix} 1 \\ e_{i,x} \\ 3e_{i,x}^{2}-2 \end{pmatrix} =\begin{pmatrix} 1 & 1 & 1 \\ 0 & 1 & -1 \\ -2 & 1 & 1 \end{pmatrix}, \qquad \qquad \textbf{S} = \begin{pmatrix} s_0 & 0 & 0 \\ 0 & s_1 & 0 \\ 0 & 0 & s_2 \end{pmatrix}, \qquad \qquad \textbf{M}^{-1}=\begin{pmatrix}\dfrac{1}{3} & 0 & -\dfrac{1}{3} \\ \dfrac{1}{3} & \dfrac{1}{2} & \dfrac{1}{6} \\ \dfrac{1}{3} & -\dfrac{1}{2} & \dfrac{1}{6} \end{pmatrix}. \end{split}\]

The equilibrium distribution function is defined by

\[ f^{eq}_{i} = w_{i}\left( u + \frac{e_{i,\alpha} B_{\alpha}}{c_{s}^{2}} + \frac{ \left( c_{s}^{2}C\delta_{\alpha\beta} + \Lambda_{\alpha\beta} - c_{s}^{2}u\delta_{\alpha\beta} \right) \left( e_{i,\alpha}e_{i,\beta} - c_{s}^{2}\delta_{\alpha\beta} \right)}{2c_{s}^{4}} \right), \]

where $\Lambda_{\alpha\beta}$ is a correction term defined by

\[ \Lambda_{\alpha\beta} = \displaystyle\int B_{\alpha}^{'}B_{\beta}^{'} du = \displaystyle\int \left(\displaystyle\frac{\partial B_{\alpha}}{\partial u} \displaystyle\frac{\partial B_{\beta}}{\partial u} \right) du. \]

The different orders of equilibrium moments can be recovered by:

\[\displaystyle \underbrace{\sum_{i=0} f_{i}^{eq} =\sum_{i=0} f_{i} }_{\textrm{Zero-Order Moment}} =u,\quad \quad \underbrace{\sum_{i=0} f_{i}^{eq}e_{i,x} }_{\textrm{x-First-Order Moment}} =B_{\alpha},\quad \quad \underbrace{\sum_{i=0} f_{i}^{eq} e_{i,x}e_{i,x} }_{\textrm{xx-Second-Order Moment}} =\frac{C}{3} + \Lambda_{\alpha\beta},\quad \quad \underbrace{\sum_{i=0} f_{i}^{eq} e_{i,x}e_{i,x}e_{i,x} }_{\textrm{xxx-Third-Order Moment}} =B_{\alpha}\]

\[\displaystyle \underbrace{\sum_{i=0} f_{i}^{eq} =\sum_{i=0} f_{i} }_{\textrm{Zero-Order Hermite Moment}} =u,\quad \quad \underbrace{\sum_{i=0} f_{i}^{eq}e_{i,x} }_{\textrm{x-First-Order Hermite Moment}} =B_{\alpha},\quad \quad \underbrace{\sum_{i=0} f_{i}^{eq} \left(e_{i,x}e_{i,x} - \frac{1}{c_{s}^{2}}\right) }_{\textrm{xx-Second-Order Hermite Moment}} =\frac{C}{3} + \Lambda_{\alpha\beta} - \frac{u}{3},\quad \quad \underbrace{\sum_{i=0} f_{i}^{eq} \left(e_{i,x}e_{i,x} - \frac{3}{c_{s}^{2}}\right) e_{i,x} }_{\textrm{xxx-Third-Order Hermite Moment}} =0\]

FDM Interpretation in MRT Framework

Reinterpreting the Eq. (62) with a shift in distribution function direction:

\[ f_{i,x}^{t+1} = f_{i,x-e_{i}}^{t} - (\textbf{M}^{-1}\textbf{S}\textbf{M})_{ik}\left( f_{k,x-e_{i}}^{t} - f_{k,x-e_{i}}^{eq,t} \right) = (\textbf{M}^{-1} (\textbf{I} -\textbf{S}) \textbf{M})_{ik} f_{k,x-e_{i}}^{t} + (\textbf{M}^{-1}\textbf{S}\textbf{M})_{ik} f_{k,x-e_{i}}^{eq,t} , \]

where $f_i( x_{\alpha} + e_{i,\alpha} \delta t, t+\delta t)=f_{i,x+c_{i}}^{t+1}$ for simplificate the notation and $\textbf{I}$ is identity matrix. Defining the shift operator $T_{\Delta x}^{c_{i}}$, where $f_{k,x-e_{i}}^{t}= T_{\Delta x}^{c_{i}}[f_{k,x}^{t}]$, we can rewrite the above equation inside of the moments space:

(63)\[ \textbf{m}_{k,x}^{t+1} = \textbf{P}\textbf{m}_{k,x}^{t} + \textbf{Q}\textbf{m}_{k,x}^{eq,t}, \]

where $\textbf{m}_{k,x}^{t}=M_{ik}f_{i,x}^{t}$, $k$ is the index of moments in the space moment, $\textbf{T}:= \textbf{M} (diag\left( T_{\Delta x}^{c_{0}}, T_{\Delta x}^{c_{1}}, T_{\Delta x}^{c_{2}} \right)) \textbf{M}^{-1}$, $\textbf{P}:=\textbf{T}(\textbf{I}-\textbf{S})$, and , $\textbf{Q}:=\textbf{T}\textbf{S}$. The matrices are given by:

\[\begin{split}\displaystyle \textbf{T} =\left[\begin{matrix}\frac{T_{\Delta_{x}}^{c_{0}}}{3} + \frac{T_{\Delta_{x}}^{c_{1}}}{3} + \frac{T_{\Delta_{x}}^{c_{2}}}{3} & \frac{T_{\Delta_{x}}^{c_{1}}}{2} - \frac{T_{\Delta_{x}}^{c_{2}}}{2} & - \frac{T_{\Delta_{x}}^{c_{0}}}{3} + \frac{T_{\Delta_{x}}^{c_{1}}}{6} + \frac{T_{\Delta_{x}}^{c_{2}}}{6}\\\frac{T_{\Delta_{x}}^{c_{1}}}{3} - \frac{T_{\Delta_{x}}^{c_{2}}}{3} & \frac{T_{\Delta_{x}}^{c_{1}}}{2} + \frac{T_{\Delta_{x}}^{c_{2}}}{2} & \frac{T_{\Delta_{x}}^{c_{1}}}{6} - \frac{T_{\Delta_{x}}^{c_{2}}}{6}\\- \frac{2 T_{\Delta_{x}}^{c_{0}}}{3} + \frac{T_{\Delta_{x}}^{c_{1}}}{3} + \frac{T_{\Delta_{x}}^{c_{2}}}{3} & \frac{T_{\Delta_{x}}^{c_{1}}}{2} - \frac{T_{\Delta_{x}}^{c_{2}}}{2} & \frac{2 T_{\Delta_{x}}^{c_{0}}}{3} + \frac{T_{\Delta_{x}}^{c_{1}}}{6} + \frac{T_{\Delta_{x}}^{c_{2}}}{6}\end{matrix}\right],\end{split}\]

\[\begin{split}\displaystyle \textbf{P} =\left[\begin{matrix}- \frac{\left(s_{0} - 1\right) \left(T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{3} & - \frac{\left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right) \left(s_{1} - 1\right)}{2} & - \frac{\left(s_{2} - 1\right) \left(- 2 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{6}\\- \frac{\left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right) \left(s_{0} - 1\right)}{3} & - \frac{\left(T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right) \left(s_{1} - 1\right)}{2} & - \frac{\left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right) \left(s_{2} - 1\right)}{6}\\- \frac{\left(s_{0} - 1\right) \left(- 2 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{3} & - \frac{\left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right) \left(s_{1} - 1\right)}{2} & - \frac{\left(s_{2} - 1\right) \left(4 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{6}\end{matrix}\right],\end{split}\]

\[\begin{split}\displaystyle \textbf{Q} =\left[\begin{matrix}\frac{s_{0} \left(T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{3} & \frac{s_{1} \left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right)}{2} & \frac{s_{2} \left(- 2 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{6}\\\frac{s_{0} \left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right)}{3} & \frac{s_{1} \left(T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{2} & \frac{s_{2} \left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right)}{6}\\\frac{s_{0} \left(- 2 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{3} & \frac{s_{1} \left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right)}{2} & \frac{s_{2} \left(4 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{6}\end{matrix}\right].\end{split}\]

By the Eq. (63) can be shifted in time, we recursively over $n$ timestep replace the shifted time function in it:

(64)\[ \underbrace{ \textbf{m}_{k,x}^{t+1} = \textbf{P}\textbf{m}_{k,x}^{t} + \textbf{Q}\textbf{m}_{k,x}^{eq,t} \quad \rightarrow \quad \textbf{m}_{k,x}^{t} = \textbf{P}\textbf{m}_{k,x}^{t-1} + \textbf{Q}\textbf{m}_{k,x}^{eq,t-1} }_{\textrm{Time Shift}} \qquad \textrm{Recursive Substitution $n$ times:} \qquad \textbf{m}_{k,x}^{t+1} = \textbf{P}^{n}\textbf{m}_{k,x}^{t-n+1} + \displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,t-l}. \]

The next steps follow in orde to employ the Calay-Hamilton theorem [7]. In htis way, we need to find the characteristical polinomial $\mathcal{X}_{P}$ of the matrix $\textbf{P}$ (the definition of characteristical polinomial are presented in toogle below), given by:

The characteristic polynomial of a square matrix $(P \in \mathbb{R}^{n \times n})$ (or $(\mathbb{C}^{n \times n}))$ is defined as:

\[ \boxed{\mathcal{X}_{P}(X) = \det(X I - P)} \qquad \rightarrow \qquad \mathcal{X}_{P}(X) = \sum_{i=0}^{n} X^{i}\gamma^{i} \]

where $\gamma_{i}$ are the polynomial coefficients.

\[ \mathcal{X}_{P} = \gamma_{3}X^{3} + \gamma_{2}X^{2} + \gamma_{1}X^{1} + \gamma_{0}, \]

where all the polynomial coefficients $\gamma_{i}$ are shifted by $T_{\Delta x}^{c_{0}}$:

\[\displaystyle \gamma_{3} =T_{\Delta_{x}}^{c_{0}},\]

\[\displaystyle \gamma_{2} =- T_{\Delta_{x}}^{c_{0}} - T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}} + s_{0} \left(\frac{T_{\Delta_{x}}^{c_{0}}}{3} + \frac{T_{\Delta_{x}}^{c_{1}}}{3} + \frac{T_{\Delta_{x}}^{c_{2}}}{3}\right) + s_{1} \left(\frac{T_{\Delta_{x}}^{c_{1}}}{2} + \frac{T_{\Delta_{x}}^{c_{2}}}{2}\right) + s_{2} \left(\frac{2 T_{\Delta_{x}}^{c_{0}}}{3} + \frac{T_{\Delta_{x}}^{c_{1}}}{6} + \frac{T_{\Delta_{x}}^{c_{2}}}{6}\right),\]

\[\displaystyle \gamma_{1} =T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}} + \frac{s_{1} s_{2} \left(T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{3} - \frac{s_{1} \left(2 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{2} - \frac{s_{2} \left(2 T_{\Delta_{x}}^{c_{0}} + 5 T_{\Delta_{x}}^{c_{1}} + 5 T_{\Delta_{x}}^{c_{2}}\right)}{6},\]

\[\displaystyle \gamma_{0} =- T_{\Delta_{x}}^{c_{0}} \left(s_{1} - 1\right) \left(s_{2} - 1\right).\]

To reach the above coeficient results, were considered: $T_{\Delta x}^{c_{1}}T_{\Delta x}^{c_{2}}=T_{\Delta x}^{c_{0}}$; $T_{\Delta x}^{c_{1}}T_{\Delta x}^{c_{0}}=T_{\Delta x}^{c_{1}}$; $T_{\Delta x}^{c_{2}}T_{\Delta x}^{c_{0}}=T_{\Delta x}^{c_{2}}$.Defined the characteristical polynomial, we multiply both sides of Eq. (64) by $\sum_{n=0}^{q}\gamma_{n}$ (where $q$ is the number of lattice velocities) and perform a time shift by change the time variable $t$ to $\tilde{t}+n=t+1$:

\[\begin{split} \begin{array}{c} \textbf{m}_{k,x}^{t+1} = \textbf{P}^{n}\textbf{m}_{k,x}^{t-n+1} + \displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,t-l} \qquad \rightarrow \qquad \textbf{m}_{k,x}^{\tilde{t}+n} = \textbf{P}^{n}\textbf{m}_{k,x}^{\tilde{t}} + \displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,\tilde{t}+n-1-l} \qquad \rightarrow \\ \displaystyle\sum_{n=0}^{q}\gamma_{n} \textbf{m}_{k,x}^{\tilde{t}+n} = \displaystyle\sum_{n=0}^{q}\gamma_{n}\textbf{P}^{n}\textbf{m}_{k,x}^{\tilde{t}} + \displaystyle\sum_{n=0}^{q}\gamma_{n} \left(\displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,\tilde{t}+n-1-l}\right), \end{array} \end{split}\]

by applying the Calay-Hamilton theorem, that imlpie $\sum_{n=0}^{q}\gamma_{n}\textbf{P}^{n}=0$:

\[ \displaystyle\sum_{n=0}^{q}\gamma_{n} \textbf{m}_{k,x}^{\tilde{t}+n} = \underbrace{\left(\displaystyle\sum_{n=0}^{q}\gamma_{n}\textbf{P}^{n}\right)}_{=0}\textbf{m}_{k,x}^{\tilde{t}} + \displaystyle\sum_{n=0}^{q}\gamma_{n}\left( \displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,\tilde{t}+n-1-l} \right) \qquad \rightarrow \qquad \displaystyle\sum_{n=0}^{q}\gamma_{n} \textbf{m}_{k,x}^{\tilde{t}+n} = \displaystyle\sum_{n=0}^{q}\gamma_{n}\left(\displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,\tilde{t}+n-1-l}\right), \]

applying a second time shift by change tha variable back to $\tilde{t}+q=t+1$ and expanding the left hand-side term to isolate the time $t+1$:

\[\begin{split} \begin{array}{c} \displaystyle\sum_{n=0}^{q}\gamma_{n} \textbf{m}_{k,x}^{\tilde{t}+n} = \displaystyle\sum_{n=0}^{q}\gamma_{n}\left(\displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,\tilde{t}+n-1-l}\right) \qquad \rightarrow \qquad \displaystyle\sum_{n=0}^{q}\gamma_{n} \textbf{m}_{k,x}^{t+1-q+n} = \displaystyle\sum_{n=0}^{q}\gamma_{n}\left(\displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,t-q+n-l}\right) \qquad \rightarrow \\ \displaystyle\sum_{n=0}^{q}\gamma_{n} \textbf{m}_{k,x}^{t+1-q+n} = \displaystyle\sum_{n=0}^{q}\gamma_{n}\left(\displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,t-q+n-l}\right) \qquad \rightarrow \qquad \gamma_{q} \textbf{m}_{k,x}^{t+1} = -\displaystyle\sum_{n=0}^{q-1}\gamma_{n} \textbf{m}_{k,x}^{t+1-q+n} + \displaystyle\sum_{n=0}^{q}\gamma_{n}\left(\displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,t-q+n-l}\right), \end{array} \end{split}\]

the last sum term can be rewritten by

\[ \gamma_{q} \textbf{m}_{k,x}^{t+1} = -\displaystyle\sum_{n=0}^{q-1}\gamma_{n} \textbf{m}_{k,x}^{t+1-q+n} + \displaystyle\sum_{n=0}^{q}\gamma_{n}\left(\displaystyle\sum_{l=0}^{n-1}\textbf{P}^{l}\textbf{Q}\textbf{m}_{k,x}^{eq,t-q+n-l}\right) \qquad \rightarrow \qquad \gamma_{q} \textbf{m}_{k,x}^{t+1} = -\displaystyle\sum_{n=0}^{q-1}\gamma_{n} \textbf{m}_{k,x}^{t+1-q+n} + \displaystyle\sum_{n=0}^{q-1}\left(\displaystyle\sum_{l=0}^{n}\gamma_{q+l-n}\textbf{P}^{l}\right)\textbf{Q}\textbf{m}_{k,x}^{eq,t-n}. \]

\[\displaystyle m^{t+1}_{k,x} {\gamma}_{3} = P^{2} Q m^{eq,t-2}_{k,x} {\gamma}_{3} + P Q m^{eq,t-1}_{k,x} {\gamma}_{3} + P Q m^{eq,t-2}_{k,x} {\gamma}_{2} + Q m^{eq,t-1}_{k,x} {\gamma}_{2} + Q m^{eq,t-2}_{k,x} {\gamma}_{1} + Q m^{eq,t}_{k,x} {\gamma}_{3} - m^{t-1}_{k,x} {\gamma}_{1} - m^{t-2}_{k,x} {\gamma}_{0} - m^{t}_{k,x} {\gamma}_{2}\]

Considering the solution for the moment $m^{t+1}_{0,x}$, the above equation is rewritten to:

(65)\[ m^{t+1}_{0,x} \gamma_{3} = \left[P^{2} Q m^{eq,t-2}_{k,x} \gamma_{3}\right]_{k=0} + \left[P Q m^{eq,t-1}_{k,x} \gamma_{3}\right]_{k=0} + \left[P Q m^{eq,t-2}_{k,x} \gamma_{2}\right]_{k=0} + \left[Q m^{eq,t-1}_{k,x} \gamma_{2}\right]_{k=0} + \left[Q m^{eq,t-2}_{k,x} \gamma_{1}\right]_{k=0} + \left[Q m^{eq,t}_{k,x} \gamma_{3}\right]_{k=0} - m^{t-1}_{0,x} \gamma_{1} - m^{t-2}_{0,x} \gamma_{0} - m^{t}_{0,x} \gamma_{2} \]

\[\displaystyle \left[ P Q m^{eq,t-1}_{k,x} \gamma_{3} \right]_{k=0} + \left[ Q m^{eq,t-1}_{k,x} \gamma_{2} \right]_{k=0} =m^{eq,t-1}_{1,x} \left(\frac{T_{\Delta_{x}}^{c_{1}} s_{1} s_{2}}{2} - \frac{T_{\Delta_{x}}^{c_{1}} s_{1}}{2} - \frac{T_{\Delta_{x}}^{c_{2}} s_{1} s_{2}}{2} + \frac{T_{\Delta_{x}}^{c_{2}} s_{1}}{2}\right) + m^{eq,t-1}_{2,x} \left(\frac{T_{\Delta_{x}}^{c_{0}} s_{1} s_{2}}{3} - \frac{T_{\Delta_{x}}^{c_{0}} s_{2}}{3} - \frac{T_{\Delta_{x}}^{c_{1}} s_{1} s_{2}}{6} + \frac{T_{\Delta_{x}}^{c_{1}} s_{2}}{6} - \frac{T_{\Delta_{x}}^{c_{2}} s_{1} s_{2}}{6} + \frac{T_{\Delta_{x}}^{c_{2}} s_{2}}{6}\right),\]

\[\displaystyle \left[ P^{2} Q m^{eq,t-2}_{k,x} \gamma_{3} \right]_{k=0} + \left[ PQ m^{eq,t-2}_{k,x} \gamma_{2} \right]_{k=0} + \left[ Q m^{eq,t-2}_{k,x} \gamma_{1} \right]_{k=0} =0,\]

\[\displaystyle \left[Q m^{eq,t}_{k,x} \gamma_{3} \right]_{k=0} =\frac{m^{eq,t}_{1,x} s_{1} \left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right)}{2} + \frac{m^{eq,t}_{2,x} s_{2} \left(- 2 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{6}.\]

Replacing the terms in the Eq. (65):

\[\displaystyle T_{\Delta_{x}}^{c_{0}} m^{t+1}_{0,x} = T_{\Delta_{x}}^{c_{0}} m^{t-2}_{0,x} \left(s_{1} - 1\right) \left(s_{2} - 1\right) + \frac{m^{eq,t-1}_{1,x} s_{1} \left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right) \left(s_{2} - 1\right)}{2} + \frac{m^{eq,t-1}_{2,x} s_{2} \left(s_{1} - 1\right) \left(2 T_{\Delta_{x}}^{c_{0}} - T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right)}{6} + \frac{m^{eq,t}_{1,x} s_{1} \left(T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}}\right)}{2} + \frac{m^{eq,t}_{2,x} s_{2} \left(- 2 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{6} - m^{t-1}_{0,x} \left(T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}} + \frac{s_{1} s_{2} \left(T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{3} - \frac{s_{1} \left(2 T_{\Delta_{x}}^{c_{0}} + T_{\Delta_{x}}^{c_{1}} + T_{\Delta_{x}}^{c_{2}}\right)}{2} - \frac{s_{2} \left(2 T_{\Delta_{x}}^{c_{0}} + 5 T_{\Delta_{x}}^{c_{1}} + 5 T_{\Delta_{x}}^{c_{2}}\right)}{6}\right) - m^{t}_{0,x} \left(- T_{\Delta_{x}}^{c_{0}} - T_{\Delta_{x}}^{c_{1}} - T_{\Delta_{x}}^{c_{2}} + s_{1} \left(\frac{T_{\Delta_{x}}^{c_{1}}}{2} + \frac{T_{\Delta_{x}}^{c_{2}}}{2}\right) + s_{2} \left(\frac{2 T_{\Delta_{x}}^{c_{0}}}{3} + \frac{T_{\Delta_{x}}^{c_{1}}}{6} + \frac{T_{\Delta_{x}}^{c_{2}}}{6}\right)\right)\]

\[\displaystyle m^{t+1}_{0,x} = m^{eq,t-1}_{1,x+1} \left(- \frac{s_{1} s_{2}}{2} + \frac{s_{1}}{2}\right) + m^{eq,t-1}_{1,x-1} \left(\frac{s_{1} s_{2}}{2} - \frac{s_{1}}{2}\right) + m^{eq,t-1}_{2,x+1} \left(- \frac{s_{1} s_{2}}{6} + \frac{s_{2}}{6}\right) + m^{eq,t-1}_{2,x-1} \left(- \frac{s_{1} s_{2}}{6} + \frac{s_{2}}{6}\right) + m^{eq,t-1}_{2,x} \left(\frac{s_{1} s_{2}}{3} - \frac{s_{2}}{3}\right) - \frac{m^{eq,t}_{1,x+1} s_{1}}{2} + \frac{m^{eq,t}_{1,x-1} s_{1}}{2} + \frac{m^{eq,t}_{2,x+1} s_{2}}{6} + \frac{m^{eq,t}_{2,x-1} s_{2}}{6} - \frac{m^{eq,t}_{2,x} s_{2}}{3} + m^{t-1}_{0,x+1} \left(- \frac{s_{1} s_{2}}{3} + \frac{s_{1}}{2} + \frac{5 s_{2}}{6} - 1\right) + m^{t-1}_{0,x-1} \left(- \frac{s_{1} s_{2}}{3} + \frac{s_{1}}{2} + \frac{5 s_{2}}{6} - 1\right) + m^{t-1}_{0,x} \left(- \frac{s_{1} s_{2}}{3} + s_{1} + \frac{s_{2}}{3} - 1\right) + m^{t-2}_{0,x} \left(s_{1} - 1\right) \left(s_{2} - 1\right) + m^{t}_{0,x+1} \left(- \frac{s_{1}}{2} - \frac{s_{2}}{6} + 1\right) + m^{t}_{0,x-1} \left(- \frac{s_{1}}{2} - \frac{s_{2}}{6} + 1\right) + m^{t}_{0,x} \left(1 - \frac{2 s_{2}}{3}\right)\]

Replacing the moments:

\[ m_{0,x}^{t}=u_{x}^{t}, \qquad m_{0,x}^{eq,t}=u_{x}^{t}, \qquad m_{1,x}^{eq,t}=B_{x}^{t}, \qquad \textrm{and} \qquad m_{1,x}^{eq,t}=C_{x}^{t} + 3\Lambda_{x}^{t} - 2u_{x}^{t}, \]

in the above equation:

\[\displaystyle u^{t+1}_{x} = B^{t-1}_{x+1} \left(- \frac{s_{1} s_{2}}{2} + \frac{s_{1}}{2}\right) + B^{t-1}_{x-1} \left(\frac{s_{1} s_{2}}{2} - \frac{s_{1}}{2}\right) - \frac{B^{t}_{x+1} s_{1}}{2} + \frac{B^{t}_{x-1} s_{1}}{2} + \frac{s_{2} \left(C^{t}_{x+1} + 3 \Lambda^{t}_{x+1} - 2 u^{t}_{x+1}\right)}{6} + \frac{s_{2} \left(C^{t}_{x-1} + 3 \Lambda^{t}_{x-1} - 2 u^{t}_{x-1}\right)}{6} - \frac{s_{2} \left(C^{t}_{x} + 3 \Lambda^{t}_{x} - 2 u^{t}_{x}\right)}{3} + u^{t-1}_{x+1} \left(- \frac{s_{1} s_{2}}{3} + \frac{s_{1}}{2} + \frac{5 s_{2}}{6} - 1\right) + u^{t-1}_{x-1} \left(- \frac{s_{1} s_{2}}{3} + \frac{s_{1}}{2} + \frac{5 s_{2}}{6} - 1\right) + u^{t-1}_{x} \left(- \frac{s_{1} s_{2}}{3} + s_{1} + \frac{s_{2}}{3} - 1\right) + u^{t-2}_{x} \left(s_{1} - 1\right) \left(s_{2} - 1\right) + u^{t}_{x+1} \left(- \frac{s_{1}}{2} - \frac{s_{2}}{6} + 1\right) + u^{t}_{x-1} \left(- \frac{s_{1}}{2} - \frac{s_{2}}{6} + 1\right) + u^{t}_{x} \left(1 - \frac{2 s_{2}}{3}\right) + \left(- \frac{s_{1} s_{2}}{6} + \frac{s_{2}}{6}\right) \left(C^{t-1}_{x+1} + 3 \Lambda^{t-1}_{x+1} - 2 u^{t-1}_{x+1}\right) + \left(- \frac{s_{1} s_{2}}{6} + \frac{s_{2}}{6}\right) \left(C^{t-1}_{x-1} + 3 \Lambda^{t-1}_{x-1} - 2 u^{t-1}_{x-1}\right) + \left(\frac{s_{1} s_{2}}{3} - \frac{s_{2}}{3}\right) \left(C^{t-1}_{x} + 3 \Lambda^{t-1}_{x} - 2 u^{t-1}_{x}\right)\]

rewriting the above equation, we have:

\[\displaystyle u^{t+1}_{x} = - B^{t-1}_{x+1} \eta_{2} + B^{t-1}_{x-1} \eta_{2} - B^{t}_{x+1} \eta_{1} + B^{t}_{x-1} \eta_{1} + \alpha_{1} u^{t}_{x} + \alpha_{2} u^{t}_{x+1} + \alpha_{2} u^{t}_{x-1} + \beta_{1} u^{t-1}_{x} + \beta_{2} u^{t-1}_{x+1} + \beta_{2} u^{t-1}_{x-1} + \gamma u^{t-2}_{x} + \kappa_{1} \left(C^{t}_{x+1} + 3 \Lambda^{t}_{x+1} - 2 u^{t}_{x+1}\right) + \kappa_{1} \left(C^{t}_{x-1} + 3 \Lambda^{t}_{x-1} - 2 u^{t}_{x-1}\right) - 2 \kappa_{1} \left(C^{t}_{x} + 3 \Lambda^{t}_{x} - 2 u^{t}_{x}\right) + \kappa_{2} \left(C^{t-1}_{x+1} + 3 \Lambda^{t-1}_{x+1} - 2 u^{t-1}_{x+1}\right) + \kappa_{2} \left(C^{t-1}_{x-1} + 3 \Lambda^{t-1}_{x-1} - 2 u^{t-1}_{x-1}\right) - 2 \kappa_{2} \left(C^{t-1}_{x} + 3 \Lambda^{t-1}_{x} - 2 u^{t-1}_{x}\right)\]

where coefficients are given by:

\[\displaystyle \alpha_{1} =1 - \frac{2 s_{2}}{3},\qquad \alpha_{2}=- \frac{s_{1}}{2} - \frac{s_{2}}{6} + 1,\qquad \beta_{1}=- \frac{s_{1} s_{2}}{3} + s_{1} + \frac{s_{2}}{3} - 1,\qquad \beta_{2}=- \frac{s_{1} s_{2}}{3} + \frac{s_{1}}{2} + \frac{5 s_{2}}{6} - 1\]

\[\displaystyle \eta_{1}=\frac{s_{1}}{2},\qquad \eta_{2}=\frac{s_{1} s_{2}}{2} - \frac{s_{1}}{2},\qquad \kappa_{1}=\frac{s_{2}}{6},\qquad \kappa_{1}=- \frac{s_{1} s_{2}}{6} + \frac{s_{2}}{6},\qquad \gamma=\left(s_{1} - 1\right) \left(s_{2} - 1\right)\]

Linear Convective-Diffusive Equation

Considering linear convective-diffusive equation:

\[ \partial_t u + \partial_{\alpha}c_{\alpha}u = \partial_{\alpha}\left(\nu \partial_{\alpha} u \right) \]

we have $B_{\alpha}=c_{\alpha}u$, $c_{\alpha}$ is a constant velocity, and $C=u$. Replacing $B_{\alpha}$ and $C$ in FD scheme obtained

\[\displaystyle u^{t+1}_{x} = u^{t-1}_{x+1} \left(\overline{c}_{\alpha} \left(- \frac{s_{1} s_{2}}{2} + \frac{s_{1}}{2}\right) - \frac{s_{1} s_{2}}{3} + \frac{s_{1}}{2} + \frac{5 s_{2}}{6} + \left(3 \overline{c}_{\alpha}^{2} - 1\right) \left(- \frac{s_{1} s_{2}}{6} + \frac{s_{2}}{6}\right) - 1\right) + u^{t-1}_{x-1} \left(\overline{c}_{\alpha} \left(\frac{s_{1} s_{2}}{2} - \frac{s_{1}}{2}\right) - \frac{s_{1} s_{2}}{3} + \frac{s_{1}}{2} + \frac{5 s_{2}}{6} + \left(3 \overline{c}_{\alpha}^{2} - 1\right) \left(- \frac{s_{1} s_{2}}{6} + \frac{s_{2}}{6}\right) - 1\right) + u^{t-1}_{x} \left(- \frac{s_{1} s_{2}}{3} + s_{1} + \frac{s_{2}}{3} + \left(3 \overline{c}_{\alpha}^{2} - 1\right) \left(\frac{s_{1} s_{2}}{3} - \frac{s_{2}}{3}\right) - 1\right) + u^{t-2}_{x} \left(s_{1} - 1\right) \left(s_{2} - 1\right) + u^{t}_{x+1} \left(- \frac{\overline{c}_{\alpha} s_{1}}{2} - \frac{s_{1}}{2} + \frac{s_{2} \left(3 \overline{c}_{\alpha}^{2} - 1\right)}{6} - \frac{s_{2}}{6} + 1\right) + u^{t}_{x-1} \left(\frac{\overline{c}_{\alpha} s_{1}}{2} - \frac{s_{1}}{2} + \frac{s_{2} \left(3 \overline{c}_{\alpha}^{2} - 1\right)}{6} - \frac{s_{2}}{6} + 1\right) + u^{t}_{x} \left(- \frac{s_{2} \left(3 \overline{c}_{\alpha}^{2} - 1\right)}{3} - \frac{2 s_{2}}{3} + 1\right)\]

where $\overline{c}_{\alpha} = c_{\alpha}\Delta x/\Delta t$ is the constant velocity in the lattice Boltzman scale. Rewriting the above equation with coeficients, we have:

\[\displaystyle u^{t+1}_{x} = \alpha_{1} u^{t}_{x} + \alpha_{2} u^{t}_{x-1} + \alpha_{3} u^{t}_{x+1} + \beta_{1} u^{t-1}_{x} + \beta_{2} u^{t-1}_{x-1} + \beta_{3} u^{t-1}_{x+1} + \gamma u^{t-2}_{x}\]

where coefficients are given by:

\[\displaystyle \alpha_{1} =- \frac{s_{2} \left(3 \overline{c}_{\alpha}^{2} - 1\right)}{3} - \frac{2 s_{2}}{3} + 1,\qquad \alpha_{2}=\frac{\overline{c}_{\alpha} s_{1}}{2} - \frac{s_{1}}{2} + \frac{s_{2} \left(3 \overline{c}_{\alpha}^{2} - 1\right)}{6} - \frac{s_{2}}{6} + 1,\qquad \alpha_{3}=- \frac{\overline{c}_{\alpha} s_{1}}{2} - \frac{s_{1}}{2} + \frac{s_{2} \left(3 \overline{c}_{\alpha}^{2} - 1\right)}{6} - \frac{s_{2}}{6} + 1\]

\[\displaystyle \beta_{1}=- \frac{s_{1} s_{2}}{3} + s_{1} + \frac{s_{2}}{3} + \left(3 \overline{c}_{\alpha}^{2} - 1\right) \left(\frac{s_{1} s_{2}}{3} - \frac{s_{2}}{3}\right) - 1,\qquad\beta_{2}=\overline{c}_{\alpha} \left(\frac{s_{1} s_{2}}{2} - \frac{s_{1}}{2}\right) - \frac{s_{1} s_{2}}{3} + \frac{s_{1}}{2} + \frac{5 s_{2}}{6} + \left(3 \overline{c}_{\alpha}^{2} - 1\right) \left(- \frac{s_{1} s_{2}}{6} + \frac{s_{2}}{6}\right) - 1\]

\[\displaystyle \beta_{3}=\overline{c}_{\alpha} \left(- \frac{s_{1} s_{2}}{2} + \frac{s_{1}}{2}\right) - \frac{s_{1} s_{2}}{3} + \frac{s_{1}}{2} + \frac{5 s_{2}}{6} + \left(3 \overline{c}_{\alpha}^{2} - 1\right) \left(- \frac{s_{1} s_{2}}{6} + \frac{s_{2}}{6}\right) - 1,\qquad \gamma=\left(s_{1} - 1\right) \left(s_{2} - 1\right)\]

Fourth-order Convective-Diffusive Equation

Expanding in Taylor series the FD form LB equation, we can truncate the equation with errors isolated in third and fourth-order of scale of $\Delta x$:

\[\begin{split} \begin{array}{c} 0 = \displaystyle\frac{\Delta_{t}^{4} \left(- 15 s_{1} s_{2} + 14 s_{1} + 14 s_{2} - 12\right)}{24} \frac{\partial^{4}{u}}{\partial{t^{4}}} - \frac{\Delta_{t}^{4} c_{\alpha} s_{1} \left(s_{2} - 1\right)}{6}\frac{\partial^{4}{u}}{\partial{x}\partial{t^{3}}} + \frac{\Delta_{t}^{3} \left(7 s_{1} s_{2} - 6 s_{1} - 6 s_{2} + 6\right)}{6}\frac{\partial^{3}{u}}{\partial{t^{3}}} + \frac{\Delta_{t}^{3} c_{\alpha} s_{1} \left(s_{2} - 1\right)}{2} \frac{\partial^{3}{u}}{\partial{x}\partial{t^{2}}} \\ + \displaystyle\frac{\Delta_{t}^{2} \Delta_{x}^{2} \left(\frac{3 \Delta_{t}^{2} c_{\alpha}^{2} s_{1} s_{2}}{\Delta_{x}^{2}} - \frac{3 \Delta_{t}^{2} c_{\alpha}^{2} s_{2}}{\Delta_{x}^{2}} + s_{1} s_{2} - 3 s_{1} - 4 s_{2} + 6\right)}{12}\frac{\partial^{4}{u}}{\partial{x^{2}}\partial{t^{2}}} - \frac{\Delta_{t}^{2} \Delta_{x}^{2} c_{\alpha} s_{1} \left(s_{2} - 1\right)}{6}\frac{\partial^{4}{u}}{\partial{x^{3}}\partial{t}} + \frac{\Delta_{t}^{2} \left(- 3 s_{1} s_{2} + 2 s_{1} + 2 s_{2}\right)}{2} \frac{\partial^{2}{u}}{\partial{t^{2}}} - \Delta_{t}^{2} \frac{\partial^{2}{u}}{\partial{x}\partial{t}} c_{\alpha} s_{1} \left(s_{2} - 1\right) \\ - \displaystyle\frac{\Delta_{t} \Delta_{x}^{2} \left(\frac{3 \Delta_{t}^{2} c_{\alpha}^{2} s_{1} s_{2}}{\Delta_{x}^{2}} - \frac{3 \Delta_{t}^{2} c_{\alpha}^{2} s_{2}}{\Delta_{x}^{2}} + s_{1} s_{2} - 3 s_{1} - 4 s_{2} + 6\right)}{6} \frac{\partial^{3}{u}}{\partial{x^{2}}\partial{t}} + \frac{\Delta_{t} \Delta_{x}^{2} c_{\alpha} s_{1} s_{2}}{6} \frac{\partial^{3}{u}}{\partial{x^{3}}} + \Delta_{t} \frac{\partial{u}}{\partial{t}} s_{1} s_{2} + \Delta_{t} \frac{\partial{u}}{\partial{x}} c_{\alpha} s_{1} s_{2} + \frac{\Delta_{x}^{4} s_{2} \left(\frac{3 \Delta_{t}^{2} c_{\alpha}^{2} s_{1}}{\Delta_{x}^{2}} - \frac{6 \Delta_{t}^{2} c_{\alpha}^{2}}{\Delta_{x}^{2}} + s_{1} - 2\right)}{72} \frac{\partial^{4}{u}}{\partial{x^{4}}} \\ + \displaystyle\frac{\Delta_{x}^{2} s_{2} \left(\frac{3 \Delta_{t}^{2} c_{\alpha}^{2} s_{1}}{\Delta_{x}^{2}} - \frac{6 \Delta_{t}^{2} c_{\alpha}^{2}}{\Delta_{x}^{2}} + s_{1} - 2\right)}{6} \frac{\partial^{2}{u}}{\partial{x^{2}}} , \end{array} \end{split}\]

isolating $\frac{\partial u}{\partial t}$:

(66)\[\begin{split} \begin{array}{c} \displaystyle\frac{\partial u}{\partial t} + \frac{\partial{u}}{\partial{x}} c_{\alpha} = \frac{\partial^{2}{u}}{\partial{x^{2}}} \nu + \Delta_{t}^{3} \frac{\partial^{4}{u}}{\partial{t^{4}}} \left(\frac{5}{8} - \frac{7}{12 s_{2}} - \frac{7}{12 s_{1}} + \frac{1}{2 s_{1} s_{2}}\right) - \frac{\Delta_{t}^{3} c_{\alpha}^{2}}{4}\frac{\partial^{4}{u}}{\partial{x^{2}}\partial{t^{2}}} + \frac{\Delta_{t}^{3} c_{\alpha}^{2}}{4 s_{1}}\frac{\partial^{4}{u}}{\partial{x^{2}}\partial{t^{2}}} \\ + \Delta_{t}^{3} \displaystyle\frac{\partial^{4}{u}}{\partial{x}\partial{t^{3}}} c_{\alpha} \left(\frac{1}{6} - \frac{1}{6 s_{2}}\right) + \Delta_{t}^{2} \frac{\partial^{3}{u}}{\partial{t^{3}}} \left(- \frac{7}{6} + \frac{1}{s_{2}} + \frac{1}{s_{1}} - \frac{1}{s_{1} s_{2}}\right) + \frac{\Delta_{t}^{2} c_{\alpha}^{2}}{2} \frac{\partial^{3}{u}}{\partial{x^{2}}\partial{t}} - \frac{\Delta_{t}^{2} c_{\alpha}^{2}}{2 s_{1}} \frac{\partial^{3}{u}}{\partial{x^{2}}\partial{t}} + \Delta_{t}^{2} \frac{\partial^{3}{u}}{\partial{x}\partial{t^{2}}} c_{\alpha} \left(- \frac{1}{2} + \frac{1}{2 s_{2}}\right) \\ + \Delta_{t} \Delta_{x}^{2} \displaystyle\frac{\partial^{4}{u}}{\partial{x^{2}}\partial{t^{2}}} \left(- \frac{1}{12} + \frac{1}{4 s_{2}} + \frac{1}{3 s_{1}} - \frac{1}{2 s_{1} s_{2}}\right) + \Delta_{t} \Delta_{x}^{2} \frac{\partial^{4}{u}}{\partial{x^{3}}\partial{t}} c_{\alpha} \left(\frac{1}{6} - \frac{1}{6 s_{2}}\right) - \frac{\Delta_{t} \Delta_{x}^{2} c_{\alpha}^{2}}{24} \frac{\partial^{4}{u}}{\partial{x^{4}}} + \frac{\Delta_{t} \Delta_{x}^{2} c_{\alpha}^{2}}{12 s_{1}} \frac{\partial^{4}{u}}{\partial{x^{4}}} \\ + \Delta_{t} \displaystyle\frac{\partial^{2}{u}}{\partial{t^{2}}} \left(\frac{3}{2} - \frac{1}{s_{2}} - \frac{1}{s_{1}}\right) - \frac{\Delta_{t} c_{\alpha}^{2}}{2} \frac{\partial^{2}{u}}{\partial{x^{2}}} + \frac{\Delta_{t} c_{\alpha}^{2}}{s_{1}} \frac{\partial^{2}{u}}{\partial{x^{2}}} + \Delta_{t} \frac{\partial^{2}{u}}{\partial{x}\partial{t}} c_{\alpha} \left(1 - \frac{1}{s_{2}}\right) + \Delta_{x}^{2} \frac{\partial^{3}{u}}{\partial{x^{2}}\partial{t}} \left(\frac{1}{6} - \frac{1}{2 s_{2}} - \frac{2}{3 s_{1}} + \frac{1}{s_{1} s_{2}}\right) - \frac{\Delta_{x}^{2} c_{\alpha}}{6}\frac{\partial^{3}{u}}{\partial{x^{3}}} + \frac{\Delta_{x}^{4} \left(- \frac{1}{72} + \frac{1}{36 s_{1}}\right)}{\Delta_{t}}\frac{\partial^{4}{u}}{\partial{x^{4}}} \end{array} \end{split}\]

where $\nu=c_{s}^{2}\left(\displaystyle\frac{2-s1}{2s1}\right)\displaystyle\frac{\Delta x^{2}}{\Delta t}$

To correct reduce the above equation and isolate the third and fourth order of spatial derivative terms, we employ the relations:

\[ \displaystyle\frac{\partial^{2}u}{\partial x\,\partial t} = \nu \displaystyle\frac{\partial^{3}u}{\partial x^{3}} - c_{\alpha} \frac{\partial^{2}u}{\partial x^{2}} + \mathcal{O}\!\left(\Delta^{3}x + \Delta t\,\Delta x\right), \]

\[ \displaystyle\frac{\partial^{2}u}{\partial t^{2}} = \nu^{2} \displaystyle\frac{\partial^{4}u}{\partial x^{4}} - 2\nu c_{\alpha} \frac{\partial^{3}u}{\partial x^{3}} + c_{\alpha}^{2} \frac{\partial^{2}u}{\partial x^{2}} + \mathcal{O}(\Delta t), \]

\[ \displaystyle\frac{\partial^{3}u}{\partial x\,\partial t^{2}} = -2\nu c_{\alpha} \displaystyle \frac{\partial^{4}u}{\partial x^{4}} + c_{\alpha}^{2} \frac{\partial^{3}u}{\partial x^{3}} + \mathcal{O}(\Delta x), \]

\[ \displaystyle\frac{\partial^{3}u}{\partial x^{2}\partial t} = \nu \displaystyle\frac{\partial^{4}u}{\partial x^{4}} - c_{\alpha} \frac{\partial^{3}u}{\partial x^{3}} + \mathcal{O}\!\left(\frac{\Delta^{4}x}{\Delta t} + \Delta^{2}x\right), \]

\[ \displaystyle\frac{\partial^{3}u}{\partial t^{3}} = 3\nu c_{\alpha}^{2} \displaystyle\frac{\partial^{4}u}{\partial x^{4}} - c_{\alpha}^{3} \frac{\partial^{3}u}{\partial x^{3}} + \mathcal{O}\!\left(\frac{\Delta^{2}x}{\Delta t}\right), \]

\[ \displaystyle\frac{\partial^{4}u}{\partial x\,\partial t^{3}} = - c_{\alpha}^{3} \displaystyle\frac{\partial^{4}u}{\partial x^{4}} + \mathcal{O}\!\left(\frac{1}{\Delta x}\right), \]

\[ \displaystyle\frac{\partial^{4}u}{\partial x^{3}\partial t} = - c_{\alpha} \displaystyle\frac{\partial^{4}u}{\partial x^{4}} + \mathcal{O}(\Delta x), \]

\[ \displaystyle\frac{\partial^{4}u}{\partial x^{2}\partial t^{2}} = c_{\alpha}^{2} \displaystyle\frac{\partial^{4}u}{\partial x^{4}} + \mathcal{O}\!\left(\frac{\Delta^{2}x}{\Delta t}\right), \]

\[ \displaystyle\frac{\partial^{4}u}{\partial t^{4}} = c_{\alpha}^{4} \displaystyle\frac{\partial^{4}u}{\partial x^{4}} + \mathcal{O}\!\left(\frac{1}{\Delta t}\right), \]

notice that in the above relations, the derivativer higher or equat to fifth order are neglected. Replacing them in Eq. (66):

\[\begin{split} \begin{array}{c} \displaystyle\frac{\partial u}{\partial t} + \frac{\partial{u}}{\partial{x}} c_{\alpha} = \frac{\partial^{2}{u}}{\partial{x^{2}}} \nu + \frac{\partial^{3}{u}}{\partial{x^{3}}} \left(\frac{\Delta_{t}^{2} c_{\alpha}^{3}}{6} - \frac{\Delta_{t}^{2} c_{\alpha}^{3}}{2 s_{2}} - \frac{\Delta_{t}^{2} c_{\alpha}^{3}}{2 s_{1}} + \frac{\Delta_{t}^{2} c_{\alpha}^{3}}{s_{1} s_{2}} - 2 \Delta_{t} \nu c_{\alpha} + \frac{\Delta_{t} \nu c_{\alpha}}{s_{2}} + \frac{2 \Delta_{t} \nu c_{\alpha}}{s_{1}} - \frac{\Delta_{x}^{2} c_{\alpha}}{3} + \frac{\Delta_{x}^{2} c_{\alpha}}{2 s_{2}} + \frac{2 \Delta_{x}^{2} c_{\alpha}}{3 s_{1}} - \frac{\Delta_{x}^{2} c_{\alpha}}{s_{1} s_{2}}\right) + \\ \displaystyle\frac{\partial^{4}{u}}{\partial{x^{4}}} \left(\frac{5 \Delta_{t}^{3} c_{\alpha}^{4}}{24} - \frac{5 \Delta_{t}^{3} c_{\alpha}^{4}}{12 s_{2}} - \frac{\Delta_{t}^{3} c_{\alpha}^{4}}{3 s_{1}} + \frac{\Delta_{t}^{3} c_{\alpha}^{4}}{2 s_{1} s_{2}} - 2 \Delta_{t}^{2} \nu c_{\alpha}^{2} + \frac{2 \Delta_{t}^{2} \nu c_{\alpha}^{2}}{s_{2}} + \frac{5 \Delta_{t}^{2} \nu c_{\alpha}^{2}}{2 s_{1}} - \frac{3 \Delta_{t}^{2} \nu c_{\alpha}^{2}}{s_{1} s_{2}} - \frac{7 \Delta_{t} \Delta_{x}^{2} c_{\alpha}^{2}}{24} + \frac{5 \Delta_{t} \Delta_{x}^{2} c_{\alpha}^{2}}{12 s_{2}} + \frac{5 \Delta_{t} \Delta_{x}^{2} c_{\alpha}^{2}}{12 s_{1}} - \frac{\Delta_{t} \Delta_{x}^{2} c_{\alpha}^{2}}{2 s_{1} s_{2}} - \frac{\Delta_{x}^{2} \nu}{12} - \frac{\Delta_{x}^{2} \nu}{3 s_{2}} + \frac{2 \Delta_{x}^{2} \nu}{3 s_{1} s_{2}} - \frac{\Delta_{x}^{2} \nu}{3 s_{1}^{2}} - \frac{\Delta_{x}^{4}}{72 \Delta_{t}} + \frac{\Delta_{x}^{4}}{36 \Delta_{t} s_{1}}\right) \end{array} \end{split}\]

\[ \begin{array}{c} \displaystyle\frac{\partial u}{\partial t} + \frac{\partial{u}}{\partial{x}} c_{\alpha} = \frac{\partial^{2}{u}}{\partial{x^{2}}} \nu + \frac{c_{\alpha}\Delta x^{2}TR_{3}}{6 s_{1}^{2}s_{2}}\frac{\partial^{3}{u}}{\partial{x^{3}}} + \frac{\Delta x^{4}TR_{4}}{72\Delta t s_{1}^{3}s_{2}} \displaystyle\frac{\partial^{4}{u}}{\partial{x^{4}}} \end{array} \]

where $TR_{3}$ and $TR_{4}$ are given by:

\[ TR_{3}=\overline{c}_{\alpha}^{2} s_{1} \left(s_{1} s_{2} - 3 s_{1} - 3 s_{2} + 6\right) + 2 \left(s_{1} - 2\right) \left(s_{1} - s_{2}\right) \]

\[ TR_{4}=\overline{c}_{\alpha}^{4} \left(15 s_{1}^{3} s_{2} - 30 s_{1}^{3} - 24 s_{1}^{2} s_{2} + 36 s_{1}^{2}\right) + \overline{c}_{\alpha}^{2} \left(3 s_{1}^{3} s_{2} + 6 s_{1}^{3} - 48 s_{1}^{2} s_{2} + 48 s_{1}^{2} + 60 s_{1} s_{2} - 72 s_{1}\right) + 4 s_{1}^{3} - 16 s_{1}^{2} + 4 s_{1} s_{2} + 16 s_{1} - 8 s_{2} \]

Fourth-Order Convergence

Making the term $TR_{3}$ and $TR_{4}$ equal zero and solving the equation, we have:

\[ TR_{3}:\qquad s_{2}=\frac{s_{1} \left(3 \bar{c}_{\alpha}^{2} s_{1} - 6 \bar{c}_{\alpha}^{2} - 2 s_{1} + 4\right)}{\bar{c}_{\alpha}^{2} s_{1}^{2} - 3 \bar{c}_{\alpha}^{2} s_{1} - 2 s_{1} + 4}, \]

\[ TR_{4}:\qquad s_{2}=\frac{2 s_{1} \left(15 \bar{c}^{4} s_{1}^{2} - 18 \bar{c}^{4} s_{1} - 3 \bar{c}^{2} s_{1}^{2} - 24 \bar{c}^{2} s_{1} + 36 \bar{c}^{2} - 2 s_{1}^{2} + 8 s_{1} - 8\right)}{15 \bar{c}^{4} s_{1}^{3} - 24 \bar{c}^{4} s_{1}^{2} + 3 \bar{c}^{2} s_{1}^{3} - 48 \bar{c}^{2} s_{1}^{2} + 60 \bar{c}^{2} s_{1} + 4 s_{1} - 8} \]

Plotting both solution for different values of $\bar{c}_{\alpha}$:

Show code cell source Hide code cell source

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'cm'
# ---------------- Symbols ----------------
s1, s2, cbar = sp.symbols('s1 s2 cbar')
A = cbar**2

# ---------------- Eq 1: solve for s2 ----------------
s2_expr_1 = sp.solve(sp.Eq(A*s1*(s1*s2 - 3*s1 - 3*s2 + 6) + 2*(s1-2)*(s1-s2), 0), s2)[0]

# ---------------- Eq 2: s2(s1) (your big equation) ----------------
s2_expr_2 = (2*s1*(15*cbar**4*s1**2 - 18*cbar**4*s1
          - 3*cbar**2*s1**2 - 24*cbar**2*s1
          + 36*cbar**2 - 2*s1**2 + 8*s1 - 8)
) / (15*cbar**4*s1**3 - 24*cbar**4*s1**2
    + 3*cbar**2*s1**3 - 48*cbar**2*s1**2
    + 60*cbar**2*s1 + 4*s1 - 8)

# difference g(s1) = f1 - f2  (roots => intersections)
g_expr = sp.simplify(s2_expr_1 - s2_expr_2)

# lambdify for fast evaluation
f1 = sp.lambdify((s1, cbar), s2_expr_1, "numpy")
f2 = sp.lambdify((s1, cbar), s2_expr_2, "numpy")
g  = sp.lambdify((s1, cbar), g_expr,  "numpy")

# numeric root refinement (mpmath via sympy)
def refine_root(cb, x0):
    x = sp.nsolve(g_expr.subs(cbar, cb), s1, x0)
    return float(x)

def unique_sorted(xs, tol=1e-6):
    xs = sorted(xs)
    out = []
    for x in xs:
        if not out or abs(x - out[-1]) > tol:
            out.append(x)
    return out

# ---------------- Parameters ----------------
cbar_vals = [0.005, 0.01, 0.05, 0.1]
s1_min, s1_max = 0.01, 3.0
Nscan = 4000
s1_grid = np.linspace(s1_min, s1_max, Nscan)

# ---------------- Plot ----------------
plt.figure(figsize=(9, 5))

for cb in cbar_vals:
    y1 = f1(s1_grid, cb)
    y2 = f2(s1_grid, cb)
    yg = g(s1_grid, cb)

    # mask blow-ups / poles
    y1 = np.where(np.isfinite(y1) & (np.abs(y1) < 10), y1, np.nan)
    y2 = np.where(np.isfinite(y2) & (np.abs(y2) < 10), y2, np.nan)
    yg = np.where(np.isfinite(yg) & (np.abs(yg) < 1e6), yg, np.nan)

    # draw both curves
    plt.plot(s1_grid, y1, linestyle='-',  label=rf"$TR_{3}$, $\bar c_\alpha={cb}$")
    plt.plot(s1_grid, y2, linestyle='--', label=rf"$TR_{4}$, $\bar c_\alpha={cb}$")

    # --- find sign changes of g on the scan grid ---
    roots_guess = []
    for i in range(len(s1_grid) - 1):
        a, b = yg[i], yg[i+1]
        if np.isnan(a) or np.isnan(b):
            continue
        if a == 0.0:
            roots_guess.append(s1_grid[i])
        elif a * b < 0:  # sign change
            roots_guess.append(0.5 * (s1_grid[i] + s1_grid[i+1]))

    # refine guesses with nsolve
    roots = []
    for x0 in roots_guess:
        try:
            xr = refine_root(cb, x0)
            if s1_min <= xr <= s1_max:
                roots.append(xr)
        except Exception:
            pass

    roots = unique_sorted(roots, tol=1e-5)

    # compute s2 at intersection points and plot markers
    for xr in roots:
        yr = float(sp.N(s2_expr_1.subs({cbar: cb, s1: xr})))
        if np.isfinite(yr) and abs(yr) < 10:
            plt.plot([xr], [yr], marker='o', linestyle='None')

    if roots:
        print(f"cbar={cb}: intersections at s1 =", roots)

plt.xlabel(r"$s_1$",fontsize=16)
plt.ylabel(r"$s_2$",fontsize=16)
plt.xlim(0, 2)
plt.ylim(0, 3)
plt.grid(True)
plt.legend(ncol=2, fontsize=9)
plt.tight_layout()
plt.show()

cbar=0.005: intersections at s1 = [0.9999937477340344]
cbar=0.01: intersections at s1 = [0.9999749637281587]
cbar=0.05: intersections at s1 = [0.9993519974731944]
cbar=0.1: intersections at s1 = [0.9971142700951519]

../../_images/46faa380191cd7a01d2f0b4046e9d5b9030f895e7651eb393a7bcbdae3c44880.png

Show code cell source Hide code cell source

import warnings
warnings.filterwarnings("ignore")
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'cm'
matplotlib.rcParams['mathtext.fontset'] = 'cm'
#--------------------------------------- Parameters (physical units) ------------------------------------------
A   = 1.0                         # Amplitude
phi = 0.0                         # Wave initial shift
k_w   = 1.0                       # Wave number
c   = 0.0125                         # wave speed
nu  = 0.066666666666666666*np.pi # Diffusive coefficient
m_waves = 1                       # Wave Period
L = m_waves * (2.0*np.pi / k_w)   # Domain Lentgh
t_end = 24.0                       # Final time
#--------------------------------------- Lattice-Properties-D1Q3 ----------------------------------------------
w = np.array([4.0/6.0, 1.0/6.0, 1.0/6.0],dtype="float64")
cx = np.array([0, 1, -1],dtype="int16")  
cs=1.0/np.sqrt(3.0);
#--------------------------------- Initialization - Save Data Arrays -----------------------------------------------
cases=4
Nx0 = np.array([10, 20, 40, 80],dtype="int64")
r0 = 2**(0)*np.array([2.0**(1), 2.0**(2), 2.0**(3), 2.0**(4)],dtype="float64")
u_cases = np.empty(cases, dtype=object) # array field used to same data over time
for i in range(cases):
    u_cases[i] = np.zeros((Nx0[i]), dtype="float64")
for case in range(0,cases):
    #------------------------------------- Parameters (numerical units) -------------------------------------------
    Nx=Nx0[case]  # Numerical Length
    x  = np.linspace(0.0, L, Nx, endpoint=False,dtype="float64") # Numerical Domain
    dx = L / Nx # Grid size
    r = 1*r0[case] # dx/dt relation
    dt = dx/r
    ce = c/r
    nue = nu * dt / (dx*dx)
    tau1 = 0.5 + nue / cs**2
    om1=1.0/tau1
    om2=om1*(3.0*ce*ce*om1 - 6.0*ce*ce - 2.0*om1 +4.0)/(ce*ce*om1*om1 - 3.0*ce*ce*om1 - 2.0*om1 +4.0)
    tau2=1.0/om2
    nt = int(np.ceil(t_end / dt))
    print(f"dx={dx:.4e},\t dt={dt:.4e}\t nt={nt:d},\t tn={nt*dt:.2f}") # Print values for check
    print(f"c_lbm={ce:.3f},\t tau1={tau1:.4f},\t tau2={tau2:.4f}")
    print(f"s1={1.0/tau1:.4f},\t s2={1.0/tau2:.4f}")
    print(f"nt={nt:.3f}\t, nue={nue:.4f}")
    #--------------------------------- Initialization - LBM - Numerical Arrays -----------------------------------------
    u=np.zeros((Nx),dtype="float64")
    u = A * np.sin(k_w*(x+dx/2) + phi)
    mx=np.zeros((Nx),dtype="float64")
    mxx=np.zeros((Nx),dtype="float64")
    f=np.zeros((3,Nx),dtype="float64")
    feq=np.zeros((3,Nx),dtype="float64")
    fp=np.zeros((3,Nx),dtype="float64")
    for k in range(0,3):
        f[k,:]=w[k]*u[:]*(1.0 + cx[k]*ce/cs**2 + ce*ce*(cx[k]**2-cs**2)/(2.0*cs**4))
        fp[k,:]=w[k]*u[:]*(1.0 + cx[k]*ce/cs**2 + ce*ce*(cx[k]**2-cs**2)/(2.0*cs**4))
    #----------------------------------------- Init Loop --------------------------------------------------------
    for t in range(nt*10):
        #========================================= LBM - Solution =============================================
        #--------------------Collision----------------
        for k in range(0,3):
            feq[k,:]= w[k]*u[:]*(1.0 + cx[k]*ce/cs**2 + ce*ce*(cx[k]**2-cs**2)/(2.0*cs**4))
        mx = np.einsum('i,ix->x', cx, f-feq)
        mxx = np.einsum('i,ix->x', 3.0*cx*cx-2.0, f-feq)
        for k in range(0,3):
            fp[k,:]= feq[k,:] + w[k]*( (1.0-1.0/tau1)*cx[k]*mx[:]/cs**2 + (1.0-1.0/tau2)*(cx[k]**2-cs**2)*mxx[:]/(2.0*cs**2) )
        #-----------------streaming-------------------
        for k in range(0,3):
            f[k,:]=np.roll(fp[k,:], cx[k], axis=0)
    #----------------------------------------- Maind Loop --------------------------------------------------------
    for t in range(nt):
        #========================================= LBM - Solution =============================================
        #--------------------Collision----------------
        for k in range(0,3):
            feq[k,:]= w[k]*(u[:]+cx[k]*(ce*u[:])/cs**2 + u[:]*ce*ce*(cx[k]**2-cs**2)/(2.0*cs**4))
        mx = np.einsum('i,ix->x', cx, f-feq)
        mxx = np.einsum('i,ix->x', 3.0*cx*cx-2.0, f-feq)
        for k in range(0,3):
            # fp[k,:]= f[k,:] - w[k]*( cx[k]*mx[:]/(tau1*cs**2) + (cx[k]**2-cs**2)*mxx[:]/(2.0*cs**2*tau2) )
            fp[k,:]= feq[k,:] + w[k]*( (1.0-1.0/tau1)*cx[k]*mx[:]/cs**2 + (1.0-1.0/tau2)*(cx[k]**2-cs**2)*mxx[:]/(2.0*cs**2) )
        #-----------------streaming-------------------
        for k in range(0,3):
            f[k,:]=np.roll(fp[k,:], cx[k], axis=0)
        #----------------------Macro------------------
        # u[:]=f[0,:]+f[1,:]+f[2,:]
        u=np.einsum('ix->x', f)
    #---------------------save-field-cases--------------
    u_cases[case][:]=u[:]

../../_images/583c87ad0e1fb91fc644591f181d814a8dc25291cf4934460e7edb21d65f03b0.png

../../_images/47c93b24577a7318891c2625a5ba62935abbe84cd0e1485958ab6894d209dc52.png

Non-Linear Convective Equation

Expanding in Taylor series the FD form LB equation, we can truncate the equation in second-order errors of $\Delta x$-scale:

\[ 0=- \Delta_{t}^{2} \frac{\partial^{2}{B}}{\partial{x}\partial{t}} s_{1} \left(s_{2} - 1\right) + \frac{\Delta_{t}^{2} \frac{\partial^{2}{u}}{\partial{t^{2}}} \left(- 3 s_{1} s_{2} + 2 s_{1} + 2 s_{2}\right)}{2} + \Delta_{t} \frac{\partial{B}}{\partial{x}} s_{1} s_{2} + \Delta_{t} \frac{\partial{u}}{\partial{t}} s_{1} s_{2} + \frac{\Delta_{x}^{2} s_{2} \left(\frac{3 \Delta_{t}^{2} \frac{\partial^{2}{\Lambda}}{\partial{x^{2}}} s_{1}}{\Delta_{x}^{2}} - \frac{6 \Delta_{t}^{2} \frac{\partial^{2}{\Lambda}}{\partial{x^{2}}}}{\Delta_{x}^{2}} \right)}{6} \]

isolating $\frac{\partial u}{\partial t}$:

\[ \frac{\partial u}{\partial t}=\Delta_{t} \frac{\partial^{2}{B}}{\partial{x}\partial{t}} \left(1 - \frac{1}{s_{2}}\right) + \Delta_{t} \frac{\partial^{2}{\Lambda}}{\partial{x^{2}}} \left(- \frac{1}{2} + \frac{1}{s_{1}}\right) + \Delta_{t} \frac{\partial^{2}{u}}{\partial{t^{2}}} \left(\frac{3}{2} - \frac{1}{s_{2}} - \frac{1}{s_{1}}\right) - \frac{\partial{B}}{\partial{x}} \]

To correct reduce the above equation and isolate the second-order of spatial derivative terms, we employ the relations:

\[ \displaystyle\frac{\partial u}{\partial t} = - \frac{\partial B}{\partial x} \]

\[ \displaystyle\frac{\partial^{2}u}{\partial t^{2}} = - \frac{\partial^{2}B}{\partial x\partial t} = -\frac{\partial }{\partial x}\left( \frac{\partial u}{\partial t}B_{u} \right) = - \frac{\partial }{\partial x}\left[ B_{u} \frac{\partial B}{\partial x} \right] \]

replacing the above relations in the second-order expanded equation:

\[ \frac{\partial u}{\partial t}=\Delta_{t} \frac{\partial^{2}{B}}{\partial{x}\partial{t}} \left(- \frac{1}{2} + \frac{1}{s_{1}}\right) + \Delta_{t} \frac{\partial^{2}{\Lambda}}{\partial{x^{2}}} \left(- \frac{1}{2} + \frac{1}{s_{1}}\right) - \frac{\partial{B}}{\partial{x}} \]

import numpy as np
import matplotlib.pyplot as plt
from pylab import *
matplotlib.rcParams['mathtext.fontset'] = 'cm'
#***************************************InputParameters************************************************
Nx=26 #Square Domain Length
#***************************************Lattice-Properties-D1Q3*************************************************
cs=1.0/sqrt(3.0);
w = np.array([4.0/6.0, 1.0/6.0, 1.0/6.0],dtype="float64")
cx = np.array([0, 1, -1],dtype="int16")  
#***************************************LBM-Scale************************************************
L = 20.0            # Length of the reservoir
T = 10.0             # Total simulation time
dx = L / (Nx-1)       # Spatial step size
c=2**(2)            # c=dx/dt
dt=dx/c             # Time step size
nt = int(T / dt)    # Time step number
uo=1.0              # Constant Fluid Flux
m=1.0               # Dynamic Viscosity ration
print('dx=',dx,'\t dt=',dt)
#***************************************Initialization************************************************
ue=uo/c
sk=np.zeros((Nx),dtype="float64")
Bx=np.zeros((Nx),dtype="float64")
f=np.zeros((3,Nx),dtype="float64")
for k in range(0,3):
    f[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]*3.0)
#***************************************Loop************************************************
for t in range(nt):
    #----------------------Macro------------------
    sk[:]=f[0,:]+f[1,:]+f[2,:]
    Bx=ue*sk*sk/(sk*sk+(1.0-sk)**(2.0))
    #--------------------Collision----------------
    for k in range(0,3):
        f[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]/cs**2)
    #-----------------streaming-------------------
    for k in range(0,3):
        f[k,:]=np.roll(f[k,:], cx[k], axis=0)
    #-----------------Boundaries-----------------------
    f[1,0]= 1.0 - f[0,0]-f[2,0]
    f[0,Nx-1]=f[0,Nx-2]
    f[1,Nx-1]=f[1,Nx-2]
    f[2,Nx-1]=f[2,Nx-2]
x = np.linspace(0, L, Nx)
plt.plot(x, sk, 'k-' ,label='Water Saturation $(t=10)$')
plt.xlabel(r"$x$",fontsize=18)
plt.ylabel(r"$u(x,t)$",fontsize=18)
plt.xlim(0,20)
plt.ylim(0,1)
plt.legend()
plt.grid(True)
plt.tight_layout()

dx= 0.8 	 dt= 0.2

../../_images/33c44822e5ed1ddeff022104943fda7da875fbd07b23306be2bb17439b4c6e82.png

%pylab inline
#***************************************Allocating Arrys************************************************
#--------------------------------------------Int--------------------------------------------------------
n=int(Nx);
cx=np.ones(3,dtype="int8")*0;
#------------------------------------------Double-------------------------------------------------------
w=np.ones(3,dtype="float64")*0
f=np.ones((3,n),dtype="float64")*0;
feq=np.ones((3,n),dtype="float64")*0;
fneq=np.ones((3,n),dtype="float64")*0;
sk=np.ones((n),dtype="float64")*0
Bx=np.ones((n),dtype="float64")*0
Pix=np.ones((n),dtype="float64")*0
La=np.ones((n),dtype="float64")*0
dsdx=np.ones((n),dtype="float64")*0
#********************************************Init-Parameters*********************************************
cs=1.0/sqrt(3.0);
w[0] = 4.0 / 6.0; w[1] = 1.0 / 6.0; w[2] = 1.0 / 6.0;
cx[0] = 0; cx[1] = 1; cx[2] = -1;
uo=ue
#---------------------Initialization-----------------------
for i in range(0,int(Nx/3)):
    sk[i]=1.0-3.0*i/Nx
    #sk[i]=0.0
Bx=uo*sk**(2.0)/(sk**(2.0)+(1.0-sk)**(2.0))
aa=sk*sk+(1.0-sk)**2
Gamma=uo*uo*( (sk-0.5)*(0.5-1.0/(6.0*aa)*(1.0-2.0/aa))/aa + 0.5*np.arctan(2.0*sk-1.0) )
Psi=-4.0*(sk**(3.0)/3.0-sk**(2.0)/2.0)
La=Gamma+cs*cs*Psi
for k in range(0,3):
    #f[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]*3.0)
    f[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]*3.0+4.5*(La[:]-sk[:]/3.0)*(cx[k]*cx[k]-1.0/3.0) )
    #f[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]*3.0+4.5*(La[:])*(cx[k]*cx[k]-1.0/3.0) )
#--------------------------------MainLoop------------------------------
import time
mstep=160
start=time.time()
for kk in range(0,mstep):
#for kk in range(0,1):
    #----------------------Macro------------------
    sk[:]=f[0,:]+f[1,:]+f[2,:]
    Bx=uo*sk*sk/(sk*sk+(1.0-sk)**(2.0))
    aa=sk*sk+(1.0-sk)**2
    Gamma=uo*uo*((sk-0.5)*(0.5-1.0/(6.0*aa)*(1.0-2.0/aa))/aa + 0.5*np.arctan(2.0*sk-1.0))
    Psi=-4.0*(sk**(3.0)/3.0-sk**(2.0)/2.0)
    La=Gamma+cs*cs*Psi
    #--------------------Collision----------------
    for k in range(0,3):
        feq[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]*3.0+4.5*(La[:]-sk[:]/3.0)*(cx[k]*cx[k]-1.0/3.0) )
        #feq[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]*3.0+4.5*(La[:])*(cx[k]*cx[k]-1.0/3.0) )
        #feq[k,:]=w[k]*(sk[:]+cx[k]*Bx[:]*3.0 )
    Pix= (f[1,:]-feq[1,:])*cx[1]+(f[2,:]-feq[2,:])*cx[2]
    for k in range(0,3):
        fneq[k,:]=w[k]*(cx[k]*Pix[:]*3.0)
    for k in range(0,3):
        #f[k,:]=feq[k,:] +(1.0-1.0/tau)*fneq[k,:]
        f[k,:]=f[k,:] -(1.0/tau)*(f[k,:]-feq[k,:])
        #f[k,:]=feq[k,:] 
    #-----------------streaming-------------------
    for k in range(0,3):
        f[k,:]=np.roll(f[k,:], cx[k], axis=0)
    #-----------------Boundaries-----------------------
    f[1,0]= 1.0 - f[0,0]-f[2,0]
    f[0,n-1]=f[0,n-2]
    f[1,n-1]=f[1,n-2]
    f[2,n-1]=f[2,n-2]

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib

/usr/local/lib/python3.10/dist-packages/IPython/core/magics/pylab.py:166: UserWarning: pylab import has clobbered these variables: ['f']
`%matplotlib` prevents importing * from pylab and numpy
  warn("pylab import has clobbered these variables: %s"  % clobbered +

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 57
     54     fneq[k,:]=w[k]*(cx[k]*Pix[:]*3.0)
     55 for k in range(0,3):
     56     #f[k,:]=feq[k,:] +(1.0-1.0/tau)*fneq[k,:]
---> 57     f[k,:]=f[k,:] -(1.0/tau)*(f[k,:]-feq[k,:])
     58     #f[k,:]=feq[k,:] 
     59 #-----------------streaming-------------------
     60 for k in range(0,3):

NameError: name 'tau' is not defined

Non-Linear Convective-Diffusive Equation

Expanding in Taylor series the FD form LB equation, we can truncate the equation with errors isolated in third and fourth-order of scale of $\Delta x$:

\[\begin{split} \begin{array}{c} 0 = \frac{\Delta_{t}^{4} \frac{\partial^{4}{u}}{\partial{t^{4}}} \left(- 15 s_{1} s_{2} + 14 s_{1} + 14 s_{2} - 12\right)}{24} - \frac{\Delta_{t}^{3} \Delta_{x} \frac{\partial^{4}{B}}{\partial{x}\partial{t^{3}}} s_{1} \left(s_{2} - 1\right)}{6} + \frac{\Delta_{t}^{3} \frac{\partial^{3}{u}}{\partial{t^{3}}} \left(7 s_{1} s_{2} - 6 s_{1} - 6 s_{2} + 6\right)}{6} \\ + \frac{\Delta_{t}^{2} \Delta_{x}^{2} \left(\frac{\partial^{4}{C}}{\partial{x^{2}}\partial{t^{2}}} s_{1} s_{2} - \frac{\partial^{4}{C}}{\partial{x^{2}}\partial{t^{2}}} s_{2} + 3 \frac{\partial^{4}{\Lambda}}{\partial{x^{2}}\partial{t^{2}}} s_{1} s_{2} - 3 \frac{\partial^{4}{\Lambda}}{\partial{x^{2}}\partial{t^{2}}} s_{2} - 3 \frac{\partial^{4}{u}}{\partial{x^{2}}\partial{t^{2}}} s_{1} - 3 \frac{\partial^{4}{u}}{\partial{x^{2}}\partial{t^{2}}} s_{2} + 6 \frac{\partial^{4}{u}}{\partial{x^{2}}\partial{t^{2}}}\right)}{12} \\+ \frac{\Delta_{t}^{2} \Delta_{x} \frac{\partial^{3}{B}}{\partial{x}\partial{t^{2}}} s_{1} \left(s_{2} - 1\right)}{2} + \frac{\Delta_{t}^{2} \frac{\partial^{2}{u}}{\partial{t^{2}}} \left(- 3 s_{1} s_{2} + 2 s_{1} + 2 s_{2}\right)}{2} - \frac{\Delta_{t} \Delta_{x}^{3} \frac{\partial^{4}{B}}{\partial{x^{3}}\partial{t}} s_{1} \left(s_{2} - 1\right)}{6} + \frac{\Delta_{t} \Delta_{x}^{2} \left(- \frac{\partial^{3}{C}}{\partial{x^{2}}\partial{t}} s_{1} s_{2} + \frac{\partial^{3}{C}}{\partial{x^{2}}\partial{t}} s_{2} - 3 \frac{\partial^{3}{\Lambda}}{\partial{x^{2}}\partial{t}} s_{1} s_{2} + 3 \frac{\partial^{3}{\Lambda}}{\partial{x^{2}}\partial{t}} s_{2} + 3 \frac{\partial^{3}{u}}{\partial{x^{2}}\partial{t}} s_{1} + 3 \frac{\partial^{3}{u}}{\partial{x^{2}}\partial{t}} s_{2} - 6 \frac{\partial^{3}{u}}{\partial{x^{2}}\partial{t}}\right)}{6} - \Delta_{t} \Delta_{x} \frac{\partial^{2}{B}}{\partial{x}\partial{t}} s_{1} \left(s_{2} - 1\right) \\ + \Delta_{t} \frac{\partial{u}}{\partial{t}} s_{1} s_{2} + \frac{\Delta_{x}^{4} s_{2} \left(\frac{\partial^{4}{C}}{\partial{x^{4}}} s_{1} - 2 \frac{\partial^{4}{C}}{\partial{x^{4}}} + 3 \frac{\partial^{4}{\Lambda}}{\partial{x^{4}}} s_{1} - 6 \frac{\partial^{4}{\Lambda}}{\partial{x^{4}}}\right)}{72} + \frac{\Delta_{x}^{3} \frac{\partial^{3}{B}}{\partial{x^{3}}} s_{1} s_{2}}{6} + \frac{\Delta_{x}^{2} s_{2} \left(\frac{\partial^{2}{C}}{\partial{x^{2}}} s_{1} - 2 \frac{\partial^{2}{C}}{\partial{x^{2}}} + 3 \frac{\partial^{2}{\Lambda}}{\partial{x^{2}}} s_{1} - 6 \frac{\partial^{2}{\Lambda}}{\partial{x^{2}}}\right)}{6} + \Delta_{x} \frac{\partial{B}}{\partial{x}} s_{1} s_{2} \end{array} \end{split}\]

isolating $\frac{\partial u}{\partial t}$:

(67)\[\begin{split} \begin{array}{c} \displaystyle\frac{\partial u}{\partial t} + \frac{\Delta_{x}}{\Delta_{t}}\frac{\partial{B}}{\partial{x}} = \nu\frac{\partial^{2}{C}}{\partial{x^{2}}} + \Delta_{t}^{3} \frac{\partial^{4}{u}}{\partial{t^{4}}} \left(\frac{5}{8} - \frac{7}{12 s_{2}} - \frac{7}{12 s_{1}} + \frac{1}{2 s_{1} s_{2}}\right) + \Delta_{t}^{2} \Delta_{x} \frac{\partial^{4}{B}}{\partial{x}\partial{t^{3}}} \left(\frac{1}{6} - \frac{1}{6 s_{2}}\right) + \Delta_{t}^{2} \frac{\partial^{3}{u}}{\partial{t^{3}}} \left(- \frac{7}{6} + \frac{1}{s_{2}} + \frac{1}{s_{1}} - \frac{1}{s_{1} s_{2}}\right) + \Delta_{t} \Delta_{x}^{2} \frac{\partial^{4}{C}}{\partial{x^{2}}\partial{t^{2}}} \left(- \frac{1}{12} + \frac{1}{12 s_{1}}\right)\\ + \Delta_{t} \Delta_{x}^{2} \frac{\partial^{4}{\Lambda}}{\partial{x^{2}}\partial{t^{2}}} \left(- \frac{1}{4} + \frac{1}{4 s_{1}}\right) + \Delta_{t} \Delta_{x}^{2} \frac{\partial^{4}{u}}{\partial{x^{2}}\partial{t^{2}}} \left(\frac{1}{4 s_{2}} + \frac{1}{4 s_{1}} - \frac{1}{2 s_{1} s_{2}}\right) + \Delta_{t} \Delta_{x} \frac{\partial^{3}{B}}{\partial{x}\partial{t^{2}}} \left(- \frac{1}{2} + \frac{1}{2 s_{2}}\right) + \Delta_{t} \frac{\partial^{2}{u}}{\partial{t^{2}}} \left(\frac{3}{2} - \frac{1}{s_{2}} - \frac{1}{s_{1}}\right) + \Delta_{x}^{3} \frac{\partial^{4}{B}}{\partial{x^{3}}\partial{t}} \left(\frac{1}{6} - \frac{1}{6 s_{2}}\right) \\+ \Delta_{x}^{2} \frac{\partial^{3}{C}}{\partial{x^{2}}\partial{t}} \left(\frac{1}{6} - \frac{1}{6 s_{1}}\right) + \Delta_{x}^{2} \frac{\partial^{3}{\Lambda}}{\partial{x^{2}}\partial{t}} \left(\frac{1}{2} - \frac{1}{2 s_{1}}\right) + \Delta_{x}^{2} \frac{\partial^{3}{u}}{\partial{x^{2}}\partial{t}} \left(- \frac{1}{2 s_{2}} - \frac{1}{2 s_{1}} + \frac{1}{s_{1} s_{2}}\right) + \Delta_{x} \frac{\partial^{2}{B}}{\partial{x}\partial{t}} \left(1 - \frac{1}{s_{2}}\right) + \frac{\Delta_{x}^{4} \frac{\partial^{4}{C}}{\partial{x^{4}}} \left(- \frac{1}{72} + \frac{1}{36 s_{1}}\right)}{\Delta_{t}} + \frac{\Delta_{x}^{4} \frac{\partial^{4}{\Lambda}}{\partial{x^{4}}} \left(- \frac{1}{24} + \frac{1}{12 s_{1}}\right)}{\Delta_{t}} - \frac{\Delta_{x}^{3} \frac{\partial^{3}{B}}{\partial{x^{3}}}}{6 \Delta_{t}} + \frac{\Delta_{x}^{2} \frac{\partial^{2}{\Lambda}}{\partial{x^{2}}} \left(- \frac{1}{2} + \frac{1}{s_{1}}\right)}{\Delta_{t}} \end{array} \end{split}\]

where $\nu=c_{s}^{2}\left(\displaystyle\frac{2-s1}{2s1}\right)\displaystyle\frac{\Delta x^{2}}{\Delta t}$

To correct reduce the above equation and isolate the third and fourth order of spatial derivative terms, we employ the relations:

\[ \displaystyle\frac{\partial u}{\partial t} = \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial B}{\partial x} \]

\[ \displaystyle\frac{\partial^{3}u}{\partial x^{2}\partial t} = \nu \displaystyle\frac{\partial^{4}C}{\partial x^{4}} - \frac{\partial^{3}B}{\partial x^{3}} \]

\[ \displaystyle\frac{\partial^{2}u}{\partial t^{2}} = \nu \displaystyle\frac{\partial^{3}C}{\partial x^{2}\partial t} - \frac{\partial^{2}B}{\partial x\partial t} = \nu \displaystyle \frac{\partial^{2}}{\partial x^{2}}\left(\frac{\partial u}{\partial t}C_{u} \right) - \frac{\partial }{\partial x}\left( \frac{\partial u}{\partial t}B_{u} \right) = \nu \displaystyle \frac{\partial^{2}}{\partial x^{2}}\left[C_{u}\left( \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial B}{\partial x} \right) \right] - \frac{\partial }{\partial x}\left[ B_{u}\left( \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial B}{\partial x} \right) \right] \]

\[ \displaystyle\frac{\partial^{2}u}{\partial t^{2}} = \nu \displaystyle\frac{\partial^{3}C}{\partial x^{2}\partial t} - \frac{\partial^{2}B}{\partial x\partial t} = \nu \displaystyle \frac{\partial^{2}}{\partial x^{2}}\left(\frac{\partial u}{\partial t}C_{u} \right) - \frac{\partial }{\partial x}\left( \frac{\partial u}{\partial t}B_{u} \right) = \nu \displaystyle \frac{\partial^{2}}{\partial x^{2}}\left[C_{u}\left( \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial B}{\partial x} \right) \right] - \frac{\partial }{\partial x}\left[ \nu B_{u} \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial \Lambda}{\partial x} \right] \]

\[ \displaystyle\frac{\partial^{4}u}{\partial x^{2}\partial t^{2}} = \nu \displaystyle\frac{\partial^{5}C}{\partial x^{4}\partial t} - \frac{\partial^{4}B}{\partial x^{3}\partial t} = \nu \displaystyle \frac{\partial^{4}}{\partial x^{4}}\left(\left( \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial B}{\partial x} \right)C_{u} \right) - \frac{\partial^{3} }{\partial x^{3}}\left( \left( \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial B}{\partial x} \right)B_{u} \right) \]

\[ \displaystyle\frac{\partial^{3}u}{\partial t^{3}} = \nu \displaystyle\frac{\partial^{4}C}{\partial x^{2}\partial t^{2}} - \frac{\partial^{3}B}{\partial x\partial t^{2}} = \nu \frac{\partial^{2}}{\partial x^{2}}\left[ C_{uu} \left(\frac{\partial u}{\partial t}\right)^{2} + C_u \frac{\partial^{2}u}{\partial t^{2}} \right] - \frac{\partial}{\partial x}\left[B_{uu}\left(\frac{\partial u}{\partial t}\right)^{2} + B_u \frac{\partial^{2}u}{\partial t^{2}} \right] \]

\[ \displaystyle\frac{\partial^{4}u}{\partial t^{4}} = \nu \displaystyle\frac{\partial^{5}C}{\partial x^{2}\partial t^{3}} - \frac{\partial^{4}B}{\partial x\partial t^{3}} = \nu \frac{\partial^{2}}{\partial x^{2}} \left[ C_{uuu}\left(\frac{\partial u}{\partial t}\right)^{3} + 3C_{uu}\left(\frac{\partial u}{\partial t}\right)\left(\frac{\partial^{2}u}{\partial t^{2}}\right) + C_u \frac{\partial^{3}u}{\partial t^{3}} \right] - \frac{\partial}{\partial x} \left[B_{uuu}\left(\frac{\partial u}{\partial t}\right)^{3} + 3B_{uu}\left(\frac{\partial u}{\partial t}\right)\left(\frac{\partial^{2}u}{\partial t^{2}}\right) + B_u \frac{\partial^{3}u}{\partial t^{3}} \right] \]

\[ \displaystyle\frac{\partial B}{\partial t} = B_{u}\frac{\partial u}{\partial t} = B_{u}\left( \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial B}{\partial x} \right) = \nu B_{u} \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - B_{u}^{2}\frac{\partial u}{\partial x} = \nu B_{u} \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial \Lambda}{\partial x} \]

\[ \displaystyle\frac{\partial^{2} B}{\partial x \partial t} = \frac{\partial }{\partial x}\left(B_{u}\frac{\partial u}{\partial t}\right) = \frac{\partial B_{u}}{\partial x} \frac{\partial u}{\partial t} + B_{u}\frac{\partial u}{\partial x\partial t} = \frac{\partial B_{u}}{\partial x} \left( \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial B}{\partial x} \right) + B_{u}\left( \nu \displaystyle\frac{\partial^{3}C}{\partial x^{3}} - \frac{\partial^{2} B}{\partial x^{2}} \right) = \frac{\partial B_{u}}{\partial x} \left( \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}}\right) + B_{u}\left( \nu \displaystyle\frac{\partial^{3}C}{\partial x^{3}} \right) - \frac{\partial^{2} \Lambda}{\partial x^{2}} \]

\[ \displaystyle\frac{\partial C}{\partial t} = C_{u}\frac{\partial u}{\partial t} = C_{u}\left( \nu \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - \frac{\partial B}{\partial x} \right) = \nu C_{u} \displaystyle\frac{\partial^{2}C}{\partial x^{2}} - C_{u}\frac{\partial B}{\partial x} \]