Numerical Analysis for Time-Averaged Equilibrium of Lattice Boltzmann Method in Heat Diffusion and Sod Shock Tube Problem

1. Introduction

The computational method for solving engineering problems has become a powerful tool to explore and investigate chemical and physical phenomena. In 1956, the finite element method (FEM) was applied to solve a structural problem by Turner et al., [1]. At the same time the FDM and FEM were used to solve partial differential problems, heat transfer and fluid dynamics problems [2]. The finite volume method was introduced by Imperial College in 1980 to solve transport phenomena problems [3]. All of the aforementioned methods belong to the same family of weighted residual methods, and the only difference between these methods is the base and weight function [4]. The LBM emerged in 1988 to overcome the limitation of the lattice gas cellular automata [5], since then the LBM evolved as an alternative powerful method for solving fluid dynamic problems due to its kinetic-theory foundation, parallelizability, and algorithm simplicity [6,7], the main difficulty in LBM is the collision operator. The simplest approximation of the collision operator without producing significant error outcomes was introduced by Bhatnagar, Gross and Krook.

Over the past two decades, many modifications to the standard LBM have been proposed to improve accuracy, stability and applicability to capture complex flows such as the multiple relaxation time (MRT) has been introduced to reduce the numerical instability in compressible and turbulent flows, as demonstrated in previous studies [8,9,10]. the two-relaxation time model (TRT) was introduced as a simplified variant to balance between performance and efficiency by using two relaxation parameters symmetric and antisymmetric components [11]. In 1999, the Entropic LBM (ELBM) was introduced by Karlin et al., [12], to improve the stability in turbulent and high-speed flows. Cascade LBM, which is known as a central moment method, reformulates the collision using central rather than raw moments to improve Galilean invariance and accuracy. This advancement was achieved by altering the physical model of LBM and the collision structure [13].

There remains a lack of focused research on modifying the equilibrium distribution function in the BGK model by using a time-averaging technique, and on investigating the effect of this averaging on the behavior of the heat diffusion and more broadly, diffusion-type equations, the idea of transforming the FVM into a BGK-like LBM structure, combined with time averaging of the equilibrium distribution function to enable efficient parallelization using the Message Passing Interface with Compute Unified Device Architecture MPI-CUDA, has not yet been explored.

The significance of this study lies in introducing a hybrid modelling approach that bridges numerical formulation and computational implementation. By combining time averaging with BGK and leveraging structural transformation, the work offers a new direction for scalable and stable simulation in both diffusive and compressible regimes.

Therefore, the objective of this study is to propose a novel time-averaged formulation of the equilibrium distribution function within the BGK collision model of the LBM, the aim is to reduce the number of time steps required for the outcomes while maintaining the accuracy, this can be validated through two domains. First, heat diffusion to assess the performance in parabolic transport regime, second the sod shock tube is employed to assess the method’s behavior in compressible flow regime by transforming the mathematical formulation of FVM into an equivalent LBM framework, focusing on preserving the mathematical structure and numerical behavior rather than replicating the physical modeling.

2. Methodology

This study starts from the conventional heat diffusion equation without source term. The equation is first discretized using an explicit FDM approach, employing a forward difference in time and a central difference in space. The discrete form is then compared with LBM formulation with BGK approximation, to establish a connection between the two methods, a mathematical technique is applied to resemble the LBM with the BGK approximation from the FDM, once this resemblance is achieved, we introduce the new formulation of the equilibrium distribution function based on time averaging.

This approach needs information from a future time step, which we don’t have it directly, to overcome this, the standard BGK model is used to estimate the future state, which is then used to recalculate the time-averaging equilibrium distribution function.

A flow diagram is provided to illustrate the proposed computational process, this approach is first tested on the heat diffusion problem and then extended to Sod shock tube, where specific mathematical transformations are applied to the FVM to transform it into BGK-like relaxation form.

Finally, numerical results are generated and compared against conventional methods to assess the performance of the proposed model. Figures are used to visualize the difference and evaluate the effectiveness of the time-averaged equilibrium approach in both diffusive and compressible regimes.

2.1. Governing Equation

We start with the conventional heat diffusion equation in one dimension, without a source term:

$${\displaystyle \frac{\partial T}{\partial t}}=k{\displaystyle \frac{{\partial }^{2}T}{\partial x^2}}$$

(1)

where $T$ is the temperature, $t$ is time, $x$ is the spatial coordinate, and $k$ is the thermal diffusivity. This equation is discretized using an explicit FDM, employing a forward difference in time and a central difference in space:

$$T_i^{t+1}=T_i^t+{\displaystyle \frac{k\Delta t}{{\Delta x}^2}}(T_{i+1}^t-2T_i^t+T_{i-1}^t)$$

(2)

where $i$ is the spatial node index.

The LBM for heat diffusion equation using BGK approximation can be written as follows:

$$f_i\left(x+,c_i\Delta t, t+\Delta t\right)=(1-w)f_i\left(x, t\right)+w{f}_i^{eq}\left(x, t\right)$$

(3)

where $f_i$ is the distribution function in direction $i$ , $c_i$ is the lattice velocity, $w$ is the relaxation time, $f_i^{eq}$ is the equilibrium distribution function, $\Delta t$ is the time step. We rearrange the terms in Eq. (2) to match the structure in Eq. (3) and extract the analogy between the two formulations (Mohamad 2011).

2.2. Resembling Procedure

The Eq. (2) can be written as follows:

$$T_i^{t+1}=T_i^t(1-{\displaystyle \frac{2k\Delta t}{{\Delta x}^2}})+{\displaystyle \frac{2k\Delta t}{{\Delta x}^2}}\left({\displaystyle \frac{T_{i+1}^t+T_{i-1}^t}{2}}\right)$$

(4)

Assume:

$$B={\displaystyle \frac{2k\Delta t}{{\Delta x}^2}}$$

(5)

Substituting Eq. (5) into Eq. (4) gives:

$$T_i^{t+1}=(1-B)T_i^t+B\left({\displaystyle \frac{T_{i+1}^t+T_{i-1}^t}{2}}\right)$$

(6)

This mathematical manipulation to compare the FDM with LBM approximation of heat diffusion, the last term in Eq. (5) $\left({\displaystyle \frac{T_{i+1}^t+T_{i-1}^t}{2}}\right)$ represents the spatially averaged temperature around $T_i^t$, because $T_i^t$ lies at the mid-distance between $T_{i-1}^t$ and $T_{i+1}^t$, in other word this is the equilibrium value of the $T_i$, Eq. (6) can be rewritten as follows:

$$T_i^{t+1}=(1-B)T_i^t+B\left(T_{i,t}^{eq}\right)$$

(7)

where $T_{i,t}^{eq}$ is the equilibrium value of the temperature at specific location and time $i,t$. By comparing Eq. (7) and (3) we can assume the $f_i^{eq}$ is an average value around $f_i\left(x, t\right)$ and lies at the mid-distance between ${\displaystyle \frac{f_i\left(x+\Delta x, t\right)+f_i\left(x-\Delta x, t\right)}{2}}$ , this analogy is very important for the current study, as it supports the next assumption, which will be the pillar of this work.

2.3. Time-Averaged Equilibrium Distribution function

We propose a new form of the equilibrium distribution function, averaged overtime:

$$f_{x,t}^{eq}={\displaystyle \frac{f_i\left(x, t+1\right)+f_i\left(x, t-1\right)}{2}}$$

(8)

This approach assumes knowledge of future value of $f_{i,t}^{eq}$, which are not directly available. To resolve this, we first compute the standard BGK solution at $t+\Delta t$, then use the result to reconstruct the time-averaged equilibrium.

2.4. Computational Flow Process

A flow diagram is constructed to summarize the step-by-step process of implementing the time-averaged equilibrium in the BGK model, the process starts with computing the BGK result at $t+\Delta t$, followed by calculating the time-averaged equilibrium, updating the solution and iteration.

Figure 1. Computational flowchart illustrating the time-averaged BGK procedure used in this study.

Figure 1. Computational flowchart illustrating the time-averaged BGK procedure used in this study.

2.5. Time-Averaged BGK Formulation for Heat Diffusion

We consider the heat diffusion equation modeled using the LBM with the BGK approximation. The chosen lattice structure is D1Q2, the simulation is conducted in lattice units, where both the spatial and temporal resolutions are normalized to unity. Under this assumption, the discrete velocity set for D1Q2 model is defined as:

$$\Delta x= \Delta t=1, c_1=1, c_2=-1, i=\mathrm{1,2}$$

(9)

Substituting Eq. (9) into Eq. (3) leads to:

$$f_1\left(x+1, t+1\right)=(1-w)f_1\left(x, t\right)+w{f}_1^{eq}\left(x, t\right)$$

(10)

$$f_2\left(x-1, t+1\right)=(1-w)f_2\left(x, t\right)+w{f}_2^{eq}\left(x, t\right)$$

(11)

where:

$$w={\displaystyle \frac{1}{k+0.5}}$$

(12)

The macroscopic temperature is calculated as the sum of two particles populations:

$$T\left(x,t\right)=f_1\left(x,t\right)+ f_2\left(\mathrm{x},\mathrm{t}\right)$$

(13)

For D1Q2 the equilibrium functions are defined as:

$$f_1^{eq}\left(x, t\right)=0.5\ast T\left(x,t\right), f_2^{eq}\left(x, t\right)=0.5\ast T\left(x,t\right)$$

(14)

The reformulated LBM with time-averaged can be expressed as follows:

$$f_1\left(x+1, t+1\right)=a{f}_1\left(x, t\right)+b{\displaystyle \frac{f_1\left(x, t+1\right)+f_1\left(x, t-1\right)}{2}}$$

(15)

$$f_2\left(x-1, t+1\right)=a{f}_2\left(x, t\right)+b{\displaystyle \frac{f_2\left(x, t+1\right)+f_2\left(x, t-1\right)}{2}}$$

(16)

The effective values of $a and b$ in the reformulated LBM cannot be determined analytically; instead, they serve as numerically determined relaxation parameters that governs the rate of convergence, similar in role to the conventional relaxation time in the standard BGK model. To reflect the structural similarity of the reformulated approach with FDM and LBM, we assume

$$b=1-a$$

(17)

Since future values of the distribution functions ${f_1\left(x, t+1\right), f}_2\left(x, t+1\right)$ are not known in advance, we first use Eq. (10) and Eq. (11) to calculate the future values, then, we substitute these values into Eq. (15) and Eq. (16) to rebuild the updated distribution functions $f_1\left(x+1, t+1\right),f_2\left(x-1, t+1\right)$, this sequence is illustrated in the computational flowchart Figure 1.

The D1Q2 lattice is selected in this study due to its minimal complexity, which allows for isolating and evaluating the effect of the proposed time- averaged equilibrium strategy without the influence of higher-order lattice dynamics. The simulation domain is initialized with a uniform temperature of 0 (in lattice unit). Dirichlet boundary conditions are applied at both ends of the domain: the left boundary is fixed at a temperature of 1, while the right boundary is fixed at a temperature of 0 unity. The thermal diffusivity $k=0.25$, the length of bar is $1 m$ and $\Delta x=0.01$ for FDM and 100 lattice nodes for LBM.

2.6. Mathematical Formulation and Finite Volume Discretization of the Sod Shock Tube

The Sod shock tube problem is a classical benchmark for evaluating the performance of numerical methods in compressible flows. It is governed by the Euler equations, which explain the conservation of mass, momentum, and energy in an inviscid fluid. In this study, gravitational force and other body forces are neglected, as the flow is assumed to be driven purely by pressure discontinuity.

For the two-dimensional case, the Euler equations can be written in conservative form as:

$${\displaystyle \frac{\partial U}{\partial t}}+{\displaystyle \frac{\partial F_x}{\partial x}}+{\displaystyle \frac{\partial F_y}{\partial y}}=0$$

(18)

where the conserved variable vector$U$, and the flux vectors $F_x$ and $F_y$, are:

$$U=\left(\begin{array}{c}\rho \\ \begin{array}{c}\rho u_x\\ \rho u_y\end{array}\\ E\end{array}\right), F_x=\left(\begin{array}{c}\rho u_x\\ \begin{array}{c}{{\rho u}_x}^2+p\\ {\rho u}_x{u}_y\end{array}\\ \left(E+P\right)u_x\end{array}\right) ,{ F}_y=\left(\begin{array}{c}\rho u_y\\ \begin{array}{c}{\rho u}_x{u}_y\\ {{\rho u}_y}^2+p\end{array}\\ \left(E+P\right)u_y\end{array}\right)$$

(19)

The system consists of four variables: $\rho$ represents the density, $u_x$ is the velocity in the x direction, $u_y$ represents the velocity in y direction, $E$ represents the energy, $p$ represents the pressure. to make it solvable, an additional relation is needed, this is provided by the equation of state.

$$E={\displaystyle \frac{P}{\gamma -1}}+{\displaystyle \frac{{\rho u}^2}{2}}$$

(20)

where $\gamma$ represents the specific heat ratio.

The standard FVM discretization of this equation for cell $x,y$ is given by:

$$U_{x,y}^{n+1}=U_{x,y}^n-{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x+{\displaystyle \frac{1}{2}}}^n-F_{x-{\displaystyle \frac{1}{2}}}^n\right)\mathrm{br}-\mathrm{to}-\mathrm{break} -{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{y+{\displaystyle \frac{1}{2}}}^n-F_{y-{\displaystyle \frac{1}{2}}}^n\right)$$

(21)

Where n represents the time step, in this study, the numerical flux at cell interface is evaluated using HLL approximate Riemann solver, combined with a MUSCL scheme, as proposed in previous studies [14,15], the resulting flux at the cell interface is expressed as:

$$F_{ i+{\displaystyle \frac{1}{2}}}^n\left(U_{HLL}\right)=HLL(U_{i+{\displaystyle \frac{1}{2}}}^L,U_{i+{\displaystyle \frac{1}{2}}}^R)$$

(22)

where $U_{i+{\displaystyle \frac{1}{2}}}^L$ and $U_{i+{\displaystyle \frac{1}{2}}}^R$ are reconstructed left and right states at the interface and $i$ is the centre of the cell, obtained by MUSCL scheme. Eq. (21) can be approximated using the second-order Runge-Kutta method. First, the slope is calculated at the first-time step using the conservative variable values. This slope is then used to extrapolate the conservative variables forward by one time step. Using the extrapolated values, a second slope is computed at the new conservative variable values, as follows:

$$K_1^n=-{\displaystyle \frac{1}{\Delta x}}\left(F_{x +{\displaystyle \frac{1}{2}}}^n \left(U_{HLL}\right)-F_{x -{\displaystyle \frac{1}{2}}}^n\left(U_{HLL}\right)\right)-{\displaystyle \frac{1}{\Delta y}}\left(F_{y +{\displaystyle \frac{1}{2}}}^n \left(U_{HLL}\right)-F_{y -{\displaystyle \frac{1}{2}}}^n\left(U_{HLL}\right)\right)$$

(23)

$$U_{x,y}^{n+1}=U_{x,y}^n+\Delta t{K}_1^n$$

(24)

Eq. (23) represent the first order Runge-Kutta method, where $K_1^0$ represents the slope at the initial time, this equivalent to the explicit Euler method and serves as a basic time integration scheme, to extrapolate the conservative variable at next time step:

$$U_{x,y}^{n+1}=U_{x,y}^n+\Delta t{K}_1^n$$

(25)

The slope at next time step can be obtained by:

$$\begin{gathered} K_2^{n+1}=-{\displaystyle \frac{1}{\Delta x}}\left(F_{x+{\displaystyle \frac{1}{2}}}^{n+1} \left({U_{x,y}^{n+1}}_{HLL}\right)-F_{x-{\displaystyle \frac{1}{2}}}^{n+1}\left({U_{x,y}^{n+1}}_{HLL}\right)\right) \\ -{\displaystyle \frac{1}{\Delta y}}\left(F_{y+{\displaystyle \frac{1}{2}}}^{n+1} \left({U_{x,y}^{n+1}}_{HLL}\right)-F_{y i-{\displaystyle \frac{1}{2}}}^{n+1}\left({U_{x,y}^{n+1}}_{HLL}\right)\right) \end{gathered}$$

(26)

Reupdate the conservative variable by averaging the slopes as follows:

$$U_{x,y}^{n+1}=U_{x,y}^n+\Delta t{\displaystyle \frac{K_1^n+K_2^{n+1}}{2}}$$

(27)

Eq. (26) represents the second order Runge-Kutta method, improving the accuracy of the time integration. The second order scheme along with the first order will be used as a base line for comparison, in this study the result from the four approaches will be compared: the analytical solution, the first order finite volume method, the second order Runge-Kutta method and the proposed time-averaged BGK.

2.7. Transformation of FVM-HLL-MUSCL into a BGK-like Relaxation Form

To enable a BGK-style relaxation approach in the context of FVM, we start from the classical 2D explicit finite volume update for the conserved variables. The update equation for a control volume at $x+1$is given by:

$$U_{x+1,y}^{n+1}=U_{x+1,y}^n-{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x +1+{\displaystyle \frac{1}{2}},y}^n-F_{x+1-{\displaystyle \frac{1}{2}},y}^n\right)-{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x+1,y+{\displaystyle \frac{1}{2}}}^n-F_{x+1,y-{\displaystyle \frac{1}{2}}}^n\right)$$

(28)

Sum Eq. (28) with Eq. (21), rearrange it in term of $U_{x+1,y}^{n+1}$:

$$\begin{gathered} U_{x+1,y}^{n+1}=U_{x,y}^n+ U_{x+1,y}^n- U_{x,y}^{n+1} -{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x+1+{\displaystyle \frac{1}{2}},y}^n-F_{x-{\displaystyle \frac{1}{2}},y}^n\right)-{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x+1,y+{\displaystyle \frac{1}{2}}}^n-F_{x+1,y-{\displaystyle \frac{1}{2}}}^n \right) \\ -{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x,y+{\displaystyle \frac{1}{2}}}^n-F_{x,y-{\displaystyle \frac{1}{2}}}^n\right) \end{gathered}$$

(29)

Assume:

$$\begin{gathered} U_{x,y}^{eq1}\left(x,y, t\right)= U_{x+1,y}^n- U_{x,y}^{n+1} -{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x+1+{\displaystyle \frac{1}{2}},y}^n-F_{x-{\displaystyle \frac{1}{2}},y}^n\right)-{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x+1,y+{\displaystyle \frac{1}{2}}}^n-F_{x+1,y-{\displaystyle \frac{1}{2}}}^n \right) \\ -{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x,y+{\displaystyle \frac{1}{2}}}^n-F_{x,y-{\displaystyle \frac{1}{2}}}^n\right) \end{gathered}$$

(30)

Leads to:

$$U_{x+1,y}^{n+1}=U_{x,y}^n+ U_{x,y}^{eq1}\left(x, t\right)$$

(31)

We can repeat the same procedure for $x-1$, the resultant equations are:

$$\begin{gathered} U_{x-1,y}^{n+1}=U_{x,y}^n+U_{x-1,y}^n- U_{x,y}^{n+1}-{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x+{\displaystyle \frac{1}{2}},y}^n-F_{x-1-{\displaystyle \frac{1}{2}},y}^n\right) -{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x-1,y+{\displaystyle \frac{1}{2}}}^n-F_{x-1,y-{\displaystyle \frac{1}{2}}}^n \right) \\ -{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x,y+{\displaystyle \frac{1}{2}}}^n-F_{x,y-{\displaystyle \frac{1}{2}}}^n\right) \end{gathered}$$

(32)

$$\begin{gathered} U_{x,y}^{eq2}\left(x, y,t\right)=U_{x-1,y}^n-U_{x,y}^{n+1}-{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x+{\displaystyle \frac{1}{2}},y}^n-F_{x-1-{\displaystyle \frac{1}{2}},y}^n\right)-{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x-1,y+{\displaystyle \frac{1}{2}}}^n-F_{x-1,y-{\displaystyle \frac{1}{2}}}^n \right) \\ -{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x,y+{\displaystyle \frac{1}{2}}}^n-F_{x,y-{\displaystyle \frac{1}{2}}}^n\right) \end{gathered}$$

(33)

$$U_{x-1,y}^{n+1}=U_x^n+U_x^{eq2}\left(x,y, t\right)$$

(34)

The same procedure for $y$:

$$U_{x,y}^{eq3}\left(x,y, t\right)= U_{x,y+1}^n- U_{x,y}^{n+1} -{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x,y+1+{\displaystyle \frac{1}{2}}}^n+F_{x,y-{\displaystyle \frac{1}{2}}}^n\right)-{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x+{\displaystyle \frac{1}{2}},y+1}^n-F_{x-{\displaystyle \frac{1}{2}},y+1}^n \right)-{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x+{\displaystyle \frac{1}{2}},y}^n-F_{x-{\displaystyle \frac{1}{2}},y}^n\right)$$

(35)

$$U_{x,y+1}^{n+1}=U_{x,y}^n+U_{x,y}^{eq3}\left(x,y, t\right)$$

(36)

$$\begin{gathered} U_{x,y}^{eq4}\left(x,y, t\right)=U_{x,y-1}^n-U_{x,y}^{n+1}-{\displaystyle \frac{\Delta t}{\Delta y}}\left(F_{x,y+{\displaystyle \frac{1}{2}}}^n-F_{x,y-1-{\displaystyle \frac{1}{2}}}^n\right)-{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x+{\displaystyle \frac{1}{2}},y-1}^n-F_{x-{\displaystyle \frac{1}{2}},y-1}^n \right) \\ -{\displaystyle \frac{\Delta t}{\Delta x}}\left(F_{x+{\displaystyle \frac{1}{2}},y}^n-F_{x-{\displaystyle \frac{1}{2}},y}^n\right) \end{gathered}$$

(37)

$$U_{x,y-1}^{n+1}=U_{x,y}^n+U_{x,y}^{eq4}\left(x,y, t\right)$$

(38)

To construct a BGK-like analogy D2Q4, we define pseudo distribution functions:

$$\begin{gathered} f_1^{eq}\left(x,y, t\right)=0.25\ast U_{x,y,t}^{eq1},\hspace{1em}\hspace{1em}{f}_2^{eq}\left(x,y, t\right)=0.25\ast U_{x,y,t}^{eq2} \\ f_3^{eq}\left(x,y, t\right)=0.25\ast U_{x,y,t}^{eq3},{ f}_4^{eq}\left(x, y,t\right)=0.25\ast U_{x,y,t}^{eq4} \end{gathered}$$

(39)

And:

$$U_{x,y}^n=f_1\left(x,y,t\right)+ f_2\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)+ f_3\left(x,y,t\right)+ f_4\left(\mathrm{x},\mathrm{y},\mathrm{t}\right)$$

(40)

The transformation of the LBM with time-averaged formulation can be expressed as follows:

$$f_1\left(x+1,y, t+1\right)=a{f}_1\left(x,y, t\right)+b{f}_1^{eq}\left(x,y, t\right)$$

(41)

$$f_2\left(x-1,y, t+1\right)=a{f}_2\left(x,y, t\right)+b{f}_2^{eq}\left(x,y ,t\right)$$

(42)

$$f_3\left(x,y+1, t+1\right)=a{f}_3\left(x, y,t\right)+b{f}_3^{eq}\left(x,y, t\right)$$

(43)

$$f_4\left(x,y-1, t+1\right)=a{f}_4\left(x, y,t\right)+b{f}_4^{eq}\left(x,y, t\right)$$

(44)

where $a$ and $b$ are numerical coefficients. In this study, we set $a=1$ and $b=1$, to maintain structural consistency. These reconstructed forms of the update equations enable the application of the time-averaged equilibrium formulation, as previously discussed.

2.8. Parallelization Compatibility

The reformulated equations, expressed in a BGK-like relaxation form with a time-averaged equilibrium function, is particularly well-suited for parallel implementation. Since the update at each lattice node relies only on its immediate neighbours and the equilibrium distribution, the method permits efficient domain decomposition and parallel execution. This locality makes it highly compatible with MPI for distributed memory systems and CUDA for GPU acceleration, enabling scalable performance on high-performance computing platforms.

2.9. Sod Shock Tube Configuration and Analytical Reference

The Sod shock tube is initialized with a partition at the centre of the domain. The left state is defined as $\mathsf{\rho }\mathrm{L}=1.0,\mathrm{u}\mathrm{L}=0,\mathrm{p}\mathrm{L}=1$ and the right state is defined as $\mathsf{\rho }\mathrm{R}=0.125,\mathrm{u}\mathrm{R}=0,\mathrm{p}\mathrm{R}=0.1$, the ratio of the specific heat is set to $\gamma =1.4$, where R and L represent right and left respectively.

The analytical solution for the sod shock tube is available from Riemann problem for the Euler equations, as described by Toro [16], this solution is used as a benchmark to validate the numerical results presented in this study.

3. Results

3.1. Heat Diffusion

To evaluate the performance of the proposed time-averaged equilibrium LBM (TAE-LBM) on heat diffusion equation, two benchmark simulations were used and compared with the FDM, the focus is on two parameters: the number of time steps needed to simulate a fixed physical time and the solution accuracy, calculated by root mean square error RMSE. Additionally, observations from additional runs using extreme values of relaxation coefficient $a$ to assess the stability of the TAE-LBM under varying relaxation conditions.

3.1.1. Simulation 1

In the first simulation, the temperature distribution was observed over a short physical duration, as shown in Table 1 and Figure 2, FDM required a 400 time-steps with $\mathsf{\Delta }\mathrm{t}=0.5$ while the LBM $\mathsf{\Delta }\mathrm{t}=1$ completed the same simulation in 200 time-steps.

TAE-LBM, calculated for two coefficients $a=0.75$ and $a=1$, needed 35 time-steps and 45 time-steps, respectively, the standard LBM produced the lowest RMSE, indicating the highest accuracy among all methods, TAE-LBM with $a=0.75$ produced the highest RMSE, while LBM with $a=1$ showed better accuracy than LBM with $a=0.75$, indicating that the accuracy improve as $a$ increases.

3.1.2. Simulation 2

In the second simulation, the temperature distribution was observed over a longer physical duration, as shown in Table 1 and Figure 2, FDM required 800 time-steps with $\mathsf{\Delta }\mathrm{t}=0.5$ and LBM $\mathsf{\Delta }\mathrm{t}=1$ required 400 time-steps, TAE-LBM performed consistently, finished the simulation in 70 time-steps for $a=0.75$ and 90 for $a=1$.

Accuracy stayed similar, LBM produced the lowest RMSE, followed by TAE- LBM with $a=1$ and $a=0.75$, again, the accuracy improved as $a$ increased.

Table 2. Convergence and accuracy in simulation 2.

Simulation 2	FDM dt=0.5	LBM	TAE-LBM a=0.75	TAE-LBM a=1
Number of time steps	800	400	70	90
RMSE		0.00031	0.006159	0.00303

Table 2. Convergence and accuracy in simulation 2.

Simulation 2	FDM dt=0.5	LBM	TAE-LBM a=0.75	TAE-LBM a=1
Number of time steps	800	400	70	90
RMSE		0.00031	0.006159	0.00303

3.1.3. Visual Analysis

Figure 2 presents the temperature distribution results from Simulations 1 and 2. In both simulations, the TAE-LBM curves closely follow the FDM and LBM results. The TAE-LBM with a=1 consistently gives better convergence and lower RMSE than TAE-LBM with a=0.75, as it reduces to a purely time-averaged equilibrium that depends only on the equilibrium distribution function, both simulations, however, use fewer time-steps than LBM and FDM, lowering the value of a results in faster convergence, but at the cost of accuracy.

Figure 2. Comparison of spatial temperature distributions using LBM, FDM, and time-averaged equilibrium method with its proposed coefficient (a) are compared at two different simulation times.

Figure 2. Comparison of spatial temperature distributions using LBM, FDM, and time-averaged equilibrium method with its proposed coefficient (a) are compared at two different simulation times.

3.1.4. Extreme Coefficient Values

Additional test with extreme values of $a$ showed limitations of TAE-LBM at $a=0$, the solution remains static, showing no diffusion. At $a=1.22$ the solution becomes unstable with non-physical values. these observations lead to the definition of the valid domain as 1.2<a>0.

3.2. Sod Shock Tube

The Sod shock tube problem was simulated to evaluate the accuracy of different numerical methods at t=0.2 s. Figure 3 presents the density distribution for three different numerical methods FVM-HLL-MUSCL scheme with first order Runge-Kutta, FVM-HLL-MUCSLE with second order Runge-Kutta and TAE-BGK-like relaxation form. These are compared against the analytical solution. The first-order FVM produces noticeable diffusion and oscillations near the shock discontinuity. In contrast, the second-order FVM and the BGK-like relaxation form solved by TAE-LBM capture the discontinuity with improved accuracy.

Table 1 lists the RMSE for each method. The TAE-BGK-like relaxation form achieved the best result with the analytical solution.

Table 3. RMSE of density for different methods.

Simulation 3	Analytical Solution	FVM first order	FVM second order	BGK-like relaxation form
RMSE		0.02809	0.01639	0.015134

Table 3. RMSE of density for different methods.

Simulation 3	Analytical Solution	FVM first order	FVM second order	BGK-like relaxation form
RMSE		0.02809	0.01639	0.015134

In addition to the 1D results, a 2D simulation was produced using the TAE-LBM method:

Figure 4. Figure 4 shows the density profile comparison at $t=0.2 s$ for the Sod shock tube problem in 2D.

Figure 4. Figure 4 shows the density profile comparison at $t=0.2 s$ for the Sod shock tube problem in 2D.

4. Conclusions

A new method, the Time-Averaged Equilibrium Lattice Boltzmann Method (TAE-LBM), was introduced and investigated by two benchmark simulations: the heat diffusion and Sod shock tube problem. In heat diffusion TAE-LBM achieved faster convergence with significantly fewer number of time steps compared with LBM and FDM, but with less accuracy than other methods, increasing the value of a leads to increased accuracy and a decrease in the number of time steps to converge, a trade-off between stability and precision, while the extreme value of a highlights the method’s limitations, determining a valid operational range. In Sod shock tube, when TAE-LBM was applied on the BGK-like relaxation form, it showed a good agreement in terms of accuracy, no reduction or increase in the number of time steps. Overall, TAE-LBM presents an efficient direction for efficient and accurate numerical simulations, with potential for further optimization and extension.