A Modified Enthalpic Lattice Boltzmann Method for Simulating Conjugate Heat Transfer Problems in Non-Homogeneous Media

Abstract

In this study, we introduced modifications to a prior existing enthalpic lattice Boltzmann method (LBM) tailored for simulating the conjugate heat transfer phenomena in non-homogeneous media with time-dependent thermal properties. Our approach is based upon the incorporation of the remaining terms of a conservative energy equation, excluding only the terms regarding flow compressibility and viscous dissipation, thereby accounting for the local and transient variations in the thermophysical properties. The solutions of verification tests, comprising assessments of both transient and steady-state solutions, validated the accuracy of the proposed model, further bolstering its reliability for analyzing heat transfer processes. The modified model was then used to perform an analysis on structured cavities under free convection, revealing compelling insights, particularly regarding transient regimes, demonstrating that the structured cavities exhibit a beneficial impact on enhancing the heat transfer processes, hence providing insights for potential design enhancements in heat exchangers. These results demonstrate the potential of our modified enthalpic LBM approach for simulating complex heat transfer phenomena in non-homogeneous media and structured geometries, offering valuable results for heat exchanger engineering and optimization.

1. Introduction

Natural convection processes have demonstrated numerous applications in refrigerating electronic devices and technology, including the cooling of printed circuit boards, memory devices, processors, and more [1,2]. This heat transfer phenomenon has also been utilized in the cooling mechanisms for buildings, electric motors, air conditioning and refrigeration systems, nuclear reactors, among others [3,4]. Typically, in these applications, the cooling process by natural convection is facilitated through heat exchangers, which come in various geometrical configurations and types. Given its widespread use and low operational cost, enhancing this heat transfer mechanism is of significant interest.

Various methods employed for this purpose involve surface modifications, such as adding fins or utilizing small cavities or grooves. Several studies found in the literature [5,6,7,8,9] investigate the impact of surface modifications on enhancing free convection heat transfer processes. Generally, these studies examine surface modifications with configurations resembling waves, employing sinusoidal functions [10,11] or evenly spaced geometries like triangular [8,12] or square [9,13] patterns. These works are often categorized based on the modifications made to horizontal [11,12], vertical [6,9,10], or inclined surfaces [7,13], typically employing experimental and Computational Fluid Dynamics (CFD) techniques to analyze the problem.

By its nature, the free-convection process occurs under conjugate heat transfer conditions, initially defined by [14]. Consequently, the numerical simulation of such a problem proves to be challenging, with many works resorting to directly imposing the boundary conditions at the fluid nodes [6,8]. Traditionally, computational procedures primarily rely on finite differences, finite volume, or finite element methods [15,16,17,18]. These methods typically involve creating coupling schemes at the solid–fluid interface, where Dirichlet and Neumann boundary conditions are imposed on different sides and the energy conservation equation is solved separately. Despite their applicability, implementing these methods on more complex interfaces can become challenging, particularly in ensuring the continuity at the interface.

On the other hand, the lattice Boltzmann method (LBM) is a mesoscopic simulation technique that solves the Boltzmann transport equation on a discrete lattice [19], offering advantages over the traditional methods. Specifically regarding the LBM models for conjugate heat transfer, the two most accurate LBM models are [20] the direct interface treatment, originally proposed by Li et al. [21], and the “enthalpy like” models, such as the one developed in Karani and Huber [22], Rihab et al. [23], Chen et al. [24].

In general, the models based on the direct interface treatment have good precision for any boundary position in relation to the boundary lattices. However, in the presence of structured surfaces or complicated geometries, the implementation of these interface conditions becomes considerably difficult. Alternatively, the “enthalpy-like” models work through the addition of specific source terms, tailored to correctly represent the target energy equation, thus not requiring special treatment for the interfaces. In this line, several works have been proposed [22,24,25,26,27].

Chen et al. [24] incorporated through a source term the effects of the local variation regarding the thermophysical properties. Additionally, the authors suggested the use of the forward Euler time discretization scheme to include the time dependence of the thermophysical properties. On the other hand, Hosseini et al. [25] proposed a method with an iterative scheme for dealing with the temperature-dependent heat capacity, which resorts to the use of a root-finding algorithm to correctly account for its variations. More recently, Kiani-Oshtorjani et al. [27], by working with a simplified conservative form of the energy equation and disregarding the flow compressibility, viscous dissipation, and time dependence of the thermophysical properties, introduced a new source term based on the variation in the thermal conductivity between the different media.

As observed in our literature review, several works employed the use of a simplified conservative energy equation [24,26,27], often disregarding the time dependence of the thermophysical properties, like the system capacitance ($\rho c_p$). To address this issue, the authors of this paper proposed the derivation of a new source term, considering its derivation from the conservative form of the energy conservation equation [28] for a system without a volumetric heat source, written in terms of temperature, T, and constant-pressure specific heat, $c_P$. Based on the work of Chen et al. [24], we derived a simple source term, which naturally accounts for the time and temperature dependence of the system’s thermophysical properties. It should be noted that the reference temperature employed for computing the body force by the Boussinesq approximation is also changed with time. In fact, this is the temperature employed for calculating the thermophysical properties. The proposal of these modifications is the main novelty of this paper.

The paper is structured into three main sections. The first section (Section 2) outlines the employed methodology, which includes the presentation of the momentum and thermal LBM models, as well as the derivation and integration of a source term to address the conjugate heat transfer. Following this, the second section (Section 3.1) presents the results obtained from the benchmark tests. These tests involve comparisons with analytical solutions for pure diffusion problems, both transient and stationary, as well as convection diffusion problems. Additionally, a mesh refinement study is conducted using natural convection results sourced from the existing literature. The final section (Section 3.2) examines the effects of the structural geometries on the natural convection problem. Finally, the principal conclusions drawn from the study are presented in the conclusion (Section 4).

2. Materials and Methods

In this section, the mathematical models employed for performing the simulations of this paper are presented. We consider the dimensional approach proposed by Martins et al. [29] given its ability to deal directly with real fluid properties.

2.1. Lattice Boltzmann Method for Fluid Flow

For the fluid flow modeling, the traditional LBE with the Bhatnagar–Gross–Krook (BGK) collision operator was used [30], which considers only a single relaxation time, $\tau$. Given time and space discrete intervals $\Delta t$ and $\Delta x$, the LBE for the fluid flow is provided by Equation (1) [31]. In this equation, $f_i$ is the density distribution function, $f_i^{eq}$ represents the equilibrium distribution function, and $S_{f_i}$ is the forcing term, which depends on the forcing scheme chosen. The sub-index i represents each discrete velocity direction, while ${\mathbf{c}}_{\mathbf{i}}$ is the particle velocity vector, according to the selected velocity scheme [32].

$$f_i(\mathbf{x}+{\mathbf{c}}_{\mathbf{i}}\Delta t,t+\Delta t)-f_i(\mathbf{x},t)=-\frac{\Delta t}{\tau }\left[f_i(\mathbf{x},t)-f_i^{eq}(\mathbf{x},t)\right]+S_{f_i}(\mathbf{x},t)\Delta t.$$

(1)

The equilibrium distribution functions for the fluid flow LBE are calculated by Equation (2) [33], where $\mathbf{u}$ and $\rho$ are the macroscopic velocity and density, respectively, $c_s$ is the sound speed, and $w_i$ indicates the weights related to each velocity direction i, both depending on the velocity scheme considered.

$$f_i^{eq}(\mathbf{x},t)=w_i\rho (\mathbf{x},t){\left[1+\frac{{\mathbf{c}}_i·\mathbf{u}}{c_s^2}+\frac{{({\mathbf{c}}_i·\mathbf{u})}^2}{2c_s^4}-\frac{\mathbf{u}·\mathbf{u}}{2c_s^2}\right]}_{(\mathbf{x},t)}.$$

(2)

In the simulations performed here, we used a two-dimensional approach with the $D2Q9$ velocity scheme [32]. Thus, according to [32], the particle velocity vectors and the respective weights are provided by Equations (3) and (4), being $c=\Delta x/\Delta t$ and $c_s=c/\sqrt{3}$.

$${\mathbf{c}}_{\mathbf{i}}=c\left\{\begin{array}{c}(0,0),i=0,\hfill \\ (1,0),(0,1),(-1,0),(0,-1),i=1,...,4,\hfill \\ (1,1),(-1,1),(-1,-1),(1,-1),i=5,...,8.\hfill \end{array}\right.$$

(3)

$$w_i=\left\{\begin{array}{c}4/9,i=0,\hfill \\ 1/9,i=1,...,4,\hfill \\ 1/36,i=5,...,8.\hfill \end{array}\right.$$

(4)

For the source term in Equation (1), the forcing scheme proposed by Guo et al. [34] was used. Then, the source term is defined by Equation (5), where $\mathbf{F}$ is the total specific force acting over the domain.

$$S_{f_i}=\left(1-\frac{\Delta t}{2\tau }\right)w_i{\left[\frac{{\mathbf{c}}_{\mathbf{i}}-\mathbf{u}}{c_s^2}+\frac{({\mathbf{c}}_{\mathbf{i}}·\mathbf{u}){\mathbf{c}}_{\mathbf{i}}}{c_s^4}\right]}_{(\mathbf{x},t)}·\mathbf{F}(\mathbf{x},t)$$

(5)

Additionally, the macroscopic fields can be obtained from the distribution functions as depicted in Equations (6) and (7), q being the number of discrete velocity directions (for example, $q=9$ for the $D2Q9$ scheme).

$$\rho (\mathbf{x},t)=\sum _{i=0}^{q-1}{f}_i(\mathbf{x},t)$$

(6)

$$\rho (\mathbf{x},t)\mathbf{u}(\mathbf{x},t)=\sum _{i=0}^{q-1}{f}_i(\mathbf{x},t){\mathbf{c}}_{\mathbf{i}}+\frac{\Delta t}{2}\mathbf{F}(\mathbf{x},t)$$

(7)

Through the Chapman–Enskog analysis [35], it is possible to recover the Navier–Stokes equation, Equation (8), from the LBE. This process results in a relationship between the relaxation time, $\tau$, and the kinematic viscosity, $\nu$, which is provided by $\nu =(\tau -0.5\Delta t)c_s^2$.

$$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla ·\left(\rho \mathbf{u}\mathbf{u}\right)=-\nabla p+\nabla ·\left[\nu \rho \left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right)\right]+\mathbf{F}.$$

(8)

Regarding boundary conditions (BC), the bounce-back (BB) rule was implemented for both static walls and fluid inlet BCs, provided by Equation (9), where $\overline{i}$ indicates the opposite direction of i, ${\mathbf{x}}_{\mathbf{b}}$ indicates the boundary nodes, and the variables with sub-index w stand for the values defined at the boundary. The more complexes geometries, such as the cavities for the natural convection, were also implemented using the BB rule.

$$f_{\overline{i}}({\mathbf{x}}_{\mathbf{b}},t+\Delta t)=f_i^{*}({\mathbf{x}}_{\mathbf{b}},t)-2w_i{\rho }_w\frac{{\mathbf{c}}_{\mathbf{i}}·{\mathbf{u}}_{\mathbf{w}}}{c_s^2}$$

(9)

The periodic kind of BC was simply implemented considering that the populations leaving one side are the same that enter at the opposite side of the domain. This relation can be expressed by the following expression: $f_i({\mathbf{x}}_{\mathbf{b}},t+\Delta t)=f_i^{*}({\mathbf{x}}_{\mathbf{b}}+\mathbf{L}-{\mathbf{c}}_{\mathbf{i}}\Delta t,t)$, where $\mathbf{L}$ is the size of the domain pointing at the normal direction of the boundary.

2.2. Boussinesq Approach for Natural Convection

In natural convection, the local temperature variations cause small local density changes, which create mass fluxes through the domain because of the gravitational field influence. These fluxes are commonly denominated “convection currents”, which help to transfer the heat in the entire domain.

A practical way to include the effects of density fluctuations in the LB simulations is by using the Boussinesq approach. Instead of considering a temperature-dependent density, $\rho \left(T\right)$, it is assumed that the fluid remains at a mean density, $\overline{\rho }$, measured at the reference temperature of $\overline{T}$. Then, the thermal expansion coefficient can be expressed according to Equation (10), and the buoyancy force that is provided by Equation (11) is added to the momentum LBE (Equations (1) and (5)) for taking into consideration the density fluctuations. In this equation, $\mathbf{g}$ is the gravitational acceleration.

$${\beta }_{exp}=-\frac{1}{\rho }\frac{\partial \rho }{\partial T}\approx -\frac{\left[1-\overline{\rho }/\rho (\mathbf{x},t)\right]}{T(\mathbf{x},t)-\overline{T}}$$

(10)

$$\mathbf{F}(\mathbf{x},t)={\mathbf{F}}_b(\mathbf{x},t)=-{\beta }_{exp}\mathbf{g}\overline{\rho }\left[T(\mathbf{x},t)-\overline{T}\right]$$

(11)

In this paper, we are simulating transient conjugate heat transfer problems, including transient natural convection problems, which can be driven by a constant heat flux or surface temperature. For all cases, the determination of the reference temperature ($\overline{T}$) becomes a challenge and a necessity. This is because, even if the fluid is initially at a uniform temperature, the local temperature values increase non-equally over time. Then, the mean temperature of the domain changes every time and it is needed to determine $\overline{T}$ for each time-step.

Consequently, in this work, all the fluid properties that are independent of the distribution functions are re-calculated for each new $\overline{T}$ value. These properties include the dynamic viscosity, specific heat capacity, conductivity, and thermal expansion coefficient. For this task, we propose to use polynomial approaches for determining the variation in these properties with the temperature. Here, these polynomials were obtained with the help of EES software [36]. Hence, the employed LBM model has the capability of dealing with local and temporal thermophysical property changes. The new $\overline{T}$ values are calculated explicitly, i.e., considering the average mean temperature of the fluid domain from the previous time step.

2.3. Lattice Boltzmann Method for Conjugate Heat-Transfer

As noted before, the main advantage of the “enthalpy-like” models is that no special treatment for the interfaces is needed due to the effects being naturally incorporated into the LBE through a source. Furthermore, for the link-wise approach [20], both methods have second-order accuracy when using the link-wise approach for the boundary nodes. Thus, in this paper, we propose a modification of the Chen et al. [24] model for accounting the temporal and local variations in the medium (solid and fluid) properties, considering conjugate heat transfer problems. The procedure for re-calculating the thermophysical properties as a function of the temperature was explained in the previous Section 2.2.

The conservative form of the energy conservation equation [28] for a system without volumetric heat source, written in terms of temperature, T, and constant-pressure specific heat, $c_P$, for a generic system is provided by Equation (12). In this equation, the term $\left(\frac{\partial lnV-}{\partial lnT}\right)\frac{DP}{dt}$ stands for the energy transfer due to compressible effects, while $\mathit{\tau }:\nabla \mathbf{u}$ refers to the viscous dissipation.

$${\partial }_t\left(\rho c_{p}T\right)+\nabla ·\left(\rho c_{p}T\mathbf{u}\right)=\nabla \left(k\nabla T\right)+\left(\frac{\partial lnV-}{\partial lnT}\right)\frac{DP}{Dt}+\mathit{\tau }:\nabla \mathbf{u}+T\frac{D\rho c_P}{Dt}$$

(12)

Now, the left-hand side of Equation (12) can be re-written as follows:

$${\partial }_t\left(\rho c_{p}T\right)+\nabla ·\left(\rho c_{p}T\mathbf{u}\right)=\rho c_P[{\partial }_{t}T+\nabla ·\left(T\mathbf{u}\right)]+T\frac{D\rho c_P}{Dt}.$$

(13)

Substituting Equation (13) into Equation (12), blue we have

$$\rho c_P[{\partial }_{t}T+\nabla ·\left(T\mathbf{u}\right)]=\nabla \left(k\nabla T\right)+\left(\frac{\partial lnV-}{\partial lnT}\right)\frac{DP}{Dt}+\mathit{\tau }:\nabla \mathbf{u}$$

(14)

Neglecting the viscous heat dissipation, the energy conservation equation is provided by Equation (15). It is important to note that this equation is generic enough to be valid for a non-homogeneous domain.

$$\rho c_P[{\partial }_{t}T+\nabla ·\left(T\mathbf{u}\right)]=\nabla \left(k\nabla T\right)$$

(15)

Similarly to the model developed by Chen et al. [24], a new variable $h_0={\left(\rho c_P\right)}_{0}T$ is defined, ${\left(\rho c_P\right)}_0$ being a constant reference heat capacitance. Substituting $T=h_0/{\left(\rho c_P\right)}_0$ into Equation (15) and defining $\sigma =\rho c_p/{\left(\rho c_P\right)}_0$, we obtain Equation (16), where $\alpha =k/\left(\rho c_p\right)$ is the local thermal diffusivity.

$${\partial }_t{h}_0+\nabla ·\left(\mathbf{u}{h}_0\right)=\nabla (\alpha \nabla h_0)-\frac{k}{{\left(\rho c_P\right)}_0}\nabla h_0\nabla \left(\frac{1}{\sigma }\right)$$

(16)

Defining a new distribution function $g_i$ related to $h_0$, and using the source term proposed by Seta [37], the LBE with the dimensional approach for the energy conservation equation is described by Equation (17). In this equation, ${\Omega }_{gi}(\mathbf{x},t)$ is the collision operator. Considering the traditional BGK approach, it is defined as ${\Omega }_i^{BGK}=-\left[g_i-g_i^{eq}\right]/{\tau }_T$. Also, $g_i^{eq}$ is the equilibrium distribution function, provided by Equation (18).

$$g_i(\mathbf{x}+{\mathbf{c}}_{\mathbf{i}}\Delta t,t+\Delta t)-g_i(\mathbf{x},t)={\Omega }_i(\mathbf{x},t)\Delta t+\Delta t{w}_i\left(1-\frac{\Delta t}{2{\tau }_T}\right)S_g(\mathbf{x},t)$$

(17)

$$g_i^{eq}=w_i{h}_0(\mathbf{x},t){\left[1+\frac{{\mathbf{c}}_i·\mathbf{u}}{c_s^2}+\frac{{({\mathbf{c}}_i·\mathbf{u})}^2}{2c_s^4}-\frac{\mathbf{u}·\mathbf{u}}{2c_s^2}\right]}_{(\mathbf{x},t)}$$

(18)

By Chapman–Enskog analysis, it is possible to recover Equation (16) from the thermal LBE (Equation (17)), resulting that the relaxation time for the BGK collision operator is related to the thermal diffusivity of the domain as $\alpha =k/\left(\rho c_p\right)$ by the following expression: $\alpha =({\tau }_T-0.5\Delta t)c_s^2$. The procedure for this analysis is provided in Appendix A. Additionally, it is possible to note that, by this expansion, the LBE only recovers the the first term on the right-hand side of Equation (16). Consequently, the other terms must be included in the LBE via the source term, which is defined by Equation (19).

$$S_g(\mathbf{x},t)=-\frac{k}{{\left(\rho c_P\right)}_0}\nabla h_0\nabla \left(\frac{1}{\sigma }\right)$$

(19)

The additional terms that must be added via source term to the LBE are not easy to calculate as it involves gradient of quantities depending on the distribution functions $g_i$. Thus, it is interesting to conduct some manipulations using the definition of heat flux from Karani and Huber [22], resulting in the following representation:

$$-\frac{k}{{\left(\rho c_p\right)}_0}\nabla h_0=-\sigma \frac{k}{\rho c_p}\nabla h_0=\sigma \left(1-\frac{\Delta t}{2{\tau }_T}\right)\sum _i\left[g_i-g_i^{eq}\right]{\mathbf{c}}_{\mathbf{i}}.$$

(20)

Also, the spatial derivative of $1/\sigma$ can be evaluated by an isotropic scheme as $\nabla \chi (\mathbf{x},t)={(c_s^2\Delta t)}^{-1}{\sum }_i{w}_i\chi (\mathbf{x}+{\mathbf{c}}_{\mathbf{i}}\Delta t,t+\Delta t){\mathbf{c}}_{\mathbf{i}}$ [19]. Then, for the BGK collision operator, the source term can be re-defined by Equation (21).

$$S_g^{BGK}(\mathbf{x},t)=\sigma \left(1-\frac{\Delta t}{2{\tau }_T}\right)\left[\sum _i\left(g_i-g_i^{eq}\right){\mathbf{c}}_{\mathbf{i}}\right]·\left[\frac{1}{c_s^2\Delta t}\sum _i{w}_i\frac{1}{\sigma (\mathbf{x}+{\mathbf{c}}_{\mathbf{i}}\Delta t,t+\Delta t)}{\mathbf{c}}_{\mathbf{i}}\right]$$

(21)

The local temperature values in each lattice can be calculated from the distribution functions as follows:

$$T(\mathbf{x},t)=\frac{1}{{\left(\rho c_p\right)}_0}\left[\sum _{i=0}^{q-1}{g}_i(\mathbf{x},t)+\frac{\Delta t}{2}{S}_g^{BGK}(\mathbf{x},t)\right]$$

(22)

There are cases where the BGK collision operator can suffer from instabilities. Then, the multiple-relaxation-time (MRT) collision operator can be used to prevent these issues. This operator is defined by ${\Omega }_i^{MRT}=-{\left[{\mathbf{M}}^{-1}\mathbf{\Lambda }\mathbf{M}\right]}_{ij}(g_j-g_j^{eq})$, where $\left[\mathbf{M}\right]$ is the transformation matrix, responsible for transforming the distribution functions to the moment space and vice-versa. $[\mathbf{\Lambda }]$ is the collision matrix, which is usually a diagonal matrix containing the relaxation rates related to each moment of the distribution function.

Considering the $D2Q9$ velocity scheme, the transformation matrix $\left[\mathbf{M}\right]$ for the dimensional LBM is defined in Martins et al. [29], and the collision matrix here is defined as $[\mathbf{\Lambda }]=diag({\omega }_0,{\omega }_1,...,{\omega }_8)$. The relaxation frequencies can be provided by ${\omega }_3^{-1}={\omega }_5^{-1}=\alpha /c_s^2+\Delta t/2$ [38], while the other frequencies can be arbitrarily chosen in such a way that guarantees the stability of the method, usually assuming values close to $\Delta t^{-1}$.

Additionally, it is very common to perform the collision process with the MRT collision operator directly in the moment space for facility. Thus Equation (17) is transformed to the moment space multiplying the equation by $\left[\mathbf{M}\right]$. Defining the moments of the distribution functions by $\mathbf{m}=\left[\mathbf{M}\right]\mathbf{g}$, the collision process in the moment space can be represented by Equation (23). Then, the post-collision distribution functions can be recovered applying the inverse of the transformation matrix: ${\mathbf{g}}^{*}={\left[\mathbf{M}\right]}^{-1}{\mathbf{m}}^{*}$.

$${\mathbf{m}}^{*}(\mathbf{x},t)-\mathbf{m}(\mathbf{x},t)=-\Delta t[\mathbf{\Lambda }]\left[\mathbf{m}(\mathbf{x},t)-{\mathbf{m}}^{\mathbf{eq}}(\mathbf{x},t)\right]+\Delta t\mathbf{w}\left\{\left[\mathbf{I}\right]-\frac{\Delta t}{2}[{\mathbf{\Lambda }}_{]}\right\}{S}_g^{MRT}(\mathbf{x},t)$$

(23)

It is important to note that the spatial derivative of $h_0$ was previously defined for the BGK collision operator. Then, it needs to be adapted for the MRT collision operator, resulting in Equation (24).

$$-\frac{k}{{\left(\rho c_p\right)}_0}\nabla h_0=-\sigma \frac{k}{\rho c_p}\nabla h_0=\sigma \sum _i{\left[{\left[\mathbf{M}\right]}^{-1}\left(\left[\mathbf{I}\right]-\frac{\Delta t}{2}[\mathbf{\Lambda }]\right)\left(\mathbf{g}-{\mathbf{g}}^{\mathbf{eq}}\right)\right]}_i{\mathbf{c}}_{\mathbf{i}},$$

(24)

Thus, we have that the source term for the MRT collision operator will be provided by Equation (25).

$$S_g^{MRT}(\mathbf{x},t)=\sigma \sum _i{\left[{\left[\mathbf{M}\right]}^{-1}\left(\left[\mathbf{I}\right]-\frac{\Delta t}{2}[\mathbf{\Lambda }]\right)\left(\mathbf{g}-{\mathbf{g}}^{\mathbf{eq}}\right)\right]}_i{\mathbf{c}}_{\mathbf{i}}·\left[\frac{1}{c_s^2\Delta t}\sum _i{w}_i\frac{1}{\sigma (\mathbf{x}+{\mathbf{c}}_{\mathbf{i}}\Delta t,t+\Delta t)}{\mathbf{c}}_{\mathbf{i}}\right]$$

(25)

Here, we proposed a model based on the Chen et al. [24] enthalpy-based LBM. Our method was deduced from the conservative energy conservation equation. Through this formulation, no hypotheses regarding homogeneous media, or time and local independence of thermal properties, were composed, and, as such, both the local and temporal variations in the thermodynamics properties are allowed. The spatial variations enter through the source term in the thermal LBE (Equation (17)); meanwhile, the temporal variations are naturally accounted through our formulation, granting the correct capture of the transient phenomena.

For the thermal boundary conditions, the BB rule for adiabatic walls was considered as provided by Equation (26). Additionally, for the Dirichlet type of boundary condition, the anti-BB rule was considered; see Equation (27). In this relation, $T_w$ corresponds to the imposed temperature at the boundary. For the Neumann type of boundary condition, we considered the scheme proposed in Martins et al. [39]. Please consult the cited references for all the needed details.

$$g_{\overline{i}}({\mathbf{x}}_{\mathbf{b}},t+\Delta t)=g_i^{*}({\mathbf{x}}_{\mathbf{b}},t)$$

(26)

$$g_{\overline{i}}({\mathbf{x}}_{\mathbf{b}},t+\Delta t)=2w_i{T}_w\left[1+\frac{{\left({\mathbf{c}}_{\mathbf{i}}·{\mathbf{u}}_{\mathbf{w}}\right)}^2}{2c_s^4}-\frac{{\mathbf{u}}_w·{\mathbf{u}}_w}{c_s^2}\right]-g_i^{*}({\mathbf{x}}_{\mathbf{b}},t)$$

(27)

The periodic kind of BC was simply implemented considering that the populations leaving one side are the same that enter at the opposite side of the domain, as occurred for the fluid flow: $g_i({\mathbf{x}}_{\mathbf{b}},t+\Delta t)=g_i^{*}({\mathbf{x}}_{\mathbf{b}}+\mathbf{L}-{\mathbf{c}}_{\mathbf{i}}\Delta t,t)$.

3. Results and Discussion

3.1. Benchmark Tests

Before exploring the natural convection in an enclosure of a heat exchanger structured with cavities and under conjugate heat transfer conditions, three benchmark tests are explored for validating the proposed methodology. These include a 1D conduction between three different solids, a forced convection between two fluids, and the natural convection in a partially heated enclosure considering the conjugate heat transfer conditions. The numerical results were compared against reference solutions, analytical for the first two and numerical for the third, for validation purposes. These comparisons were quantitatively evaluated using the global relative error, defined by Equation (28).

$$E_R=100\%\left[\sqrt{\frac{\sum {(T_{analytical}-T_{LBM})}^2}{\sum T_{analytical}^2}}\right]$$

(28)

3.1.1. Heat Diffusion between Three Solids

The first benchmark problem solved is the transient one-dimensional conduction between three different solids. A square domain with side $L=0.60$ m is considered, the first layer being composed of copper (starting from the left of the domain), the second silicon, and the third aluminum. Each layer has a length of $d=0.20$ m, and then the interfaces are located at $x_1=0.20$ m and $x_2=0.40$ m. Both the left and right walls are set with constant temperatures (Dirichlet’s boundary conditions) about $T_{wal{l}_l}=313.15K$ and $T_{wal{l}_r}=293.15$ K, respectively.

Here, all the solid’s properties are assumed to be constants, calculated at a temperature of $293.15K$. k is the thermal conductivity; the properties are for copper $k_{Cu}=401.2$ ${\mathrm{Wm}}^{-1}{K}^{-1}$, $c_{p_{Cu}}=384.5$ ${\mathrm{Jkg}}^{-1}{K}^{-1}$, and ${\rho }_{Cu}=8934.0$ ${\mathrm{kgm}}^{-3}$; for silicon $k_{Si}=149.6$ ${\mathrm{Wm}}^{-1}{K}^{-1}$, $c_{p_{Si}}=709.1$ ${\mathrm{Jkg}}^{-1}{K}^{-1}$, and ${\rho }_{Si}=2330.0$ ${\mathrm{kgm}}^{-3}$; and, for aluminum, $k_{Al}=237.0$ ${\mathrm{Wm}}^{-1}{K}^{-1}$, $c_{p_{Al}}=905.0$ ${\mathrm{Jkg}}^{-1}{K}^{-1}$, and ${\rho }_{Al}=2707.0$ ${\mathrm{kgm}}^{-3}$. The simulation was performed using the $D1Q3$ velocity scheme. The space and time were discretized considering intervals of $\Delta x=1.2·{10}^{-3}$ m and $\Delta t=2.5·{10}^{-3}$ s, respectively. The analytical solution for this problem is depicted in Appendix B [40].

The results for several time steps are depicted in Figure 1; as expected, the major divergences between both the analytical and the numerical solutions occur at the interfaces. Despite that, the methods are shown to capture the overall dynamics of the problem over the transient regime and also correctly predict the stationary regime.

Additionally, the global errors between the numerical and the analytical solution, estimated through Equation (28), are provided in

Table 1. Total errors for the convection–diffusion problem at different times.

t $\left[\mathbf{s}\right]$	Global Errors $·{10}^2[\%]$	t $\left[\mathbf{s}\right]$	Global Errors $·{10}^2[\%]$
$100.0$	$0.0152$	$900.0$	$0.0435$
$200.0$	$0.0337$	$1000.0$	$0.0424$
$300.0$	$0.0434$	$1100.0$	$0.0414$
$400.0$	$0.0472$	$1200.0$	$0.0407$
$500.0$	$0.0481$	$1300.0$	$0.0401$
$600.0$	$0.0474$	$1400.0$	$0.0396$
$700.0$	$0.0462$	$1500.0$	$0.0393$
$800.0$	$0.0448$	$1600.0$	$0.0391$

. From a quantitative point of view, both solutions present good agreement, resulting in relatively small errors (below $0.05·{10}^{-2}\%$) for all the time steps, thus further corroborating the results shown in Figure 1.

Analyzing both results (Figure 1 and Table 1), it can be affirmed that the employed model is shown to capture both the transient and stationary regimes referring to the heat diffusion between different media with relatively good accuracy, returning low global errors between the analytical and numerical solutions.

3.1.2. Convection–Diffusion with a Flat Interface

Next, a transient convection–diffusion problem with a thin flat interface is studied. The solution for this problem is provided in Li et al. [21]. This problem consists of a channel with fixed height H and length L, separated in two subdomains: the bottom half, with water flowing, and the top half, with carbon dioxide flowing. A periodic boundary condition is imposed on the left and right boundaries, while a Dirichlet boundary condition in the form of Equation (29) is imposed on both the top and bottom boundaries. Both fluids are assumed to have a constant uniform velocity $u=\frac{Pe{\alpha }_w}{H}$ in the x-direction, and $\Delta T$ is set to $20K$. The analytical solution for this problem is provided in Appendix C.

$$T\left(x,0,t\right)=T\left(x,L,t\right)=T_0+\Delta Tcos\left(\frac{2\pi }{L}x+\frac{2\pi St{k}^2{\alpha }_w}{H^2}t\right)=T_0+\Delta Tcos\left(Kx+\omega t\right)$$

(29)

where $Pe$ refers to the Péclet number and $Stk$ to the Stokes number, and they were set as 20 and 1, respectively. The properties were assumed as ${\rho }_w=978.16114$ kg ${\mathrm{m}}^{-3}$, $c_{P_w}=4188.10655$J ${\mathrm{Kg}}^{-1}$${\mathrm{K}}^{-1}$, and ${\alpha }_w=1.61·{10}^{-7}$${\mathrm{m}}^2{s}^{-1}$ for the water, evaluated at $P_w=\mathrm{101,325}Pa$ and $T_{ref}=343.15K$; and ${\rho }_{C{O}_2}=177.45$kg ${\mathrm{m}}^{-3}$, $c_{P_{C{O}_2}}=1688.66$J ${\mathrm{Kg}}^{-1}{K}^{-1}$, and ${\alpha }_{C{O}_2}=9.76·{10}^{-8}$${\mathrm{m}}^2$${\mathrm{s}}^{-1}$ for $C{O}_2$, evaluated at $P_{C{O}_2}=\mathrm{8,115,030}Pa$ and $T_{ref}=343.15$ K. The simulation is performed using a two-dimensional velocity scheme, the $D2Q9$. For numerical simulation, the space and time were discretized with $\Delta x=1·{10}^{-4}$ m and $\Delta t=1·{10}^{-4}$ s.

To facilitate the comparisons, the results are shown for various time fractions related to one period, $\Gamma =1/\omega$, considering that the solution is periodic in time. Figure 2 shows both the two-dimensional temperature distribution obtained with the LBM (Figure 2a) and the numerical and analytical temperature profiles along different cuts in the x-axis (Figure 2b) for a fixed time fraction $t=0.75\Gamma$. Although Figure 2 represents only a single time, we can still drawn some important conclusions about it. First, from Figure 2b, we can see that both the temperature profiles present good agreement. The numerical results (hollow circles) almost entirely coincide with the analytical solution (solid lines). Next, we can also observe that, yet again, the major divergences between the solutions will occur at the interface, as expected.

Additionally, to better study the performance of our proposed model, we also estimated the global errors (Equation (28)) between the numerical and analytical solutions. The results are provided in

Table 2. Global errors for the convection–diffusion problem at different times.

$\mathit{t}/\mathbf{\Gamma }$	Global Errors $·{10}^2[\%]$	$\mathit{t}/\mathbf{\Gamma }$	Global Errors $·{10}^2[\%]$
$0.125$	$2.1$	$0.625$	$1.65$
$0.250$	$1.58$	$0.750$	$1.65$
$0.375$	$1.65$	$0.875$	$1.65$
$0.500$	$1.66$	$1.000$	$1.65$

. Overall, the results are shown to be in good agreement, with error values below $3·{10}^{-2}\%$. It is interesting to observe that, for the initial moments, the output errors were larger, about $2·{10}^{-2}\%$. These results can be associated with difficulty in precisely matching the initial conditions for both solutions. Finally, we can also see that, once the initial times have passed, the errors will both diminish and stabilize, resulting in an almost constant error for the remaining studied times.

3.1.3. Natural Convection with a Fixed Heat Flux

In this case, a natural convection problem consisting of a closed square cavity with a partially heated bottom was tested. The problem configuration, along with its main dimensions, are presented in Figure 3; the problem was presented and simulated in Cheikh et al. [41]. It is important to note that some modification were made to the original problem: instead of imposing heat directly on the fluid, the heat flux was applied to a thin layer of copper, with a thickness equal to one $dx$, thus also serving to evaluate the conjugate model.

In the simulations, the Rayleigh number was maintained at $Ra={10}^4$ and the dimensionless length of the heat source was set as $ϵ=0.8$. Moreover, serving the purpose of validating the employed model, this test was performed with the intent of choosing a reasonable mesh size to correctly predict the natural convection phenomenon that will be studied in the proceeding chapters of this paper. For the boundary conditions, the adiabatic wall (black lines in Figure 3) and the no-slip boundary conditions were implemented via a bounce-back scheme, while the Dirichlet boundary conditions were implemented via an anti-bounce-back scheme. Lastly, for the Neumann boundary conditions, the scheme proposed in Martins et al. [39] was used.

First, the results were organized with reference to the mesh size in

Table 3. Nusselt number errors for a natural convection problem.

Reference Nusselt Number: $\mathit{Nu}=3.80$
$\Delta \mathbf{x}·{\mathbf{10}}^{\mathbf{5}}$[m]	$\overline{\mathit{Nu}}$	Errors $\overline{\mathit{Nu}}$ [%]
$2.500$	$3.72$	$2.17$
$1.250$	$3.75$	$1.39$
$0.625$	$3.80$	$0.07$

to facilitate the comparison between different meshes. The Nusselt number was calculated using the procedure presented in Appendix D. As we expect, a clear reduction in the relative errors can be observed with the mesh refinement. The model is shown to correctly predict the Nusselt number for all the meshes with relative errors below $3\%$, even with the larger meshes, thus advocating for its accuracy. Additionally, Figure 4 presents both the temperature distributions and the streamlines for the finest mesh tested here. We can observe that, as expected, the fluid presents higher temperatures above the heating section, where the heat flux is applied. This temperature difference leads to two recirculating vortexes, with the fluid rising in the central region of the enclosure.

From the results of Table 3, we can assert that the mesh with $\Delta x=6.25·{10}^{-6}$ m produced better results, consequently presenting a better estimation of the resulting Nusselt number, with errors less than $0.1\%$. Thus, it is the most reasonable choice for the simulations presented in the following Section 3.2.

3.2. Results for Natural Convection with Structured Cavities

The study examines the natural convection under conjugate heat transfer conditions on structured cavities, analyzing the transient and steady-state operations. Both Dirichlet and Neumann boundary conditions at the enclosure’s bottom are considered to assess how cavities affect heat transfer, determining their impact on the intensifying heat transfer rates at the solid surface. We conducted natural convection simulations using the structured cavities outlined in Moreira et al. [42]. Figure 5 provides a schematic representation of the tested geometries, Geometry 1 with $\theta ={30}^{°}$, Geometry 2 with $\theta ={45}^{°}$, and Geometry 3 with $\theta ={60}^{°}$. Through these simulations, we aim to accomplish two primary objectives: firstly, to exhibit the method’s proficiency in handling complex problems, and, secondly, to investigate the influence of structured cavities on the natural convection phenomenon inside a heat exchanger under conjugate heat transfer conditions.

The structured cavities examined in this study have demonstrated the ability to enhance the flow boiling phenomenon, as shown in Moreira et al. [42]. Consequently, we aim here to investigate whether natural convection can be enhanced by the same surface structures. In fact, free convection is the first step in the nucleating boiling process. Our simulations employ water at $P=1atm$ and a solid copper base characterized by the following thermal properties: $k_{Cu}=401.2$ ${\mathrm{Wm}}^{-1}$${\mathrm{K}}^{-1}$, $c_{p_{Cu}}=384.5$ ${\mathrm{Jkg}}^{-1}$${\mathrm{K}}^{-1}$, and ${\rho }_{Cu}=8934.0$ ${\mathrm{kgm}}^{-3}$.

The simulations are categorized as follows: firstly, simulations with a fixed heat flux at the base are conducted to explore the transient regime. Subsequently, simulations with a fixed base temperature are performed to investigate both transient and steady-state regimes. In all cases the Nusselt number was calculated using the procedure presented in Appendix D.

For the hydrodynamic LBE, a non-slip boundary condition is adopted for all the boundaries, implemented via the bounce-back scheme. Regarding the energy LBE, an adiabatic condition is assumed for both lateral boundaries, again via the bounce-back scheme. For the bottom and top boundaries, Dirichlet boundary conditions are implemented via the anti-bounce-back scheme, while, for the Neumann boundary condition, the boundary condition proposed in Martins et al. [39] is utilized. Finally, considering the results presented in Section 3.1.3, the simulations in this section are performed with $\Delta x=6.25\times {10}^{-6}$m and $\Delta t=6.50\times {10}^{-6}$s.

3.2.1. Geometry Impact on Natural Convection with Imposed Heat Flux

In this section, the results for the simulations where a heat flux, $\dot{Q}wall$, was imposed at the base will be presented for each geometry configuration and for a flat surface. Two different heat fluxes were tested: $\dot{Q}wall=1\times {10}^5$ ${\mathrm{Wm}}^{-2}$ and ${\dot{Q}}_{wall}=1\times {10}^6$ ${\mathrm{Wm}}^{-2}$, resulting in Rayleigh numbers of $Ra=3.8\times {10}^4$ and $Ra=3.8\times {10}^5$, respectively. Unlike the previous works [5,6,7,8,9], we specifically address the transient regime, which can be most relevant for heat exchangers operating under real cooling conditions. For example, an electronic device can be refrigerated by a tested surface in a submerged application.

For a more comprehensive understanding of the simulation results, it is interesting to analyze the temporal evolution of the variables of interest. Therefore, Figure 6 and Figure 7 present spatially averaged values for the Nusselt number (${\overline{Nu}}_{L_{eq}}$), for the heat flux transferred across the base interface (${\overline{\dot{Q}}}_{base}$), and for the temperature difference ($T_{avg,base}-T_c$). It is important to note that all the averaged values presented were estimated considering the actual cavity length. As time evolves, both the heat flux and the temperature differences increase. However, in contrast, the predicted Nusselt number experiences a significant decline from the values observed at the outset. This effect can be associated with the tendency of heat transfer to approach a steady-state condition, reaching heat transfer rate values close to the imposed heat flux, while the temperature difference demonstrates an almost linear growth pattern.

From Figure 6 and Figure 7, it is evident that the structured cavities have a positive impact on the heat transfer process, leading to an increase in the Nusselt number compared to the plane surface. Among the studied geometries, Geometry 3 (with $\theta ={60}^{°}$) exhibits a more significant influence on the process. This enhanced heat transfer efficiency is further illustrated by the base temperature differences, as shown in Figure 6c and Figure 7c. Initially, the temperature differences for the geometries are lower than those observed for the plane surface, suggesting the more efficient heat transfer from the structured surfaces (resulting in higher heat fluxes). However, as the estimated heat transfer rate for the plane interface surpasses that of the structured geometries, these geometries begin to exhibit higher temperature differences than the plane surface.

The enhancement in the heat transfer process observed in Figure 6 and Figure 7 can be further investigated by analyzing both the temperature profiles and fluid motion. To perform such an analysis, we plotted both the two-dimensional temperature profiles and the streamlines for each geometry at a given instant. Figure 8 shows the plotted results for $Q_{wall}=1\times {10}^5$ at $t=0.25$ s, which is associated with a period of time in which the structured surface is shown to enhance the heat transfer process.

Starting with the temperature distributions, we observed a concentration of isothermal lines at the interface between the fluid and the solid wall, indicating higher heat transfer in those regions, particularly at the upper parts of the structured cavities (middle of the cavity). Additionally, from left to right, the isotherms appear more closely packed, indicating enhanced heat transfer. These observations confirm the performance shown in Figure 6 and Figure 7.

Now, considering the resultant fluid flow as shown in Figure 8b, we observe the fluid flow in and out of the cavities created by the geometries, thus avoiding stagnation. Upon closer observation, we note that the fluid still exhibits relatively low speed in relation to the surfaces, especially for the sharper geometry ($\theta ={30}^{°}$). The higher velocities and the circulating fluid can be associated with the observed results in two aspects. Firstly, this observed circulation avoids the creation of stagnation zones, where the fluid would almost be stationary, effectively reducing the thermal resistance zones. Secondly, the moving fluid advects energy from the surface, thus increasing the heat transferred between the solid and the fluid.

In addition to the streamlines and temperature plots, we also estimated the vorticity at different times for the $Q_{wall}=1\times {10}^5$ condition. The selected times shown in Figure 9 correspond to regions where the geometries exhibit better performance ($t=0.25$ s), a region where the geometries and the surfaces have similar performance ($t=1.00$ s), and, finally, regions where the plane surface surpasses the cavities ($t=2.50$ s). By comparing the results in Figure 9a with Figure 6 and Figure 8, it is evident that the higher performance of the cavities is associated with the observed vorticity, where higher levels of vorticity indicate better performance. Still considering Figure 9a, it can be seen that the smoother geometries (greater $\theta$) favor the formation of vortices near the surfaces, thereby aiding energy advection.

Considering the different times plotted in Figure 9, it is observable that all the configurations favor the formation of counter-rotating symmetrical vortices. Overall, the vorticity is shown to increase, indicating a higher circulation of the fluid within the domain and consequently the overall fluid rotation speed. Given the present configurations, the smoother geometries are shown to maintain the inner vortices for longer times. By comparing the results from Figure 9 with those of Figure 6, we can see that the presence of these vortices is fundamental for the better performance of the cavities. These vortices ensure that the “extra” length provided by the surface modifications remains in contact with the moving fluid, thus avoiding the creation of stagnation zones and favoring the heat transfer process.

Alternatively, the impact of structured cavities can be assessed through their effectiveness or enhancement factor, ${\epsilon }_f$, defined as the ratio between the heat transferred through the structured geometry and the plane geometry, respectively [43,44]. Figure 10 shows the evolution of the effectiveness of the geometries over time, confirming the previous results. Geometry 3 with $\theta ={60}^{°}$ presents the best performance, especially in the initial moments. One interesting observation concerns the evolution of effectiveness. Since the maximum heat flux is limited, it is expected that, as time advances, the effectiveness reaches a value of about $100\%$ due to energy conservation. In a steady-state condition, the heat flux at the solid–fluid interface should be the same as that imposed at the base. Thus, the steady-state effectiveness may not be the best parameter to evaluate the structured cavities.

Nonetheless, a transient analysis of ${\epsilon }_f$ can provide interesting results. By analyzing the results displayed in Figure 10a, it is observable that the effectiveness of Geometry 3 takes longer to decay, followed by Geometry 2 and then by Geometry 1. This indicates that the heat fluxes at the solid–fluid interface of structured cavities grow faster than those for the plane interface.

Additionally, a time constant for each geometry can be estimated, providing a way to evaluate the transient performance of each geometry. Since we are measuring the time needed for ${\epsilon }_f$ to reach a value of $100\%$, a higher time constant means that the heat flux for the cavities reaches the value of the stationary heat flux imposed at the base faster than the plane surface. Estimated time constants are presented in

Table 4. Time constant estimation.

${\mathit{Q}}_{\mathit{wall}}\left[{\mathbf{W}/\mathbf{m}}^2\right]$	$1·{10}^5$			$1·{10}^6$
θ °	30	45	60	30	45	60
$tau\left[s\right]$	$0.22$	$0.34$	$0.48$	$0.06$	$0.10$	$0.18$

for two heat fluxes and the three cavities As expected, Geometry 3 presented a higher estimated time constant, which again serves as an indication of its better performance.

The higher effectiveness exhibited by the geometries in Figure 10, for both heat flux values, can also be observed in the evolution of the fluid temperatures, or, alternatively, of the temperature differences shown in Figure 6c and Figure 7c. Since the structured cavities demonstrate effectiveness over $100\%$, it is expected that these geometries transfer heat more efficiently than the plane surface, resulting in higher temperature differences for the more “efficient” geometries. This observation further confirms the obtained results.

3.2.2. Natural Convection with Fixed Base Temperature

In addition to the previous results, we also conducted simulations of the natural convection process for five different base temperatures, $T_{wall}$, effectively resulting in $Ra$ values ranging between $1.5\times {10}^3$ and $2.6\times {10}^4$, considering each geometry configuration and a flat surface. The results primarily concern the average Nusselt number, $\overline{Nu}$, and the average heat flux at the solid–fluid interface, $\overline{\dot{Q}}$, throughout the transient regime.

Figure 11 illustrates the evolution of the average heat transfer rate at the solid–fluid interface. The behavior observed differs drastically betweeen each condition. For $T_{wall}\le 40$ °C, all the structured cavities are shown to have a negative impact on the heat transfer process. Conversely, when considering higher wall temperatures ($T_{wall}\ge 60$ °C), the enclosures begin to exhibit regions (typically at the beginning of the simulations) where an enhancement of the free-convection process through the addition of structured surfaces is noticeable. In other words, these are regions where an intensification of the heat transfer process occurs. This simulated behavior should be a direct consequence of the $Ra$ values. For small $Ra$ values, the cavities increase the solid–fluid thermal resistance due to the low flow vorticity, leading to low heat transfer rates. This behavior will be discussed further in this section.

Figure 12 and Figure 13 depict both the temperature distributions and the streamlines for $T_{wall}=80$ °C at $t=2.50$ s (transient) and the steady-state regimes, respectively. Similar to the tests considering an imposed heat flux, the conditions of the problem favor the formation of counter-rotating symmetrical vortices. However, at the evaluated time steps, the vortices exhibit an opposite rotation compared to what was observed in the previous section, with the fluid descending at the middle section.

Despite their similarities, Figure 12 and Figure 13 show some major differences, which may help to explain the behaviors observed during the transient regime. Looking solely at the temperature distributions, a clear difference is observed in the distribution of the isothermal lines. Figure 12 presents more closely packed lines, especially near the corners of the solid–fluid interface, indicating the presence of sharper temperature gradients. Physically, this behavior translates to higher heat transfer rates. From a hydrodynamic perspective, both conditions exhibit similar behaviors, with the exception of the observed absolute fluid velocity ($\left|\mathbf{u}\right|$), which is higher when evaluated at $t=2.50$ s compared to the stationary regime. This indicates the capability to advect more energy during the transient regime.

When evaluating the impact of structured cavities on the steady-state geometry performance, we observed a negative impact compared to a plane surface, as shown in

Table 5. Geometry comparison.

	${\mathit{T}}_{\mathit{wall}}={30}^{°}$C ${\mathit{Ra}}_{\mathit{L}}=1.5·{10}^3$				${\mathit{T}}_{\mathit{wall}}={40}^{°}$C ${\mathit{Ra}}_{\mathit{L}}=4.1·{10}^3$				${\mathit{T}}_{\mathit{wall}}={60}^{°}$C ${\mathit{Ra}}_{\mathit{L}}=1.2·{10}^4$				${\mathit{T}}_{\mathit{wall}}={80}^{°}$C ${\mathit{Ra}}_{\mathit{L}}=2.4·{10}^4$				${\mathit{T}}_{\mathit{wall}}={90}^{°}$C ${\mathit{Ra}}_{\mathit{L}}=2.6·{10}^4$
θ °	0	30	45	60	0	30	45	60	0	30	45	60	0	30	45	60	0	30	45	60
${\overline{Nu}}_{L_{eq}}$	1.03	0.94	0.84	0.8	1.03	0.94	0.84	0.8	1.86	1.74	1.55	1.43	2.75	2.56	2.28	2.1	3.04	2.84	2.52	2.32
$\overline{Q}\left[{\mathrm{kWm}}^{-2}\right]$	3.14	3.01	2.89	2.63	6.36	6.1	5.86	5.33	23.29	23.0	21.97	19.41	52.46	51.86	49.4	43.75	68.25	67.46	64.25	56.83

, which leads to a decrease in the Nusselt number, in accordance with previous studies such as [5,10]. These results can be associated with the stagnation of the fluid within the cavities, as can be seen in Figure 13, where the fluid reaches temperatures as close to the solid as possible, creating regions with high thermal resistance.

To better understand the negative impact of the structured cavities at a stationary state, we also plotted the local heat transfer rates at the solid–fluid interface; see Figure 14. Due to the symmetry of the problem, the results are only shown for half of the domain (from the middle section to the wall). Two distinct regions can be identified: outside and within the cavities. In the first region, the solid is in contact with fluid that exhibits relative motion to it. This motion increases the heat transferred from the solid to the fluid. In the second region, the fluid tends to stagnate, and its temperature approaches that of the solid regions. As a result, the heat transferred in this region is reduced. This behavior is observed for both regimes (transient and steady), but it is more pronounced for the steady state, resulting in a negative effect of the structured cavities. In the transient regime, the cavities produce a positive effect, enhancing the heat transfer rate from the heated surface to the fluid.

4. Conclusions

In this paper, we introduced modifications to an existing thermal LBM model for simulating conjugate heat transfer processes [24], aimed at accounting for the local and temporal variations in the thermophysical properties. By working with the conservative form of the energy conservation equation written on the basis of the temperature and constant-pressure specific heat, we derived a source term that encompasses both the time and spatial variations in the domain properties. We also proposed the calculation of the reference temperature for each time step. These modifications are the main novelty of the paper and enabled the numerical simulation of the natural convection process in an enclosure with structured cavities as a function of time. The following are the main conclusions of the paper.

1. The proposed LBM model was validated through the solution of three benchmark problems, including purely diffusion, advection–diffusion, and free-convection problems. Overall, the proposed modified model presented consistently low relative errors compared to the reference analytical and numerical solutions.

2. Furthermore, the proposed LBM model was employed to investigate the natural convection phenomenon within an enclosure with structured cavities. Two principal simulations were undertaken: first, cases involving the imposition of a uniform heat flux at the domain base (solid medium); second, in a more traditional way, the imposition of a uniform base temperature was also studied.

3. Alternatively to previous works, we focused our analysis on the transient regime, which revealed interesting results. In both simulations, an enhancement of the free-convection process was perceived for the structured surfaces, mainly for higher $Ra$. As discussed before, these results are associated with higher vorticities, which favors the fluid circulation within the domain, thus preventing its stagnation inside the cavities. This kind of result is very scarce in the open literature.

4. Additionally, we also examined the stationary regimes (only for the second simulation case). In these circumstances, the base problem (an enclosure with a plane interface) was shown to generate higher heat transfer rates, and consequently higher Nusselt numbers, when compared to the structured surfaces, supporting the findings from prior research [5,10]. Contrary to the transient cases, in these circumstances, some portion of the fluid stagnates inside the cavities, essentially creating additional thermal resistances, ultimately reducing the heat transfer rates.

In summary, the natural convection simulation within an enclosure with a structured surface provided interesting insights regarding this heat transfer process. First, in accordance with the previous works [5,10], the presence of a structured surface was shown to have a negative impact on the stationary heat transfer process, favoring the stagnation of the fluid between the cavities. On the other hand, when considering transient regimes, where most of the practical application occurs, the structured surfaces were shown to present a positive impact, in particular for what concerns higher $Ra$, since in this condition no fluid stagnation within the cavities was perceived.