Computational Performance of Disparate Lattice Boltzmann Scenarios under Unsteady Thermal Convection Flow and Heat Transfer Simulation

: The present work highlights the capacity of disparate lattice Boltzmann strategies in simulating natural convection and heat transfer phenomena during the unsteady period of the ﬂow. Within the framework of Bhatnagar-Gross-Krook collision operator, diverse lattice Boltzmann schemes emerged from two different embodiments of discrete Boltzmann expression and three distinct forcing models. Subsequently, computational performance of disparate lattice Boltzmann strategies was tested upon two different thermo-hydrodynamics conﬁgurations, namely the natural convection in a differentially-heated cavity and the Rayleigh-B è nard convection. For the purposes of exhibition and validation, the steady-state conditions of both physical systems were compared with the established numerical results from the classical computational techniques. Excellent agreements were observed for both thermo-hydrodynamics cases. Numerical results of both physical systems demonstrate the existence of considerable discrepancy in the computational characteristics of different lattice Boltzmann strategies during the unsteady period of the simulation. The corresponding disparity diminished gradually as the simulation proceeded towards a steady-state condition, where the computational proﬁles became almost equivalent. Variation in the discrete lattice Boltzmann expressions was identiﬁed as the primary factor that engenders the prevailed heterogeneity in the computational behaviour. Meanwhile, the contribution of distinct forcing models to the emergence of such diversity was found to be inconsequential. The ﬁndings of the present study contribute to the ventures to alleviate contemporary issues regarding proper selection of lattice Boltzmann schemes in modelling ﬂuid ﬂow and heat transfer phenomena.


Introduction
The Lattice Boltzmann Method (LBM) has raised considerable interest amongst the community of computational fluid dynamics (CFD) due to its efficacy in handling multitude fluid flow problems. Relying on the statistical mechanics, LBM regards the flowing materials through the collective behaviour of the accompanying molecules [1,2]. It differs from traditional CFD methods in sense that it solves the representative Boltzmann expression of the flowing substances rather than directly handling the corresponding hydrodynamics equations in their operations. Because of this unique feature, modelling a flow problem using LBM has several advantages, including a clear algorithm [3], straightforward treatment of boundary conditions [4], and innate feasibility for parallel computing architecture [5]. As a promising numerical tool, LBM is currently a vibrant research topic in the discipline of CFD.
As far as modelling the flow problem using LBM is concerned, two properties are routinely considered: (a) the discrete lattice Boltzmann expressions, and (b) the discrete 2 of 21 forcing schemes. The former parameter has been extensively discussed in the literature. Ubertini et al. [6] investigated three distinct models of discrete lattice Boltzmann expression for hydrodynamics simulation, namely the first-order, the second-order, and the scheme derived through implementation of the Verlet discretisation, which all showed a secondorder accuracy both in spatial and time coordinates with respect to the convective system. They argued that such equivalence breaks down when the nature of the physical systems necessitates the inclusion of external forcing expression.
Silva and Semiao [7] carried out a comprehensive assessment of distinct lattice Boltzmann remarks using Chapman-Enskog analysis. Guo et al. [8] highlight the significance of the different lattice Boltzmann schemes towards the accuracy of the recovered Navier-Stokes expression from the Chapman-Enskog analysis. Additionally, they mentioned that the choice of discrete forcing model depends heavily upon the exactitude of the restored continuum hydrodynamics representation.
On the other hand, the discrete forcing model in LBM is the other prominent property when modelling the flow is sought. Several authors have proposed diverse expressions to accommodate external forcing term in the generic lattice Boltzmann equation. Luo [9] was among the first authors to propose a mathematical expression for the discrete forcing term in LBM, alongside with He et al. [10] and Guo et al. [8,11]. Later on, Kupershtokh et al. [12] introduced a forcing model based on the association of exact-difference-method upon the corresponding Boltzmann equation.
The presence of different mathematical expressions for both the discrete lattice Boltzmann equation and the forcing model offers diverse LBM strategies for hydrodynamics modelling. However, choosing a suitable approach is still a matter of debate. To alleviate such issues, Mohamad and Kuzmin [13] investigated the behaviour of three different forcing models by simulating natural convection in closed and open-ended cavities. They found that the investigated forcing models produced equivalent numerical solutions at steady-state conditions. Subsequently, Krivovichev [14] presented a comprehensive evaluation regarding stability analysis of the six widely-used forcing models based on the application of the von Neumann method to linear approximation of the system of nonlinear lattice Boltzmann expressions. They found that better stability properties prevailed upon the forcing models that are implicit in their nature. Zheng et al. [15] found that as long as the simulation comprises low Mach number flow, different forcing schemes in LBM would return equivalent steady-state solutions.
It is apparent from the available literature that few works have looked into the implications of both distinct lattice Boltzmann expressions and different forcing models in LBM simulation. Moreover, the primary focus of the published works has revolved around the discrepancy in the computational characteristics of distinct LBM scenarios during the steady-state condition of the simulation, leaving the behaviour during the unsteady period of the flow unexplained. Therefore, the present work aims to expand the evaluation of the capacity of disparate LBM schemes by incorporating the variety in both the discrete lattice Boltzmann expressions and the discrete forcing models into the workflow of the assessment. Emphasis was given upon revealing the computational characteristic of distinct LBM scenarios during the unsteady-period of the simulation. The steady-state solutions were considered for the exhibition and validation purposes.
In this study, two different embodiments of discrete Boltzmann expression and three distinct forcing models were incorporated into the analysis. For the sake of providing comprehensive evaluation, the capacity of disparate LBM schemes was tested upon simulating two non-isothermal thermo-hydrodynamics systems, namely natural convection in a differentially-heated cavity and Rayleigh-Bènard convection.
The current research article is organized as follows. Section 2 presents the foundational aspect of natural convection and heat transfer arrangement. Section 3 describes the distinct LBM strategies for simulating the concomitant physical phenomena. Section 4 deals with Computation 2021, 9, 65 3 of 21 theoretical outlook of disparate LBM scenarios, while the results of numerical testing are discussed in Section 5. Some concluding remarks are posed in Section 6.

Governing Mathematical Remarks and Principal Dimensionless Groups
Natural convection phenomena is governed by two flowing materials that propagate simultaneously and synergistically within the flow domain, namely the fluid and thermal substances. Accordingly, the natural convection is governed by three fundamental equations, namely, the continuity, Navier-Stokes, and heat equations [16], expressed as ∂ρ ∂t ∂T ∂t where ρ, u i , p, T,µ and D denote the fluid density, velocity, pressure, temperature, dynamic viscosity, and thermal diffusion coefficient, respectively. The contribution from the buoyancy effect was designated by F α .
In the present work, the natural convection was examined under the Boussinesq approximation [13,17]. As such, it neglected the contributions to the dynamical behaviour of the flowing substances from viscous heat dissipation and compression caused. Therefore, the buoyancy term occupies the following definition [13]: where G α specifies the gravitational acceleration, β T depicts the thermal expansion coefficient, and T ref represents the assigned reference temperature. It is customary within CFD to express the associated physical quantities using nondimensional units. For natural convection phenomena, the key dimensionless groups were identified as follows: where Nu, Pr, and Ra designate the Nusselt, Prandtl, and Rayleigh number, accordingly. T hot and T cold denote the hot and cold temperature conditions, respectively. Physical quantities of h, c p , and L represent, respectively, the heat transfer coefficient, fluid specific heat capacity, and the characteristic length of the domain. Furthermore, the present work uses the following generic dimensionless expressions to designate the horizontal and vertical length of the flow domain.
In Equation (6), x and y express, respectively, the horizontal and vertical length of the investigated domain.

Lattice Boltzmann Expressions from Distinct Truncation Order
In the present work, consideration was given to the double-distribution-function (DDF) approach to establish a proper representation of the dynamical entities [17]. This would mean that the prevalence of two continuous Boltzmann expressions, which should be discretised either by capturing only the first-truncation order or considering up to the second-truncation term of the expanded-Boltzmann equation. Regardless the approach considered, different representations of the discrete lattice Boltzmann equation are obtained. For the fluid component, the first-and second-order discrete lattice Bhatnagar-Gross-Krook (BGK) equations are expressed as follows: where f i indicates the discrete fluid populations, ξ iα specifies the discrete velocity of fluid particles, R i designates the discrete form of the external force term, τ υ denotes the relaxation time for fluid particles, and ∆t specifies the discrete time step. For the second-order lattice BGK model, the corresponding discrete fluid population f i and its corresponding relaxation time τ υ were defined accordingly as in [6,10]: where f eq i depicts the equilibrium condition of the fluid population. The equilibrium distribution function f eq i occupies the following remark: where c s designates the lattice speed of sound and w i represents the discrete weighting coefficient.
In the present work, the D2Q9 velocity arrangement was selected to model fluid displacements. Such configuration satisfies the following descriptions: Apart from the dichotomy in the lattice BGK representation for fluid movements, the discrete manifestation of forcing term R i also occupies diverse expressions. In this study, three of the most prominent R i models were selected for the evaluation, namely the scheme proposed by Luo [9], Guo et al. [8], and Kupershtokh et al. [12], expressed respectively as follows: For the thermal constituent, due to the absence of external force term, the discretisation operation of Boltzmann representation of thermal substance is not affected by the selected truncation order. Therefore, the lattice BGK formula representing the thermal dissemination only occupies a single mathematical expression, namely, where g i identifies the discrete thermal population, e iα denotes the discrete velocity of thermal particles, and τ D designates the relaxation time for thermal elements. The associated thermal equilibrium population g eq i was defined as follows: For thermal dissemination, D2Q5 discrete velocity set was adopted. Implementation of fewer velocity directions to represent thermal displacements was possible due to the lower isotropy nature pertinent to thermal populations in the corresponding lattice Boltzmann arrangement. As can be seen in Appendix B, the associated Chapman-Enskog analysis of thermal particles only necessitates low-order moment expansion terms in order to recover correct macroscopic heat formulation. Consequently, the isotropy requirement for thermal displacement modelling in LBM can be conducted properly by fewer lattice arrangement, such as the D2Q5 velocity set [2].
For the adopted D2Q5 configuration, the discrete velocity of thermal particles e iα and their corresponding weighting factor z i adhere to the following descriptions: z i = At this point, it is apparent that there exist distinct strategies to administer natural convection and heat transfer modelling with LBM. Table 1 outlines the plausible LBM scenarios considered in the present work, as well as the corresponding scheme. Luo [9] (Equation (14)) IIB Guo et al. [8] (Equation (15)) IIC Kupershtokh et al. [12] (Equation (16))

Recovery of the Macroscopic Thermo-Hydrodynamics Expressions from Disparate LBM Arrangements in Natural Convection Heat Transfer Modelling
As a thermo-hydrodynamics solver, LBM is closely tied to the macroscopic heat and mass transport formulations. It is therefore pivotal to uncover the capacity of each considered LBM schemes in returning the fundamental macroscopic relationships. This objective was fulfilled through evaluating the Chapman-Enskog analysis [7,8,18]. Appendices A and B provide the details of undertaken procedures for the associated fluid and thermal components, respectively.
From the corresponding Chapman-Enskog analysis, the restored macroscopic equations of mass, momentum, and heat transport were captured as follows: ∂ρ ∂t Computation 2021, 9, 65 ∂T ∂t In Equation (22), the term ∑ i ξ iα ξ iβ R i specifies the second-order moment of the forcing term. Details regarding the mathematical expressions of all the associated forcing moments are appended in Appendix A of this manuscript. In Equations (21) and (22), quantities m, σ, and ϕ occupy the following definitions: From the recovered thermo-hydrodynamics expressions of Equations (21)-(23), it was observed that each of the corresponding LBM scenarios returns the macroscopic hydrodynamics relationships with residual fractions. Such properties occupy different mathematical expressions depending upon the selected LBM scenarios. Table 2 summarizes the corresponding mathematical remarks of hydrodynamics residual fractions from distinct LBM schemes. On the other hand, the restored heat formulation of Equation (23) contains only one residual fraction, appropriately depicted by the last term on the right hand side of the formula. Nevertheless, it is important to mention that the particular thermal residual fraction incorporates momentum terms from the fluid constituents within its expression. Table 2. The captured residual fractions in the recovered continuity and Navier-Stokes equations from every considered LBM scenarios.

Residual Fractions in the Restored
Navier-Stokes Equation Consequently, the residual term in the recovered heat equation is linked indirectly to the residual terms in the restored Navier-Stokes formula of Equation (22).

Numerical Results and Discussions
To test the computational performance of disparate LBM scenarios, two distinctive thermo-hydrodynamics flow systems (the natural convection in a differentially-heated cavity and the Rayleigh-Bènard convection) were selected. The simulations of both physical phenomena were performed in a two-dimensional close system under the condition of Ra = 10 4 and Pr = 0.71. The computational workloads were administered upon the graphical processing unit (GPU) ecosystem so as to expedite the analyses. In such physical arrangements, the Mach number (Ma) was defined as where u char designates the characteristic velocity of the flowing material.
The particular property was defined as where Θ hot and Θ cold specify, respectively, the dimensionless hot and cold temperatures.
In general, the dimensionless temperature is expressed as where T ref denotes the reference temperature.
The expressions for fluid kinematic viscosity υ and thermal diffusivity D were obtained, respectively, as    The entire domain was assumed filled with an inert fluid. The vertical walls of the cavity were characterised by contrasting thermal conditions. The left-border was considered hot (T hot ) while the opposite wall was assumed cold (T cold ). On the other hand, the horizontal boundaries were set to be perfectly insulated. For the fluid substance, the non-equilibrium bounce-back (NEBB) technique [19] was adopted upon all the corresponding perimeters, including the four corners. The NEBB method defines the ered hot (T hot ) while the opposite wall was assumed cold (T cold ). On the other hand, the horizontal boundaries were set to be perfectly insulated. For the fluid substance, the nonequilibrium bounce-back (NEBB) technique [19] was adopted upon all the corresponding perimeters, including the four corners. The NEBB method defines the unknown fluid populations at the boundaries from the known populations by evaluating the non-equilibrium condition of the normal populations at local boundary sites.

Simulation of Natural Convection in a Differentially-Heated Cavity
The non-equilibrium condition at the boundary wall satisfies the following remark: where N α notifies the momentum correction factor due to inclusion of the external forcing term into the generic Boltzmann expression of the fluid populations. The subscript symbol i specifies the opposite direction of i. Parameter ζ i was defined as follows: Meanwhile, for the thermal population, the anti-bounce-back (ABB) [2] and central finite difference [20] strategies were administered upon the vertical and horizontal boundaries, respectively. ABB technique was administered by adopting the following relationship for the associated boundary walls: where g * i depicts thermal populations after collision and T wall notifies the prescribed temperature condition at particular boundary wall. To accomplish a valid comparison, numerical simulations for all the considered LBM scenarios were accomplished under similar hydrodynamic relaxation time (τ υ = 0.6) and Mach number (Ma = 0.1) conditions. Thereupon, the average Nusselt number Nu was appointed as the dimensionless property that represents the heat transfer performance of each concomitant LBM scheme. This property abides to the following description: where u * x Θ represents the average of the product between the dimensionless horizontal velocity and the dimensionless temperature over the entire simulation domain. For the sake of validation, the associated steady-state flow profile from scenario IIB was selected as the representation to exhibit the final streamlines and isotherms from the corresponding heat and mass transport system. Figure 2 exhibits the steady-state streamlines and isotherms from the selected LBM scenario. The corresponding streamlines and isotherms profiles demonstrate an excellent agreement with the literature [16,17,21,22]. It was found, however, that the results from the six associated LBM scenarios exhibit almost similar steady-state flow profiles. Therefore, every considered LBM scheme possesses similar qualification in delivering legitimate steady-state numerical solutions. corresponding heat and mass transport system. Figure 2 exhibits the steady-state streamlines and isotherms from the selected LBM scenario. The corresponding streamlines and isotherms profiles demonstrate an excellent agreement with the literature [16,17,21,22]. It was found, however, that the results from the six associated LBM scenarios exhibit almost similar steady-state flow profiles. Therefore, every considered LBM scheme possesses similar qualification in delivering legitimate steady-state numerical solutions.  Table 3 outlines the principal steady-state numerical properties from distinct LBM schemes. The parameters were compared with the outcomes of finite difference method (FDM) [23] and finite element method (FEM) [24]. A closer look on Table 3 reveals that LBM scenarios, which adopt a second-order lattice BGK model (i.e., scenario IIA, IIB, and IIC) were capable of delivering better compliance with the benchmark solutions than the ones which comprise first-order lattice BGK scheme (i.e., scenario IA, IB, and IC).
Undoubtedly, the inclusion of higher-truncation order terms empowers the secondorder lattice BGK scheme with superior numerical accuracy to the corresponding firstorder counterpart. Contrastingly, trivial disparity in computing performance was observed upon implementation of different forcing models.  Table 3 outlines the principal steady-state numerical properties from distinct LBM schemes. The parameters were compared with the outcomes of finite difference method (FDM) [23] and finite element method (FEM) [24]. A closer look on Table 3 reveals that LBM scenarios, which adopt a second-order lattice BGK model (i.e., scenario IIA, IIB, and IIC) were capable of delivering better compliance with the benchmark solutions than the ones which comprise first-order lattice BGK scheme (i.e., scenario IA, IB, and IC). Table 3. Steady-state properties of two-dimensional natural convection in a differentially-heated cavity for Ra = 10 4 , Pr = 0.71, τ υ = 0.6, and Ma = 0.1, compared with solutions of finite difference method (FDM) [23] and finite element method (FEM) [24]. Undoubtedly, the inclusion of higher-truncation order terms empowers the secondorder lattice BGK scheme with superior numerical accuracy to the corresponding first-order counterpart. Contrastingly, trivial disparity in computing performance was observed upon implementation of different forcing models.

Simulation
Recording the behaviour of specific parameters along the simulation process provides a way of uncovering additional key information regarding the computational capacity of the considered LBM schemes. Figure 3a shows the profiles of Nu from different LBM scenarios with the dimensional simulation time t appointed as the horizontal axis. Figure 3b illustrates similar profiles with dimensionless simulation time t * used as the corresponding horizontal axis. Both figures display the existence of striking contrast in the computational characteristics among distinct LBM scenarios.
Recording the behaviour of specific parameters along the simulation process provides a way of uncovering additional key information regarding the computational capacity of the considered LBM schemes. Figure 3a shows the profiles of ⟨Nu⟩ from different LBM scenarios with the dimensional simulation time t appointed as the horizontal axis. Figure 3b illustrates similar profiles with dimensionless simulation time t * used as the corresponding horizontal axis. Both figures display the existence of striking contrast in the computational characteristics among distinct LBM scenarios. As depicted in Figure 3a, the discrepancy was predominantly apparent when the simulations were performed in an unsteady-state period fashion. Altering the fashion from unsteady to steady state, the disparity either decreases gradually or was negligible (Figure 3a). Such discrepancy was in agreement with the physical properties outlined Table 3. Figure 3a unveils essential information regarding the primary factors responsible for the observed discrepancy in computational performance of distinct LBM scenarios. It As depicted in Figure 3a, the discrepancy was predominantly apparent when the simulations were performed in an unsteady-state period fashion. Altering the fashion from unsteady to steady state, the disparity either decreases gradually or was negligible (Figure 3a). Such discrepancy was in agreement with the physical properties outlined Table 3. Figure 3a unveils essential information regarding the primary factors responsible for the observed discrepancy in computational performance of distinct LBM scenarios. It appears then that the different order of lattice BGK expression is the predominant factor that generates the observed disparity in computational characteristics of different LBM schemes. On the other hand, the contribution of distinct forcing strategies upon such disparity was found to be inconsequential. Figure 3a reveals further that LBM schemes, which administer second-order lattice BGK model (i.e., scenarios IIA, IIB, and IIC) exhibited slower progression towards steady-state condition than those schemes which implement first-order lattice BGK model (i.e., scenarios IA, IB, and IC).
However, when non-dimensional time t * was assimilated, the slower progression characteristic of the second-order schemes vanished. Figure 3b shows that the computational performance of the second-order lattice BGK schemes is proportional to the first-order scenarios. The profiles of dimensionless horizontal velocity u * x at the vertical mid-plane of the enclosure (x * = 0.5) were displayed in Figure 4. Therein, the profiles of u * x demonstrate similar behaviour with the ones observed in Figure 3a.
the first-order scenarios. The profiles of dimensionless horizontal velocity u x * at the vertical mid-plane of the enclosure x * = 0.5 were displayed in Figure 4. Therein, the profiles of u x * demonstrate similar behaviour with the ones observed in Figure 3a. Figure 5 displays the corresponding computational overhead from every considered LBM scenario. The scenario IIC was identified as the particular LBM scheme with the highest computational demand in modelling fluid flow and heat transfer in a differentially-heated cavity, therein.

Simulation of Rayleigh-Bènard Convection
After undertaking evaluation regarding computational characteristics of distinct LBM scenarios upon simulating natural convection in a differentially-heated enclosure in the previous section, the current segment of the article aims at elucidating the capacity of disparate LBM scenarios while simulating the Rayleigh-Bènard convection (RBC) phenomena as illustrated in Figure 6.

Simulation of Rayleigh-Bènard Convection
After undertaking evaluation regarding computational characteristics of distinct LBM scenarios upon simulating natural convection in a differentially-heated enclosure in the previous section, the current segment of the article aims at elucidating the capacity of disparate LBM scenarios while simulating the Rayleigh-Bènard convection (RBC) phenomena as illustrated in Figure 6.  Hot temperature conditions were imposed upon the horizontal bottom wall while the opposing top margin was set to occupy cold temperature. The vertical boundaries were set to be perfectly insulated. Boundary treatments were accomplished through adopting similar strategies with the erstwhile case of natural convection in a differentiallyheated cavity. However, appropriate adjustments were necessary in order to account the appointed wall conditions in RBC configuration.
The contrasting driving force from buoyancy and gravitational attraction in the RBC system results in perpetual competition between the tendency of the flowing materials to move upward and downward, correspondingly. Such a situation enables the associated Hot temperature conditions were imposed upon the horizontal bottom wall while the opposing top margin was set to occupy cold temperature. The vertical boundaries were set to be perfectly insulated. Boundary treatments were accomplished through adopting similar strategies with the erstwhile case of natural convection in a differentially-heated cavity. However, appropriate adjustments were necessary in order to account the appointed wall conditions in RBC configuration.
The contrasting driving force from buoyancy and gravitational attraction in the RBC system results in perpetual competition between the tendency of the flowing materials to move upward and downward, correspondingly. Such a situation enables the associated thermo-hydrodynamics phenomena to exhibit a number of plausibly distinct flow behaviours with variable convection roll patterns [25,26]. To mitigate such complexity, the present study assimilates infinitesimal disturbances into the corresponding physical system. These perturbation functions are described as where N x and N y designate the associated lattice nodes in the horizontal and vertical directions of the spatial coordinate, respectively. The performance of heat transfer in RBC system was represented by the average Nusselt number at the hot wall, Nu 0 . The corresponding parameter was mathematically expressed as Here, q y specifies the local heat flux in vertical direction, defined as The last term on the right-hand side of the above formula denotes temperature gradient in a vertical direction. Figure 7 presents the final streamlines and isotherms for RBC simulation with Ra = 10 4 and Pr = 0.71. Similar to the former case of natural convection in a differentially-heated cavity, the steady-state flow profile from scenario IIB was selected for exhibition and validation purposes. The corresponding streamlines and isotherms displayed in Figure 7 were in excellent agreement with the earlier work of Ouertatani et al. [27]. Similarly to the former case of natural convection in a differentially-heated cavity, the obtained steady-state responses from distinct LBM schemes demonstrate identical flow profiles.
Nevertheless, minor discrepancy was observed in the captured Nu 0 solutions, which are summarized and compared with the outcomes of finite volume method (FVM) [27] in Table 4. Table 4. Average Nusselt number at the hot wall Nu 0 of the Rayleigh-Bènard convection system from distinct LBM schemes during the steady-state period of the flow for aspect ratio = 1, Ra = 10 4 , Pr = 0.71, τ υ = 0.6, and Ma = 0.1 compared with the outcome of the finite volume method (FVM) [27]. A better accuracy of Nu 0 solutions was obtained from the LBM scenarios which adopt a second-order lattice BGK model. The profiles of Nu 0 along the simulation process exhibited characteristics similar to those observed in the former natural convection case. Figure 8a illustrates further the behaviour of Nu 0 . A slow progression characteristic of the second-order lattice BGK schemes was observed that disappears as the horizontal axis is replaced by the associated dimensionless simulation time t * . Figure 9 shows the profiles of dimensionless vertical velocity u * y at the horizontal mid-plane of the cavity (y * = 0.5). Therein, similar computational behaviour was observed as the one prevailing in Figure 8a. gradient in a vertical direction. Figure 7 presents the final streamlines and isotherms for RBC simulation with Ra = 10 4 and Pr = 0.71. Similar to the former case of natural convection in a differentially-heated cavity, the steady-state flow profile from scenario IIB was selected for exhibition and validation purposes. The corresponding streamlines and isotherms displayed in Figure 7 were in excellent agreement with the earlier work of Ouertatani et al. [27]. Similarly to the former case of natural convection in a differentially-heated cavity, the obtained steadystate responses from distinct LBM schemes demonstrate identical flow profiles.

Simulation
Nevertheless, minor discrepancy was observed in the captured ⟨Nu⟩ 0 solutions, which are summarized and compared with the outcomes of finite volume method (FVM) [27] in Table 4.   A better accuracy of ⟨Nu⟩ 0 solutions was obtained from the LBM scenarios which adopt a second-order lattice BGK model. The profiles of ⟨Nu⟩ 0 along the simulation process exhibited characteristics similar to those observed in the former natural convection case. Figure 8a illustrates further the behaviour of ⟨Nu⟩ 0 . A slow progression characteristic of the second-order lattice BGK schemes was observed that disappears as the horizontal axis is replaced by the associated dimensionless simulation time t * . Figure 9 shows the profiles of dimensionless vertical velocity u y * at the horizontal mid-plane of the cavity y * = 0.5 . Therein, similar computational behaviour was observed as the one prevailing in Figure 8a. Figure 10 depicts the corresponding profiles of computational cost from every considered LBM schemes. Similar with the case of natural convection in a differentially-

Conclusions
Comprehensive evaluation regarding the efficacy of disparate Lattice Boltzmann Method (LBM) scenarios upon simulation of fluid flow and heat transfer phenomena were studied. The primary objective, herein addressed, was the evaluation of the plausible discrepancy in the computational characteristics of different LBM scenarios when simulating natural convection and heat transfer systems during the unsteady period of the flow. To fulfil the sought objective, the LBM schemes were tested upon two distinctive thermo-hydrodynamics systems, namely the natural convection in a differentially-heated cavity and the Rayleigh-Bènard convection. The key findings of this work are as follows: 1 The presence of considerable discrepancy in computational characteristics of disparate LBM schemes was seen during the unsteady period of the simulation, which diminished gradually as the simulation advanced towards a steady-state condition. 2 Variation in the associated discrete lattice Boltzmann expression was identified as the predominant factor inherent to discrepancy in computational characteristics. 3 The contribution of distinct forcing models upon the heterogeneity in computational behaviour was found to be trivial. 4 At a steady-state condition, the LBM schemes which administer a second-order lattice BGK model recovered better numerical accuracy than those scenarios which comprise a first-order lattice BGK model. However, the scheme is challenged by higher computational demand. α, β, γ direction of spatial coordinates (Einstein notation).

Appendix A. The Chapman-Enskog Analysis for Fluid Populations
The Chapman-Enskog analysis was commenced by defining the following expanded fractions: In the above relationships, denotes small quantity which value lies within the order of Knudsen number and ζ i satisfies the remark described in Equation (33).
Adopting Taylor series expansion upon the generalized lattice Boltzmann expression for fluid components around equilibrium condition by taking as the expansion quantity returns the mathematical expressions for the first-and second-order expansion terms of the fluid density evolution equation, correspondingly depicted as where O( ) and O 2 specify respectively the first-and second-order expansion terms of the fluid density evolution equation. Parameters σ and ϕ occupy similar definitions as those provided in Equations (25) and (26), respectively.
Thereupon, expressions of moments of O( ) and O 2 can be obtained as follows: Moments of O( ): Moments of O 2 : i .

(A4)
Subsequently, the expressions for the three particular moments, namely the moments of the discrete forcing terms R i , the moments of the equilibrium density population f eq i and the moments of the expanded population ζ i have to be configured. Table A1 summarizes the expressions of forcing moments.  (14)) 0 F α 0 Guo, et al. (Equation (15)) 0 F α F α u β + u α F β Kupershtokh, et al.
(Equation (16) The expressions for the moments of f eq i and ζ i were obtained accordingly as: Moments of f eq i : Moments of ζ i :