# A Truly Second-Order and Unconditionally Stable Thermal Lattice Boltzmann Method
## 3 Unconditionally Stable Thermal Lattice Boltzmann Method (USTLBM)
### Stability Analysis  
Theoretical analysis for the stability of lattice Boltzmann method (LBM) was firstly conducted by [PERSON] and [PERSON] [18]. It was found that standard LBM suffers from numerical instability when the single relaxation parameter \(\tau\) is close to 0.5, which indicates that standard LBM is not suitable to problems at high Reynolds numbers. The recently developed simplified thermal lattice Boltzmann method (STLBM), on the other hand, performs very well in numerical stability. Through von Neumann stability analysis, it can be proven that STLBM is unconditionally stable [12]. The USTLBM proposed in this paper can be regarded as an improved version of STLBM that ensures higher order of accuracy in time. In this section, we use von Neumann stability analysis to study its stability performance. For the ease of analysis, we take the two-dimensional formulations as the example here. Analysis on the three-dimensional formulations can be carried out through similar processes, and same conclusions can be obtained. To begin with, the following macroscopic variables are defined:  
\[y_{1}=\rho,\;\;y_{2}=\rho u,\;\;y_{3}=\rho v,\;\;y_{4}=\rho e \tag{48}\]  
The vector of these variables is then written as  
\[\mathbf{Y}=\left( \begin{array}{cccc}y_{1}&y_{2}&y_{3}&y_{4}\end{array} \right) ^{\mathsf{T}} \tag{49}\]  
The equilibrium distribution functions Equations (3) and (4) are functions of the macroscopic variables, which can be written as  
\[f_{\alpha}^{eq}=w_{\alpha}y_{1\alpha}+\frac{w_{\alpha}}{c_{s}^{2}}\big{(}e_{ \alpha x}y_{2\alpha}+e_{\alpha y}y_{3\alpha}\big{)}+\frac{w_{\alpha}}{2c_{s}^{ 4}}\frac{\big{(}e_{\alpha x}y_{2\alpha}+e_{\alpha y}y_{3\alpha}\big{)}^{2}}{y_{ 1\alpha}}-\frac{w_{\alpha}}{2c_{s}^{2}}\frac{y_{2\alpha}^{2}+y_{3\alpha}^{2}}{ y_{1\alpha}} \tag{50}\]  
\[g_{\alpha}^{eq}=\left\{ \begin{array}{ll}-\frac{2}{3}y_{4\alpha}\frac{y_{2 \alpha}^{2}+y_{3\alpha}^{2}}{y_{1\alpha}^{2}}&\alpha=0\\ \frac{1}{6}y_{4\alpha}+\frac{1}{6}y_{4\alpha}\frac{e_{\alpha x}y_{2\alpha}+e_{ \alpha y}y_{3\alpha}}{y_{1\alpha}}+\frac{1}{2}y_{4\alpha}\Big{(}\frac{e_{ \alpha x}y_{2\alpha}+e_{\alpha y}y_{3\alpha}}{y_{1\alpha}}\Big{)}^{2}-\frac{1} {6}y_{4\alpha}\frac{y_{2\alpha}^{2}+y_{3\alpha}^{2}}{y_{1\alpha}^{2}}&\alpha= 1\sim 4\\ \frac{1}{12}y_{4\alpha}+\frac{1}{6}y_{4\alpha}\frac{e_{\alpha y_{2\alpha}+e_{ \alpha y}y_{3\alpha}}}{y_{1\alpha}}+\frac{1}{8}y_{4\alpha}\Big{(}\frac{e_{ \alpha x}y_{2\alpha}+e_{\alpha y}y_{3\alpha}}{y_{1\alpha}}\Big{)}^{2}-\frac{1} {24}y_{4\alpha}\frac{y_{2\alpha}^{2}+y_{3\alpha}^{2}}{y_{1\alpha}^{2}}&\alpha= 5\sim 8\end{array} \right. \tag{51}\]  
where \(y_{1\alpha}\), \(y_{2\alpha}\), \(y_{3\alpha}\) and \(y_{4\alpha}\) denote the functional values of \(y_{1}\), \(y_{2}\), \(y_{3}\) and \(y_{4}\) at the space/time level of (**r**-**e**\({}_{\alpha}\), \(t\)-**\(\delta_{t}\)), respectively. Substituting Equations (50) and (51) into the expressions for the intermediate density, momentum and energy (Equations (39)-(41)) gives a set of nonlinear algebraic equations \(\mathbf{Y}^{*}=\mathbf{F(Y^{n})}\). To linearize the problem to perform von Neumann stability analysis, the obtained equations are differentiated with respect to \(y_{1}\), \(y_{2}\), \(y_{3}\) and \(y_{4}\), which leads to  
\[\mathbf{d}\mathbf{Y}^{*}=\left[\frac{\partial\mathbf{Y}^{*}}{\partial y_{1}} \right]^{n}\mathbf{d}y_{1}^{n}+\left[\frac{\partial\mathbf{Y}^{*}}{\partial y _{2}}\right]^{n}\mathbf{d}y_{2}^{n}+\left[\frac{\partial\mathbf{Y}^{*}}{ \partial y_{3}}\right]^{n}\mathbf{d}y_{3}^{n}+\left[\frac{\partial\mathbf{Y}^ {*}}{\partial y_{4}}\right]^{n}\mathbf{d}y_{4}^{n} \tag{52}\]where \(\mathbf{Y^{*}}=({y_{1}}^{*},{y_{2}}^{*},{y_{3}}^{*},{y_{4}}^{*})^{T}\); and the superscript \(n\) represents the time level.