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

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 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], evenly spaced geometries like triangular [8,12] or squared [9,13] patterns. These works are often categorized based on 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, numerical simulations of such a problem prove to be challenging, with many works resorting to directly imposing 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 solving the energy conservation equation separately. Despite their applicability, implementing these methods on more complex interfaces can become challenging, particularly in ensuring continuity at the interface.

On the other hand, the Lattice Boltzmann Method (LBM), a mesoscopic simulation technique that solves the Boltzmann transport equation on a discrete lattice [19], offers advantages over traditional methods. It excels in handling complex geometries, simulating multiphase flows [20,21,22,23,24], conjugate heat transfer [25,26,27,28,29], and capturing non-equilibrium effects [30]. Due to its inherent nature, the LBM does not necessitate specific coupling schemes at solid-liquid interfaces for simulating conjugate heat transfer problems. The method is capable of dealing with interfaces more seamlessly, which, in many aspects, can facilitate its implementation for more complex geometries.

Most of the studies identified in our literature review have focused on analyzing the impact of adding protuberances (such as fins) to enhance the natural convection process. Typically, these analyses are conducted under stationary regimes, overlooking the dynamics leading to such conditions. Furthermore, many studies do not consider the surface, imposing thermal boundary conditions directly at the solid-fluid interface. This approach may, once again, influence the transient regime of the problem. By its very nature, this type of problem involves heat transfer between different mediums, or, in other words, occurs under conjugate heat-transfer conditions. Therefore, to accurately capture the transient regime, both the solid and fluid must be modeled.

To address this issue, the authors of this paper propose the implementation of a modified model of the enthalpic thermal LBM to handle conjugate heat transfer. The LBM possesses intrinsic characteristics that make it well-suited for addressing such problems, as demonstrated throughout this study. Therefore, this work aims to develop an LBM model to address conjugate heat transfer conditions and simulate a natural convection problem with different surface geometries.

The considered computational domain characterizes a typical heat exchanger configuration employed in electronic cooling applications. The development of the conjugate heat transfer model is based upon the expansion of the complete energy conservation equation assuming an incompressible flow and neglecting viscous dissipation, and finally by inserting the resulting additional term to the lattice Boltzmann equation (LBE) through a thermal forcing scheme. Since the proposed model is based upon the inclusion of missing term of the energy equation into the LBE through a source term, and no special treatment of the interface is applied, the proposed model can be easily implemented when considering complex interfaces and extended to higher dimensions, while incorporating the physics associated with non-homogeneous mediums.

The paper is structured into three main sections. The first section (Sec. Section 2) outlines the employed methodology, which includes the presentation of momentum and thermal LBM models, as well as the derivation and integration of a source term to address conjugate heat transfer. Following this, the second section (Sec. Section 3) presents the results obtained from 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 existing literature. The final section (Sec. Section 4) examines the effects of structural geometries on the natural convection problem. Finally, the principal conclusions drawn from the study are presented in the concluding Section 5).

2. Mathematical Modeling

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

2.1. Lattice Boltzmann Method for Fluid Flow

For the fluid flow modeling it was considered the traditional LBE with the BGK collision operator [32], which considers only a single relaxation time, $\tau$. Giving a time and space discrete intervals $\Delta t$ and $\Delta x$, the LBE for the fluid flow is given by Eq. 1[30]. 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 depend 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 [33].

$$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 Eq. 2 [34], where $\mathbf{u}$ and $\rho$ are the macroscopic velocity and density, respectively, $c_s$ is the sound speed and $w_i$ are the weighs 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. Thus, according to [33] the particle velocity vectors and the respective weighs are given by Eq. 3 and Eq. 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 Eq. 1, the forcing scheme proposed by Guo et al. [35] was used. Then, the source term is defined by Eq. 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 Eq. 6 and Eq. 7, being q 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 [36] it is possible to recover the Navier-Stokes equation, Eq. 8, from the LBE. This process results in a relation between the relaxation time, $\tau$, and the kinematic viscosity, $\nu$, which is given 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), it was implemented the bounce-back (BB) rule for both static walls and fluid inlet BCs, given by Eq. 9, where $\overline{i}$ indicates the opposite direction of i, ${\mathbf{x}}_{\mathbf{b}}$ are the boundary nodes and the variables with sub-index w stand for the values defined at the boundary. The more complexes geometries, 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 the 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 Eq. 10, and the buoyancy force that is given by Eq. 11 is added to the momentum LBE (Eq. 1) 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}}_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 an 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 of these properties with the temperature. Here, these polynomials were obtained with the help of EES software [37]. 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 pointed by Korba et al. [25], the two more accurate LBM models to deal with conjugate heat transfer are the interface treatment, proposed by Li et al. [26], and the ”enthalpy like" models, as the developed by Chen et al. [27], Rihab et al. [38]. In general, models with a 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 become considerably difficult.

The main advantage of the “enthalpy-like” models is that they do not have to perform special treatments for the interfaces, in which effects are naturally incorporated by the LBE. Then, there is no additional difficulties related to complicated geometries between the different materials. Furthermore, for the link-wise approach [25], both methods have second order accuracy when it is used the link-wise approach for the boundary nodes. Thus, in this paper, we propose a modification of the Chen et al. [27] model for accounting the temporal and local variation of 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 [39] 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 given by Eq. 12. In this equation, the term $\left(\frac{\partial ln\sout{V}}{\partial lnT}\right)\frac{DP}{dt}$ stands for the energy transfer due to compressible effects, while $\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 ln\sout{V}}{\partial lnT}\right)\frac{DP}{Dt}+\tau :\nabla \mathbf{u}+T\frac{D\rho c_P}{Dt}$$

(12)

Now, the left-hand side of Eq. 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 Eq. 13 into Eq. 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 ln\sout{V}}{\partial lnT}\right)\frac{DP}{Dt}+\tau :\nabla \mathbf{u}$$

(14)

Neglecting the viscous heat dissipation, the energy conservation equation, is given by Eq. 15. It is important to note that this equation is generic enough for being 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. [27], a new variable $h_0={\left(\rho c_P\right)}_{0}T$ is defined, being ${\left(\rho c_P\right)}_0$ a constant reference heat capacitance. Substituting $T=h_0/{\left(\rho c_P\right)}_0$ into Eq. 15 and defining $\sigma =\rho c_p/{\left(\rho c_P\right)}_0$, we obtain Eq. 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 [40], the LBE with the dimensional approach for the energy conservation equation is described by Eq. 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, given by Eq. 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 the Chapman-Enskog analysis it is possible to recover the Eq. 16 from the thermal LBE (Eq. 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, is is possible to note that, by this expansion, the LBE only recovers the the first term in the right-hand side of Eq. 16. Consequently, the other terms must be included in the LBE via the source term, which is defined by Eq. 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 to via source term to the LBE are not easy to calculate, as it involves gradient of quantities that depends on the distribution functions $g_i$. Thus, it is interesting to make some manipulations using the definition of heat flux from Karani and Huber [41], 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 a 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 Eq. 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 of 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}\Lambda \mathbf{M}\right]}_{ij}(g_j-g_j^{eq})$, where $\left[\mathbf{M}\right]$ is the transformation matrix, responsible to transform the distribution functions to the moment space and vice-versa. $[\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. [31], and the collision matrix here is defined as $[\Lambda ]=diag({\omega }_0,{\omega }_1,...,{\omega }_8)$. The relaxation frequencies can be given by ${\omega }_3^{-1}={\omega }_5^{-1}=\alpha /c_s^2+\Delta t/2$[42], 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, all the Eq. 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 Eq. 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[\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}[{\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 Eq. 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}[\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 given by Eq. 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}[\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.’s [27] enthalpy-based LBM. Our method was deducted from the conservative energy conservation equation. Through this formulation, no hypothesis regarding homogeneous mediums, or time and local independence of thermal properties, were made, and as such, both local and temporal variations of thermodynamics properties are allowed. The spatial variations enters through the source term in the thermal LBE (Eq. 17), meanwhile, the temporal variation are naturally accounted through our formulation, granting the correct capture of the transient phenomena.

3. 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 are: a 1D conduction between three different solids, a forced convection between two fluids and the natural convection in a partially-heated enclosure considering conjugate heat transfer conditions. The numerical results were compared against reference solutions, being analytical for the first two and numerical for the third, for its validation. These comparisons were quantitatively evaluated using the global relative error, defined by Eq. 26. [r]$E_R$relative error

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

(26)

3.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.60m$ is considered, being the first layer made of copper (starting from the left of the domain), the second one is silicon and the third one, aluminum. Each layer has a length of $d=0.20m$, and then the interfaces are located at $x_1=0.20m$ and $x_2=0.40m$. Both left and right wall are set with constant temperatures (Dirichlet’s boundary conditions) about $T_{wal{l}_l}=313.15K$ and $T_{wal{l}_r}=293.15K$, respectively.

Here all the solid’s properties are assumed to be constants, calculated at a temperature of $293.15K$. Being k the thermal conductivity, the properties are: for copper $k_{Cu}=401.2W{m}^{-1}{K}^{-1}$, $c_{p_{Cu}}=384.5Jk{g}^{-1}{K}^{-1}$ and ${\rho }_{Cu}=8934.0kg{m}^{-3}$; for silicon $k_{Si}=149.6W{m}^{-1}{K}^{-1}$, $c_{p_{Si}}=709.1Jk{g}^{-1}{K}^{-1}$ and ${\rho }_{Si}=2330.0kg{m}^{-3}$; and for aluminum $k_{Al}=237.0W{m}^{-1}{K}^{-1}$, $c_{p_{Al}}=905.0Jk{g}^{-1}{K}^{-1}$ and ${\rho }_{Al}=2707.0kg{m}^{-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 [43].

The results for several time steps are depicted in Figure 1, as expect the major divergences between both the analytical and the numerical solutions, occurs at the interfaces. Despite that the method 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 Eq. 26, are given in Table 1. From a quantitative point of view, both solutions present a good agreement, resulting in relatively small errors (bellow $0.05·{10}^{-2}\%$) for all time steps, thus further corroborating with the results shown in Figure 1.

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

3.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 given in Li et al. [26]. This problem consists on a channel with fixed high 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 on the form of Eq. 27 is imposed on both 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 given 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)$$

(27)

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.16114kg{m}^{-3}$, $c_{P_w}=4188.10655JK{g}^{-1}{K}^{-1}$ and ${\alpha }_w=1.61·{10}^{-7}{m}^2{s}^{-1}$ for the water, evaluated at $P_w=101325Pa$ and $T_{ref}=343.15K$; and ${\rho }_{C{O}_2}=177.45kg{m}^{-3}$, $c_{P_{C{O}_2}}=1688.66JK{g}^{-1}{K}^{-1}$ and ${\alpha }_{C{O}_2}=9.76·{10}^{-8}{m}^2{s}^{-1}$ for $C{O}_2$, evaluated at $P_{C{O}_2}=8115030Pa$ and $T_{ref}=343.15K$. 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 2 a), and the numerical and analytical temperature profiles along different cuts in the x-axis, Figure 2 b), 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 2 b), we can see that both the temperatures profiles presents a 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 (Eq. 26) between the numerical and analytical solutions. The results are given in Table 2. Overall, the results are shown to be in good agreement, with error values bellow $3·{10}^{-2}\%$. It is interesting to observe that for the initial moments the output errors were bigger, about $2·{10}^{-2}\%$. This results can be associated with a difficult to precisely match 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 a almost constant error for the remaining studied times.

3.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. [44]. It is important to effort that some modification were made to the original problem, instead of imposing heat directly into 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 simulations the Rayleigh number was maintained at $Ra={10}^4$ and the dimensionless length of the heat source were set as $ϵ=0.8$. Besides, serving the purpose of the validating the employed model, this test was performed with the intuit of choosing a reasonable mesh size to correctly predict the natural convection phenomenon that will be studied on 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 bounce-back scheme, while the Dirichlet boundary conditions were implemented via anti-bounce-back scheme. Lastly, for the Neumann boundary conditions, the scheme proposed in Martins et al. [45] were used.

First, the results were organized with reference to the mesh size in Table 3, to facilitate the comparison between different meshes. 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 meshes with relative errors bellow $3\%$, even with the bigger meshes, thus advocating for its accuracy. Additionally, Figure 4, presents both the temperature distributions and the streamlines for each mesh. Overall all simulations produced similar results, with the fluid presenting higher temperatures close to the cavity, and thus rising from this point, creating 2 recirculating vortexes.

Considering both results (Table 3 and Figure 4), 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 inferiors to $0.1\%$. Thus, it is the most reasonable choice for the simulations presented in the following Sec. 4.

4. Results for Natural Convection with Structured Cavities

The study examines natural convection under conjugate heat transfer conditions on structured cavities, analyzing 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 intensifying heat transfer rates at the solid surface. We conducted natural convection simulations using the structured cavities outlined in Moreira et al. [46]. Figure 5 provides schematic representation of the tested geometries. 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. [46]. Consequently, we aim here to investigate whether natural convection can be enhanced by the same surface structures. In fact the 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.2W{m}^{-1}{K}^{-1}$, $c_{p_{Cu}}=384.5Jk{g}^{-1}{K}^{-1}$, and ${\rho }_{Cu}=8934.0kg{m}^{-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.

For the hydrodynamic LBE, a non-slip boundary condition is adopted for all 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. [45] is utilized. Finally, considering the results presented in Section 3.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$.

4.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}W{m}^{-2}$ and ${\dot{Q}}_{wall}=1\times {10}^{6}W{m}^{-2}$, resulting in Rayleigh numbers of $Ra=3.8\times {10}^4$ and $Ra=3.8\times {10}^5$, respectively. Unlike 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 surfaces 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{N{u}_{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 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.

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}^{\circ }$) 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 6 c) and Figure 7 c). Initially, the temperature differences for the geometries are lower than those observed for the plane surface, suggesting a 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 analysis, we plotted both the two-dimensional temperature profiles and the streamlines for each geometry at a given instant. Figure 8 and Figure 9 show the plotted results for $Q_{wall}=1\times {10}^5$ at $t=0.25s$ and $Q_{wall}=1\times {10}^6$ at $t=0.05s$, respectively. Both instants are associated with periods of time in which the structured geometry is shown to augment 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 serve as confirmation of the performance shown in Figure 6 and Figure 7.

Now, considering the resultant fluid flow as shown in Figure 8 b) and Figure 9 b), we observe that the fluid flows 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}^{\circ }$). 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 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 10 correspond to regions where the geometries exhibit better performance ($t=0.25s$), a region where the geometries and the surfaces have similar performance ($t=1.00s$), and finally, regions where the plane surface surpasses the cavities ($t=2.50s$). By comparing the results in Figure 10a) with Figure 8 and Figure 6, it is evident that the higher performance of the cavities is associated with the observed vorticity, where higher levels of vorticity indicates better performance. Still considering Figure 10a), 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 10, it is observable that all 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 in 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 10 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 moving fluid, thus avoiding 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 [47,48]. Figure 11 shows the evolution of the effectiveness of the geometries over time, confirming the previous results. Geometry 3 with $\theta ={60}^{\circ }$ 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 the one imposed at the base. Thus, the steady-state effectiveness may not be the best parameter to evaluate the structured cavities.

Nonetheless, transient analysis of ${\epsilon }_f$ can provide interesting results. By analyzing the results displayed in Figure 11a), 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. As expected, Geometry 3 presented a higher estimated time constant, which again serves as an indication of its better performance.

Table 4. Time constant estimation.

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

Table 4. Time constant estimation.

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

The higher effectiveness exhibited by the geometries in Figure 11, 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.

4.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. [r]$Nu$Nusselt number [r]$\dot{Q}$heat flux $\left[W{m}^{-2}\right]$

Figure 12, Figure 13, and Figure 14 illustrate the evolution of the average Nusselt number, heat transfer rate and temperature difference at the solid-fluid interface, respectively. The behavior observed for each condition drastically differs from one another. For $T_{wall}\le {40}^{\circ }C$, all structured cavities are shown to have a negative impact on the heat transfer process. Conversely, when considering higher wall temperatures ($T_{wall}\ge {60}^{\circ }C$), the cavities begin to exhibit regions (typically at the beginning of the simulations) where the predicted Nusselt number surpasses those predicted for the plane surface. In other words, regions where an intensification of the heat transfer process is observed. This simulated behavior should be a direct consequence of $Ra$ values. For small $Ra$ the cavities increase the solid-fluid thermal resistance due to the low flow vorticiy, leading to low heat transfer rates. This behavior will be discussed further in this section.

Figure 15 and Figure 16 depict both the temperature distributions and the streamlines for $T_{wall}={80}^{\circ }C$ at $t=2.50s$ 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 15 and Figure 16 show some major differences, which may help 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 15 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.50s$ compared to the stationary regime. This indicates a 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, 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 16, where the fluid reaches temperatures as close to the solid as possible, creating regions with high thermal resistance.

Figure 17 presents the local heat transfer rate values at the solid-fluid interface. Due to the symmetry of the problem, the results are only shown until the middle section of the problem. Two distinct regions can be identified: Outside the cavities 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.

A brief examination of Figure 18 confirms that at steady state all the studied geometries indeed exhibited lower performance compared to the plane surface, resulting in reduced heat transfer rates. However, an interesting observation is noted: the elevation of the imposed base temperature, and consequently the Rayleigh number, enhances the effectiveness of the geometries. In particular, Geometry 1 nearly matches the performance of the plane surface for higher $Ra$.

5. Conclusions

In this paper, a lattice Boltzmann method model designed for simulating conjugate heat transfer between non-heterogeneous mediums is introduced. The model incorporates a forcing term into the conventional advection-diffusion LBE, circumventing challenges associated with interface modeling. To validate the accuracy of the model, simulations of benckmark problems were conducted, demonstrating its robust performance across various conditions with consistently low relative errors compared to expected solutions.

Furthermore, the proposed LBM model was employed to investigate the natural convection phenomenon within structured cavities under the influence of conjugate heat transfer. Two principal simulations were undertaken: one involving a fixed heat flux at the domain base and another with fixed base temperature at this same domain base. Unlike previous works, emphasis was placed on the transient analysis of the problem, revealing interesting results. In both studied conditions, an augmentation in the natural convection process was perceived on the structured surfaces, mainly for higher Ra, indicating that the cavities enhance the conjugated heat transfer process in transient conditions. Specifically, regarding the simulation with a fixed heat flux, formation of secondary vortices close to the surfaces was observed. These vortices can be associated with the perceived augmentation of the Nusselt number, being responsible for enhanced advection of energy.

The results were also examined with consideration of stationary regimes. In these instances, the surface with cavities displayed lower Nusselt numbers and average heat transfer rates compared to a flat plane surface, thus supporting findings from prior research [5,10]. This occurrence is attributed to fluid stagnation between the structured cavities, which introduces an additional thermal resistance between the wall and the fluid, ultimately reducing heat transfer rates.

In summary, modifications were introduced to an existing thermal LBM model for simulating conjugate heat transfer processes [27], aimed at accounting for local and temporal variations of thermophysical properties. Subsequently, the natural convection process under conjugate heat transfer conditions on structured surfaces was investigated, revealing insights relevant to their application in heat exchangers for transient cooling of electronic devices.