# A Truly Second-Order and Unconditionally Stable Thermal Lattice Boltzmann Method
## 2 Simplified Thermal Lattice Boltzmann Method (STLBM) without Evolution of Distribution Functions  
### Thermal Lattice Boltzmann Method (TLBM) and Chapman-Enskog (C-E) Expansion Analysis  
In incompressible thermal fluids, effects of temperature and velocity are coupled. Within LBM framework, one popular approach to interpret such coupling effect is to introduce internal energy distribution function [13]. This approach is named as the thermal lattice Boltzmann method (TLBM) and the evolutions of the distribution functions can be expressed as  
\[f_{\alpha}(\mathbf{r}+\mathbf{e}_{\alpha}\delta_{t},\,t+\delta_{t})=f_{\alpha} (\mathbf{r},\,t)+\frac{f_{\alpha}^{eq}(\mathbf{r},\,t)-f_{\alpha}(\mathbf{r}, \,t)}{\tau_{v}},\,\,\,\alpha=0,\,1,\,...\,,\,M \tag{1}\]  
\[g_{\alpha}(\mathbf{r}+\mathbf{e}_{\alpha}\delta_{t},\,t+\delta_{t})=g_{\alpha} (\mathbf{r},\,t)+\frac{g_{\alpha}^{eq}(\mathbf{r},\,t)-g_{\alpha}(\mathbf{r}, \,t)}{\tau_{c}},\,\,\,\alpha=0,\,1,\,...\,,\,N \tag{2}\]  
where \(f_{\alpha}\) and \(g_{\alpha}\) represent the density distribution function and the internal energy distribution function along the \(\alpha\) direction, respectively, with the superscript \"\(eq\)\" denoting their equilibrium state; \(\tau_{v}\) and \(\tau_{c}\) are single relaxation parameters related to the kinematic viscosity and the thermal diffusivity, respectively; \(\mathbf{e}_{\alpha}\) and \(\delta_{t}\) represent the lattice velocity and the time step applied in the model; and \(M\) and \(N\) are the number of the lattice velocities used in the density and the internal energy distribution functions, respectively. The equilibrium density distribution function is given by  
\[f_{\alpha}^{eq}=\rho\omega_{\alpha}\bigg{[}1+\frac{\mathbf{e}_{\alpha}\cdot \mathbf{u}}{c_{s}^{2}}+\frac{\left(\mathbf{e}_{\alpha}\cdot\mathbf{u}\right) ^{2}-\left(c_{s}|\mathbf{u}|\right)^{2}}{2c_{s}^{4}}\bigg{]} \tag{3}\]while the equilibrium internal energy distribution functions in two and three dimensions are, respectively, expressed as  
\[g_{a}^{eq}(\textbf{r},t)=\left\{ \begin{array}{ll}-\frac{2\rho e}{3}\frac{| \textbf{u}|^{2}}{c^{2}}&\alpha=0\\ \frac{\rho e}{9}\left[\frac{3}{2}+\frac{3}{2}\frac{\textbf{e}_{\alpha}\cdot \textbf{u}}{c^{2}}+\frac{9}{2}\frac{\left(\textbf{e}_{\alpha}\cdot\textbf{u} \right)^{2}}{c^{4}}-\frac{3}{2}\frac{|\textbf{u}|^{2}}{c^{2}}\right]&\alpha=1- 4\\ \frac{\rho e}{36}\left[3+6\cdot\frac{\textbf{e}_{\alpha}\cdot\textbf{u}}{c^{2} }+\frac{9}{2}\frac{\left(\textbf{e}_{\alpha}\cdot\textbf{u}\right)^{2}}{c^{4} }-\frac{3}{2}\frac{|\textbf{u}|^{2}}{c^{2}}\right]&\alpha=5-8\end{array} \right. \tag{4}\]  
\[g_{a}^{eq}(\textbf{r},t)=\left\{ \begin{array}{ll}-\frac{\rho e}{3}\frac{| \textbf{u}|^{2}}{c^{2}}&\alpha=0\\ \frac{\rho e}{9}\left[1+\frac{\textbf{e}_{\alpha}\cdot\textbf{u}}{c^{2}}+ \frac{9}{2}\frac{\left(\textbf{e}_{\alpha}\cdot\textbf{u}\right)^{2}}{c^{4}}- \frac{3}{2}\frac{|\textbf{u}|^{2}}{c^{2}}\right]&\alpha=1-6\\ \frac{\rho e}{72}\left[3+7\cdot\frac{\textbf{e}_{\alpha}\cdot\textbf{u}}{c^{2} }+\frac{9}{2}\frac{\left(\textbf{e}_{\alpha}\cdot\textbf{u}\right)^{2}}{c^{4} }-\frac{3}{2}\frac{|\textbf{u}|^{2}}{c^{2}}\right]&\alpha=7-14\end{array} \right. \tag{5}\]  
where \(c=\delta_{x}/\delta_{t}\) is set to be 1 in this work; and \(\delta_{x}\) is the lattice spacing. The weighting coefficient \(\omega_{a}\) and the sound speed \(c_{s}\) are  
\[\left\{ \begin{array}{ll}D2Q9:&\omega_{0}=\frac{4}{9}&\omega_{1-4}=\frac{1}{ 9}&\omega_{5-8}=\frac{1}{36}\\ D3Q15:&\omega_{0}=\frac{2}{9}&\omega_{1-6}=\frac{1}{9}&\omega_{7-14}=\frac{1}{ 72}\end{array}