Coupling finite volume-lattice Boltzmann methods for advanced heat transfer simulations
Yang Zhou, Alessandro De Rosis, Alistair Revell
TL;DR
This paper introduces a new framework that combines two simulation methods to improve heat transfer modeling in complex thermal flows.
Contribution
A novel coupling framework combining FVM and LBM with improved stability and accuracy for thermal flow simulations.
Findings
The central-moments-based collision operator improves numerical stability and accuracy.
The coupling framework shows excellent numerical accuracy and convergence in benchmark and melting scenarios.
Abstract
We present a high-performance coupled framework that advances the integration of the finite volume method (FVM) and the lattice Boltzmann method (LBM) for multi-physics thermal flow simulations, including heat conduction, conjugated heat transfer, natural and forced convection, and phase change. The proposed scheme employs a central-moments-based collision operator for both velocity and temperature fields, substantially improving numerical stability and accuracy over traditional approaches within the LBM community. The reconstruction strategy, combining regularised and high-order truncated equilibrium methods, ensures smooth and accurate data exchange at FVM–LBM coupling interfaces. The implementation employs the Parallel Location and Exchange coupling library, enabling efficient and scalable communication between the FVM and LBM. Validation against standard benchmark problems and…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 1
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 2
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 36- —https://doi.org/10.13039/501100000266Engineering and Physical Sciences Research Council
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsLattice Boltzmann Simulation Studies · Advanced Numerical Methods in Computational Mathematics · Model Reduction and Neural Networks
Introduction
Multiscale thermo-fluid problems, characterised by disparate spatial and temporal scales, are widespread in practical engineering systems such as heat exchangers [1], phase-change devices and fuel cells [2, 3]. Accurately and efficiently predicting such multiscale phenomena remains a significant challenge for conventional monolithic simulation approaches, which often struggle to balance computational efficiency with physical fidelity [4]. These limitations have motivated the development of advanced coupled numerical framework for multiscale simulations, in which the broader system behaviour is captured at macroscopic level with a comprehensive but less detailed view, while the finer-scale physics are resolved using meso/microscopic methods that are imperceptible at the macroscopic resolution. Such coupling strategies enhance both physical insight and predictive capability without the prohibitive computational demands of a purely microscopic approach.
Coupling strategies for thermo-fluid simulations generally fall into two main paradigms. The first applies heterogeneous numerical methods across the entire computational domain, assigning different techniques to distinct physical fields and coupling them through macroscopic governing equations [5–8]. The second adopts a domain decomposition approach, dividing the domain into subregions, each solved by a different numerical method [9–11], with key variables (such as velocity and temperature) exchanged at the interfaces. This latter approach is particularly well-suited to multiscale problems, as it enables tailored spatial and temporal resolutions in different regions. Indeed, macroscopic solvers efficiently handle large-scale domains with broader but coarser representations, while meso/microscopic methods resolve fine-scale details within localised regions that would otherwise be lost. This synergistic integration strikes a crucial balance, achieving both the computational efficiency necessary for large-scale simulations and the physical fidelity required to capture the multiscale complexity of real-world thermo-fluid systems.
In response to their practical importance, considerable researches have focused on developing coupled numerical frameworks for multiscale thermo-fluid simulations. A number of these coupling schemes have been successfully established and validated through rigorous benchmark studies, including hybrid combinations such as molecular dynamics–finite volume method [12, 13], molecular dynamics–lattice Boltzmann method [14], direct simulation Monte Carlo–finite volume method [15], lattice Boltzmann method–finite difference method [16], and finite volume method–lattice Boltzmann method (FVM–LBM) [9, 10]. Amongst available methods, the lattice Boltzmann method (LBM) stands out as a powerful mesoscopic approach, bridging the gap between microscopic and macroscopic models for fluid flow and heat transfer. In short, the LBM computes the spatio-temporal evolution of collections (also known as populations or distributions) of fictitious particles, which carry information about the underlying flow physics. In parallel, the finite volume method (FVM) remains the leading choice for macroscopic simulations, valued for its robustness and computational efficiency in engineering applications. By combining these strengths, the coupled FVM–LBM framework offers a compelling solution for practical multiscale problems, capturing fine microscale physics with LBM while leveraging FVM’s ability to efficiently handle large-scale domains.
A central challenge in FVM–LBM coupling lies in the accurate reconstruction of LBM particle distribution functions from macroscopic variables, due to the lack of a direct mapping between the two formulations. This issue was first addressed by Xu et al. [9], who derived analytical relations between macroscopic fields and mesoscopic populations using a single-relaxation-time (SRT) LBM scheme based on the Bhatnagar-Gross-Krook (BGK) approximation. Their pioneering framework enabled coupled simulations of laminar lid-driven cavity flows up to Reynolds numbers of 1000. Building on this foundation, Luan et al. [17, 18] extended the methodology to more complex configurations, including airfoil aerodynamics and porous media transport, thereby highlighting the framework’s versatility for multiscale flow problems. Later, Tong et al. [19] applied the FVM–LBM coupling to model fouling layer formation on heat exchanger tubes, validating the approach for industrial-scale particulate flows through parametric studies on particle size and inlet velocity.
Efforts to extend coupled FVM–LBM methods for thermal flows have also seen important progress [20]. Luan et al. [21] developed the first thermal coupling algorithm based on Chapman–Enskog expansion under the SRT scheme, and validated it in two-dimensional natural convection scenarios. Tong & He [22] proposed a unified FVM–LBM framework capable of simulating coupled flow and heat transfer processes up to Rayleigh numbers (Ra) of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.4\times 10^5$$\end{document} . Salimi et al. [11] improved numerical stability by incorporating a multiple-relaxation-time (MRT) collision operator, though their reconstruction operators remained based on SRT LBM [21]. Their work also explored conjugate heat transfer over heated square cylinders with porous layers. At the cost of increased complexity, Salimi et al. [23] developed reconstruction operators tailored to the MRT framework and demonstrated successful simulations of three-dimensional natural convection at Ra = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^5$$\end{document} . Alternative approaches used D2Q9 and D2Q5 lattice models with non-equilibrium extrapolation for velocity and temperature fields, respectively [10]. However, simulations at Ra = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} still exhibited interfacial temperature discontinuities, highlighting unresolved challenges in high-Rayleigh-number coupled simulations.
Despite significant progress, three major challenges still limit current FVM–LBM coupling strategies for thermal flows. First, reconstruction algorithms often involve high computational costs, reducing the overall efficiency and practical applicability of the coupled schemes. Second, validation remains limited for complex regimes such as high-Rayleigh-number natural convection, forced convection, and conjugate heat transfer. Third, little attention has been given to large-scale parallel implementations, which are crucial for industrial applications, as highlighted by He & Tao [24].
To address these gaps, here we develop a robust two-way coupling framework integrating the industrial FVM solver \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{code\_saturne}$$\end{document} [25] and the LBM solver LUMA [26], facilitated by the Parallel Location and Exchange (PLE) library [27] for efficient data transfer. This framework provides proof-of-concept applications for multiscale thermal flow simulations, in which the mesoscopic LBM resolves the local thermal flow in one sub-domain while the macroscopic FVM is used to simulate the another sub-domain. Building upon the multiphysics groundwork by De Rosis & Coreixas [28], we advance the state-of-the-art by implementing a thermal LBM scheme based on a central-moments (CMs) collision operator. This CMs-based scheme offers enhanced stability and accuracy compared to traditional MRT and SRT models [29–32]. Additionally, we adopt regularised and high-order truncated equilibrium distribution schemes to reconstruct the particle distribution functions from macroscopic data [4, 33]. The framework is rigorously validated across several well-consolidated benchmark cases involving natural and forced convection, as well as conjugate heat transfer. Moreover, it proves particularly effective in capturing complex melting dynamics, benefiting from the complementary strengths of the FVM and LBM approaches.
The rest of this paper is organised as follows. Section 2 presents the numerical methodology underpinning the coupling framework. Section 3 evaluates its performance across a suite of benchmark problems. Section 4 concludes with key findings and outlines future research directions in FVM–LBM coupling for thermal flow simulations.
Methodology
Governing equations
Let us consider a three-dimensional Cartesian coordinate system of axes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x-y-z$$\end{document} . The macroscopic behaviour of an incompressible Newtonian fluid with thermal effects (under the Boussinesq approximation) is governed by the following set of equations:
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \displaystyle \boldsymbol{\nabla } \cdot \boldsymbol{u}= & 0 , \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} &\displaystyle \partial _t \boldsymbol{u} + \left( \boldsymbol{u} \cdot \boldsymbol{\nabla } \right) \boldsymbol{u}= -\dfrac{1}{\rho _0} \boldsymbol{\nabla } p \\&+ \nu \boldsymbol{\nabla }^2 \boldsymbol{u} +\boldsymbol{g} \beta \left( T-T_0\right) , \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\partial _t T} + \boldsymbol{u} \cdot \boldsymbol{\nabla } T= & \boldsymbol{\nabla } \cdot (\alpha \boldsymbol{\nabla } T), \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}=[u_x, u_y, u_z]$$\end{document} denotes the velocity vector, t is the time, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _0$$\end{document} is the reference mass fluid density, p indicates the pressure, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document} denotes the fluid kinematic viscosity, T denotes the temperature, the thermal diffusivity is defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha = {\lambda /( \rho c_p})$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_p$$\end{document} representing the thermal conductivity and specific heat capacity of the fluid, respectively, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_0$$\end{document} is the reference temperature, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} represents the thermal expansion coefficient. The acceleration due to gravity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{g}$$\end{document} acts in the direction y and has module \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g=-9.806$$\end{document} . Gradient and Laplacian operators are defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\nabla } = [\partial _x, \partial _y, \partial _z]$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\nabla }^2 = [\partial ^2_x+\partial ^2_y+\partial ^2_z]$$\end{document} , respectively, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial $$\end{document} representing a partial derivative.
The physics of the problem is governed by the Rayleigh number and the Prandtl number, which are defined as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textrm{Ra} = \frac{g \, \beta \, \Delta T \, L^3}{\nu \, \alpha }, \qquad \textrm{Pr} = \frac{\nu }{\alpha }, \end{aligned}$$\end{document}where L and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta T$$\end{document} are the characteristic length and temperature difference of the problem, respectively.
In the FVM implementation, the computational domain is decomposed into numerous control volumes. The governing Eqs. (1)–(3) are numerically integrated over each control volume to derive the corresponding discretised equations. For resolving the velocity-pressure coupling, the semi-implicit method for pressure-linked equation consistent [25] algorithm is employed on the FVM side.
Unlike traditional Navier–Stokes solvers, which operate at the macroscopic scale, the lattice Boltzmann method models fluid dynamics from a mesoscopic perspective. It tracks the spatio-temporal evolution of discrete particle distribution functions, whose statistical moments yield macroscopic quantities such as velocity and pressure. For thermal flow simulations, we adopt the double-distribution function (DDF) approach, which employs two distinct sets of distribution functions: one for the velocity field and another for the temperature field. This separation allows for efficient and accurate resolution of coupled thermo-fluid dynamics within the LBM framework. Let us consider the classic BGK scheme with D3Q19 spatial discretization first. The lattice Boltzmann equation (LBE) for the velocity and temperature fields are defined as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} &\displaystyle f_i(\boldsymbol{x} + \boldsymbol{c}_i \Delta t, t + \Delta t)= f_i(\boldsymbol{x}, t)\\&+ \Omega _i(\boldsymbol{x},t) + \left( 1-{\omega \over 2}\right) F_i(\boldsymbol{x}, t), \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} g_i(\boldsymbol{x} + \boldsymbol{c}_i \Delta t, t + \Delta t)&= g_i(\boldsymbol{x}, t)+ \Omega _{T,i}(\boldsymbol{x},t), \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{x}=[x, \, y, \, z]$$\end{document} represents the spatial coordinates of a lattice node, the time step \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t$$\end{document} is implicitly taken as 1. The particle distribution functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|f_i\rangle = [f_0,...,~f_{i},...,~f_{18}]^{\top }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|g_i\rangle = [g_0,...,~g_{i},...,~g_{18}]^{\top }$$\end{document} govern the evolution of the velocity and temperature fields, respectively, with the index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=0,\,1,..., \,17,\,18$$\end{document} representing the discrete velocity directions in the D3Q19 model. Here and henceforth, the Dirac notation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\bullet \rangle $$\end{document} indicates a column vector and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\top $$\end{document} represents the transpose operator. These distribution functions propagate along the prescribed lattice directions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{c}_i = \left[ |c_{ix}\rangle ,~|c_{iy}\rangle ,~|c_{iz}\rangle \right] $$\end{document} , which are defined in Eq. (11). Within the BGK approximation, the collision operators \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{T,i}$$\end{document} can be written as a relaxation of populations towards an equilibrium state,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Omega _i= \omega (f_i^{eq}-f_i), \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Omega _{T,i}= \omega _T (g_i^{eq}-g_i), \end{aligned}$$\end{document}where the superscript eq denotes the equilibrium distribution functions. The macroscopic governing Eqs. (1)–(3) are recovered through Chapman-Enskog analysis [34], where the kinematic viscosity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document} and thermal diffusivity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} are related to the relaxation frequency \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _T$$\end{document} , respectively, through the expressions
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \nu= & c_s^2\left( 1/\omega - {1 / 2}\right) , \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \alpha= & c_s^2\left( 1/\omega _T - {1/ 2}\right) , \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_s=1/\sqrt{3}$$\end{document} represents the lattice sound speed in the D3Q19 discretisation scheme. While the source term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_i$$\end{document} is commonly handled through the second-order truncated model of Guo et al. [35], in this work we instead adopt the fully Hermite-based formulation introduced by De Rosis & Coreixas [28]. In addition, the equilibrium distribution functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_i^{eq}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_i^{eq}$$\end{document} in the CMs-based LBM scheme are not approximated using the conventional second-order truncated formulation, and the incompressibility correction is not employed in the present study. In fact, exploiting the full potential of any lattice Boltzmann discretisation requires the usage of the complete allowable set of Hermite polynomials [36, 37]. Defining a tensor product notation for the discrete parameter space \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\psi ,\gamma ,\chi )\in \{\pm 1\}^3$$\end{document} and making reference to De Rosis & Coreixas [28], \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_i^{eq}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_i^{eq}$$\end{document} are rewritten in Eq. (12)
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} | c_{ix}\rangle&= [0,\, 1,\, -1,\, 0,\, 0,\, 0,\, 0,\, 1,\, \nonumber \\&-1,\, 1,\, -1,\, 1,\, -1,\, 1,\, -1,\, 0,\, 0,\, 0,\, 0 ]^{\top }, \nonumber \\ | c_{iy}\rangle&= [0,\, 0,\, 0,\, 1,\, \nonumber \\&-1,\, 0,\, 0,\, 1,\, -1,\, -1,\, 1,\, 0,\, 0,\, 0,\, 0,\, 1,\, -1,\, 1,\, -1 ]^{\top }, \nonumber \\ | c_{iz}\rangle&= [0,\, 0,\, 0,\, 0,\, 0,\, 1,\, -1,\, 0,\, 0,\, 0, \, 0,\, 1,\, \nonumber \\&-1,\, -1,\, 1,\, 1,\, -1,\, -1,\, 1 ]^{\top }. \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & f^{eq}_{(0, 0, 0)}=g^{eq}_{(0, 0, 0)}=\frac{\xi }{3} \left[ 1-(u_x^2+u_y^2+u_z^2)\right. \nonumber \\ & \left. + 3(u_x^2 u_y^2+u_x^2 u_z^2+u_y^2 u_z^2)\right] ,\nonumber \\ & f^{eq}_{(\psi , 0, 0)}=g^{eq}_{(\psi , 0, 0)}= \frac{\xi }{18} \left[ 1 +3 \psi u_x + 3 (u_x^2- u_y^2- u_z^2) \right. \nonumber \\ & -\left. 9\psi (u_x u_y^2 + u_x u_z^2)-9 (u_x^2 u_y^2 + u_x^2 u_z^2)\right] ,\nonumber \\ & f^{eq}_{(0, \lambda , 0)}=g^{eq}_{(0, \lambda , 0)}= \frac{\xi }{18} \left[ 1 +3 \lambda u_y +3 (- u_x^2+ u_y^2- u_z^2) \right. \nonumber \\ & \left. -9\lambda (u_x^2 u_y + u_y u_z^2)-9 (u_x^2 u_y^2+ u_y^2 u_z^2)\right] ,\nonumber \\ & f^{eq}_{(0, 0, \chi )}=g^{eq}_{(0, 0, \chi )}= \frac{\xi }{18} \left[ 1+3 \chi u_z +3(- u_x^2- u_y^2+ u_z^2) \right. \nonumber \\ & \left. -9\chi (u_x^2 u_z + u_y^2 u_z)-9 (u_x^2 u_z^2+ u_y^2 u_z^2)\right] , \\ & f^{eq}_{(\psi , \lambda , 0)}=g^{eq}_{(\psi , \lambda , 0)} = \frac{\xi }{36} \left[ 1+3(\psi u_x + \lambda u_y) + 3 (u_x^2+ u_y^2)+9\psi \lambda u_x u_y\right. \nonumber \\ & \left. + 9(\lambda u_x^2 u_y+\psi u_x u_y^2)+9 u_x^2 u_y^2\right] ,\nonumber \\ & f^{eq}_{(\psi , 0, \chi )}=g^{eq}_{(\psi , 0, \chi )} = \frac{\xi }{36} \left[ 1+3(\psi u_x + \chi u_z) + 3 (u_x^2+ u_z^2) \right. \nonumber \\ & \left. +9\psi \chi u_x u_z+ 9(\chi u_x^2 u_z+\psi u_x u_z^2)+9 u_x^2 u_z^2\right] ,\nonumber \\ & f^{eq}_{(0, \lambda , \chi )}=g^{eq}_{(0, \lambda , \chi )}= \frac{\xi }{36} \left[ 1+3(\lambda u_y + \chi u_z) + 3 (u_y^2+ u_z^2)\right. \nonumber \\ & \left. +9\lambda \chi u_y u_z+ 9(\chi u_y^2 u_z+\lambda u_y u_z^2)+9 u_y^2 u_z^2\right] ,\nonumber \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi $$\end{document} equals to the value of density ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} ) and temperature (T) for the velocity and temperature fields, respectively. The above equilibrium states are able to recover certain third- and fourth-order moments of the Maxwell-Boltzmann distribution. While fourth-order equilibrium moments do not affect the macroscopic behaviour but only on the numerical stability [38], the third-order moments are essential for recovering the correct viscous stress tensor [37]. Eventually, the macroscopic variables ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}$$\end{document} and T) are obtained by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \rho (\boldsymbol{x},t)= & \sum _{i=0}^{18} f_i(\boldsymbol{x}, t), \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \boldsymbol{u}(\boldsymbol{x},t)= & \frac{1}{\rho (\boldsymbol{x},t)} \left[ \sum _{i=0}^{18} \boldsymbol{c}_i f_i(\boldsymbol{x}, t) + {\boldsymbol{F}(\boldsymbol{x},t) \over 2} \right] , \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} T(\boldsymbol{x},t)= & \sum _{i=0}^{18} g_i(\boldsymbol{x}, t). \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle \boldsymbol{F}(\boldsymbol{x},t) = \rho (\boldsymbol{x},t) \, \boldsymbol{g} \, \beta \, \left[ T(\boldsymbol{x},t)-T_0 \right] $$\end{document} for buoyancy force induced by temperature difference.
Central-moments-based LBM for velocity field
When implementing the central-moments-based (CMs) LBM collision operator for fluid flow simulations, the lattice Boltzmann equation (LBE) can be expressed as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} &\displaystyle |f_i(\boldsymbol{x} + \boldsymbol{c}_i, t + 1)\rangle= |f_i(\boldsymbol{x}, t)\rangle \\&+ \boldsymbol{\Lambda } \left( |f_i^{eq}(\boldsymbol{x}, t)\rangle - |f_i(\boldsymbol{x}, t)\rangle \right) \nonumber \\&+ (\boldsymbol{\textrm{I}} - \boldsymbol{\Lambda }/2)|F_i(\boldsymbol{x}, t)\rangle , \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{I}$$\end{document} is the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$19 \times 19$$\end{document} unit tensor, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\Lambda }$$\end{document} denotes the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$19 \times 19$$\end{document} collision matrix [28, 39], and its role will be elucidated later. The term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_i$$\end{document} accounts for external body forces \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{F}=[F_x,F_y,F_z]$$\end{document} , and its prefactor accounts for discrete effects originating from the change of variables that aims at obtaining a numerical scheme explicit in time [35]. The LBE can be split into the collision and streaming steps, which are expressed as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&|f_i^\star (\boldsymbol{x}, t)\rangle = |f_i(\boldsymbol{x}, t)\rangle + \boldsymbol{\Lambda } \left( |f_i^{eq}(\boldsymbol{x}, t)\rangle - |f_i(\boldsymbol{x}, t)\rangle \right) \nonumber \\&\qquad \qquad + (\boldsymbol{\textrm{I}} - \boldsymbol{\Lambda }/2)|F_i(\boldsymbol{x}, t)\rangle , \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&|f_i(\boldsymbol{x}+\boldsymbol{c}_i, t+1)\rangle = |f_i^\star (\boldsymbol{x}, t)\rangle , \end{aligned}$$\end{document}where the superscript \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\star $$\end{document} represents the post-collision distribution functions. For a CM-based collision operator, the discrete velocities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle \boldsymbol{\overline{c}}_i = \left[ \langle \overline{c}_{ix}|, \langle \overline{c}_{iy}|, \langle \overline{c}_{iz}| \right] ^{\top }$$\end{document} are shifted by the local fluid velocity [29], which are defined as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \langle \overline{c}_{ix}|= & \langle c_{ix} - u_x|, \nonumber \\ \langle \overline{c}_{iy}|= & \langle c_{iy} - u_y|, \nonumber \\ \langle \overline{c}_{iz}|= & \langle c_{iz} - u_z|. \end{aligned}$$\end{document}where the symbol \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle \bullet |$$\end{document} represents a row vector. To apply the collision step in the CMs space, a suitable basis is used to transform CMs from populations, and transformation matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{T}$$\end{document} is given by [28]
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\textbf{T}} = \left[ \begin{array}{c} \langle |\boldsymbol{c}_i|^0| \\ \langle \bar{c}_{ix}|\\ \langle \bar{c}_{iy}| \\ \langle \bar{c}_{iz}| \\ \langle \bar{c}_{ix}^2+ \bar{c}_{iy}^2 + \bar{c}_{iz}^2 | \\ \langle \bar{c}_{ix}^2- \bar{c}_{iy}^2 | \\ \langle \bar{c}_{iy}^2- \bar{c}_{iz}^2|\\ \langle \bar{c}_{ix} \bar{c}_{iy}| \\ \langle \bar{c}_{ix} \bar{c}_{iz} | \\ \langle \bar{c}_{iy} \bar{c}_{iz}| \\ \langle \bar{c}_{ix}^2\bar{c}_{iy}|\\ \langle \bar{c}_{ix}\bar{c}_{iy}^2| \\ \langle \bar{c}_{ix}^2\bar{c}_{iz}|\\ \langle \bar{c}_{ix}\bar{c}_{iz}^2| \\ \langle \bar{c}_{iy}^2\bar{c}_{iz}| \\ \langle \bar{c}_{iy}\bar{c}_{iz}^2| \\ \langle \bar{c}_{ix}^2\bar{c}_{iy}^2| \\ \langle \bar{c}_{ix}^2\bar{c}_{iz}^2| \\ \langle \bar{c}_{iy}^2\bar{c}_{iy}^2| \end{array} \right] . \end{aligned}$$\end{document}The collision operator in the populations space is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf {\Lambda } = \mathbf {T^{-1} K T}$$\end{document} , where in the present study \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{K}=\textrm{diag}[1,~1,~1,~1,~1,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega ,~\omega ,~\omega ,~\omega ,~\omega ,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1,~1,~1,~1,~1,~1,~1,~1,~1]$$\end{document} is the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$19\times 19$$\end{document} relaxation matrix in the CM space while the omitted elements are equal to zero. Notably, if the main diagonal elements of the relaxation matrix are replaced with the relaxation frequency ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} ), the CM scheme will transform into CM-SRT scheme [40]. Substituting the relaxation matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf {\Lambda }$$\end{document} into collision process in Eq. (17), and it takes place as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} |k_i^*\rangle =(\mathbf {I-K})|k_i\rangle + \textbf{K}|k_i^{eq}\rangle + \left( \textbf{I} - {\textbf{K}\over 2} \right) |R_i\rangle . \end{aligned}$$\end{document}where the post-collision, pre-collision, equilibrium and forcing term CMs are calculated by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|k_i^{*}\rangle = \textbf{T}|f_i^*\rangle $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|k_i\rangle = \textbf{T}|f_i\rangle $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|k_i^{eq}\rangle = \textbf{T}|f_i^{eq}\rangle $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|R_i\rangle = \textbf{T}|F_i\rangle $$\end{document} , respectively. Applying the transformation matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{T}$$\end{document} to equilibrium populations in Eq. (12) generates the equilibrium CMs
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & k_0^{eq}=\rho , \nonumber \\ & k_4^{eq}=3\rho c_s^2, \nonumber \\ & k_{16}^{eq}=\rho c_s^4,\\ & k_{17}^{eq}=\rho c_s^4, \nonumber \\ & k_{18}^{eq}=\rho c_s^4,\nonumber \end{aligned}$$\end{document}while the remaining terms are equal to zero. After the collision, non-zero CMs obtain as [28]
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & k_0^*=\rho ,\nonumber \\ & k_1^*=F_x/2,\nonumber \\ & k_2^*=F_y/2,\nonumber \\ & k_3^*=F_z/2,\nonumber \\ & k_4^*=3\rho c_s^2,\nonumber \\ & k_5^*=(1-\omega ) k_5,\nonumber \\ & k_6^*=(1-\omega ) k_6,\nonumber \\ & k_7^*=(1-\omega ) k_7,\nonumber \\ & k_8^*=(1-\omega ) k_8,\\ & k_9^*=(1-\omega ) k_9,\nonumber \\ & k_{10}^*=k_{15}^*=F_y c_s^2/2,\nonumber \\ & k_{11}^*=k_{13}^*=F_x c_s^2/2,\nonumber \\ & k_{12}^*=k_{14}^*=F_z c_s^2/2,\nonumber \\ & k_{16}^*=\rho c_s^4,\nonumber \\ & k_{17}^*=\rho c_s^4,\nonumber \\ & k_{18}^*=\rho c_s^4,\nonumber \end{aligned}$$\end{document}where pre-collision CMs are given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & k_5=\sum _i^{18}f_i(\overline{c}_{ix}^2-\overline{c}_{iy}^2), \nonumber \\ & k_6=\sum _i^{18}f_i(\overline{c}_{iy}^2-\overline{c}_{iz}^2), \nonumber \\ & k_7=\sum _i^{18}f_i \overline{c}_{ix} \overline{c}_{iy}, \\ & k_8=\sum _i^{18}f_i \overline{c}_{ix} \overline{c}_{iz}, \nonumber \\ & k_9=\sum _i^{18}f_i \overline{c}_{iy} \overline{c}_{iz}. \nonumber \end{aligned}$$\end{document}The post-collision populations can be obtained by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} |f^*_i\rangle =\textbf{T}^{-1}|k_i^*\rangle . \end{aligned}$$\end{document}According to the Refs. [41, 42], the two-step approach is adopted to compute the post-collision distribution function from central moments. Eventually, the populations are streamed following the directions in Eq. (11) and the macroscopic density and velocity are obtained by Eqs. (13)–(14).
Central-moments-based LBM for temperature field
Within the CMs-based LBM framework, the LBE for the convection-diffusion process in Eq. (3) without heat source contributions can be formulated as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} |g_i(\boldsymbol{x}+\boldsymbol{c}_i, t+1)\rangle = |g_i(\boldsymbol{x},t)\rangle + \boldsymbol{\Lambda }_T \left( |g_i^{eq}(\boldsymbol{x},t)\rangle -|g_i(\boldsymbol{x},t)\rangle \right) , \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\Lambda }_T$$\end{document} represents the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$19 \times 19$$\end{document} collision matrix for temperature distribution functions. This thermal LBE is decomposed into the collision and streaming steps as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & |g_i^{\star }(\boldsymbol{x},t)\rangle = |g_i(\boldsymbol{x},t)\rangle + \boldsymbol{\Lambda }_T \left( |g_i^{eq}(\boldsymbol{x},t)\rangle -|g_i(\boldsymbol{x},t)\rangle \right) , \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & |g_i(\boldsymbol{x}+\boldsymbol{c}_i, t+1)\rangle = |g_i^{\star }(\boldsymbol{x},t)\rangle . \end{aligned}$$\end{document}For the thermal CM-based scheme, the lattice directions are also shifted by the local fluid velocity through the Eqs. (19). In addition, the collision matrix is formulated as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\Lambda }_T = \textbf{T}^{-1} \textbf{K}_T \textbf{T}$$\end{document} , where the same transformation matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{T}$$\end{document} in Eq. (20) is applied to the collision step. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{K}_T$$\end{document} represents the diagonal matrix of relaxation rates in temperature field, which is defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{K}_T=\textrm{diag}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[1,~\omega _T,~\omega _T,~\omega _T,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1,~1,~1,~1,~1,~1,~1,~1,~1,$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$~1,~1,~1,~1,~1,~1]$$\end{document} , and the off-diagonal elements remain as zero.
The collision operation is executed in the space of CMs, and it is expressed as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\begin{matrix} |k_{i,T}^{\star }\rangle =& \left( \textbf{I}- \textbf{K}_{T} \right) \textbf{T} |g_i\rangle + \textbf{K}_{T} \textbf{T} |g_i^{\textrm{eq}}\rangle \\ =& \left( \textbf{I}- \textbf{K}_{T} \right) |k_{i,T}\rangle + \textbf{K}_{T} |k_{i_,T}^{\textrm{eq}}\rangle , \end{matrix}} \end{aligned}$$\end{document}where the vectors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|k_{i,T}\rangle $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|k_{i,T}^{\textrm{eq}}\rangle $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|k_{i,T}^{\star }\rangle $$\end{document} represent the collections of pre-collision, equilibrium and post-collision moments in temperature field, respectively. Following the definition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|k_{i,T}^{\textrm{eq}}\rangle = \textbf{T} |g_i^{\textrm{eq}}\rangle $$\end{document} , only five equilibrium moments are non-zero, i.e.,
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} k_{0,T}^{eq}= & T, \nonumber \\ k_{4,T}^{eq}= & T, \nonumber \\ k_{16,T}^{eq}= & T c_{s}^4, \nonumber \\ k_{17,T}^{eq}= & T c_{s}^4, \nonumber \\ k_{18,T}^{eq}= & T c_{s}^4. \end{aligned}$$\end{document}Moreover, the equilibrium moments exhibit both elegant mathematical form and Galilean invariance, that is, they remain independent of the reference frame velocity. Such findings align with the theoretical analysis carried out by De Rosis & Luo [31], where they demonstrated that Galilean invariant equilibrium central moments can be obtained only when the transformation matrix is applied to equilibrium populations equipped with the complete basis of Hermite polynomials. The non-zero post-collision central moments are given by
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} k_{0,T}^{\star }= & T, \nonumber \\ k_{4,T}^{\star }= & T, \nonumber \\ k_{1,T}^{\star }= & \left( 1-\omega _{T} \right) k_{1,T}, \nonumber \\ k_{2,T}^{\star }= & \left( 1-\omega _{T} \right) k_{2,T}, \nonumber \\ k_{3,T}^{\star }= & \left( 1-\omega _{T} \right) k_{3,T}, \nonumber \\ k_{16,T}^{\star }= & T c_{s}^4,\nonumber \\ k_{17,T}^{\star }= & T c_{s}^4,\nonumber \\ k_{18,T}^{\star }= & T c_{s}^4, \end{aligned}$$\end{document}where the required pre-collision central moments are
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} k_{1,T}= & \sum _{i=0}^{18} g_i \bar{c}_{ix}, \nonumber \\ k_{2,T}= & \sum _{i=0}^{18} g_i \bar{c}_{iy}, \nonumber \\ k_{3,T}= & \sum _{i=0}^{18} g_i \bar{c}_{iz}. \end{aligned}$$\end{document}The post-collision distribution functions are calculated through
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} |g_i^{\star }\rangle = \textbf{T}^{-1} |k_{i, T}^{\star } \rangle , \end{aligned}$$\end{document}and streamed by Eq. (28). Similar to the works [41, 42], the two-step approach is employed to calculated the post-collision distribution functions from central moments, and the interested reader can find the resultant expression in Appendix. Eventually, the macroscopic temperature is computed by the Eq. (15).
For the FVM solver using code_saturne, identical boundary type is applied to the velocity components due to the coupled treatment of velocity solver. Besides, the constant and adiabatic thermal boundary conditions are employed in the temperature field [27]. On the LBM side solved by LUMA, the boundary conditions for the velocity and temperature fields are imposed as follows: the no-slip and adiabatic wall conditions via the bounce-back scheme [34, 43], which defines as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_{\bar{i}}(\boldsymbol{x}, t+\Delta t)= f_i^{\star }(\boldsymbol{x}, t)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g_{\bar{i}}(\boldsymbol{x}, t+\Delta t)= g_i^{\star }(\boldsymbol{x}, t)$$\end{document} where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{i}$$\end{document} is the reflected direction of i; velocity inlets using the regularised boundary method [33], i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_i(\boldsymbol{x}, t+\Delta t)=f_i^{eq}(\boldsymbol{x}, t+\Delta t)+(w_i \boldsymbol{Q:\Pi }^{(1)})/2c_s^4 $$\end{document} , and constant temperature following Ref. [44].
Coupling process between FVM and LBM
In our framework, the computational domain is partitioned into subregions, independently solved by the FVM and LBM. As shown in Fig. 1, both solvers operate on grids with identical spatial resolution: the FVM computes the solution over the solid red grid region, while the LBM operates within the blue dashed lattice region. A grey overlapping buffer zone, where both methods coexist, facilitates a smooth and stable coupling interface. Coupling is enforced at designated boundaries by exchanging macroscopic physical quantities such as velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}$$\end{document} and temperature T. The process begins by establishing a geometric correspondence between interface points, matching spatial coordinates across the FVM and LBM domains. Once identified, a two-way data exchange is performed, with macroscopic variables transferred to the corresponding solver interfaces, as indicated by the black arrows in Fig. 1. Dirichlet boundary conditions are imposed on both sides to ensure continuity. For parallel execution, the domain is further divided among multiple MPI ranks, shown in yellow, green, blue, and purple. The entire coupling procedure—including dynamic mapping generation, data exchange, and synchronisation at each coupling step—is managed by the PLE library [27], ensuring efficient and scalable communication between the FVM and LBM solvers.Fig. 1. Spatial coupling between the FVM and LBM. The yellow, green, blue and purple regions represent MPI ranks
The identical grid density across solvers enables a direct one-to-one mapping between FVM cell centres and LBM lattice nodes at the coupling interface. However, due to the staggered arrangement of the grids, FVM interface points correspond to interstitial positions relative to the LBM lattice. To ensure accurate LBM-to-FVM information transfer, a second-order Taylor-series-based interpolation scheme is employed, reconstructing macroscopic variables at the FVM coupling interface using data from neighbouring LBM nodes. The required spatial derivatives for this interpolation are computed using the gradient formulation [45]
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \boldsymbol{\nabla }\phi \left( \boldsymbol{x} \right) = \frac{1}{c_s^2} \sum _i^{18} w_i \boldsymbol{c}_i \phi \left( \boldsymbol{x} +\boldsymbol{c}_i \right) , \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} is a certain physical quantity to be interpolated. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_i$$\end{document} are the weights of the selected lattice discretisation, which are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_0=1/3$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_{1-6}=1/18$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w_{7-18}=1/36$$\end{document} .
Macroscopic variables obtained from Eqs. (13)–(15) are directly transferred from the LBM to FVM after a suitable unit conversion. In contrast, when transferring velocity and temperature from FVM to LBM, these macroscopic quantities must be transformed into the corresponding unknown distribution functions (namely \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_8$$\end{document} as shown in Fig. 1). Here, we reconstruct these by implementing the regularized boundary approach [33] to determine the distribution functions of the velocity field, which are expressed as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & f_i = f_i^{eq} + \frac{w_i}{2c_s^4} \boldsymbol{Q}: \boldsymbol{\Pi }^{(1)}, \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle \boldsymbol{Q} = \boldsymbol{c}_i \boldsymbol{c}_i - c_s^2\textbf{I}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle \boldsymbol{\Pi }^{(1)} = \sum _i \boldsymbol{c}_i \boldsymbol{c}_i \left( f_i -f_i^{eq} \right) $$\end{document} . At the coupling interface, the unknown distribution functions for the temperature field are determined using the extended equilibrium distribution functions, i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle g_i = g_i^{eq}$$\end{document} as specified in Eq. (12).
Beyond spatial coupling, careful treatment of the temporal coupling scheme is equally critical. Standard stability requirements must be satisfied by both solvers: a Courant number of unity for FVM and a lattice velocity much smaller than the speed of sound for LBM. Furthermore, while FVM simulates incompressible flows using an implicit time integration scheme, LBM adopts an explicit one, leading to inherently larger time steps for FVM compared to LBM. This temporal mismatch requires a specialised coupling strategy, illustrated in Fig. 2. The coupling process unfolds in two stages: first, interface points are identified across FVM and LBM domains during the initialisation phase, establishing a mapping relationship and performing an initial data exchange. Then, during the simulation, each FVM time step encompasses multiple LBM sub-iterations ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t_{\textrm{FVM}} = n\Delta t_{\textrm{LBM}}$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n \ge 1$$\end{document} ). Throughout these LBM sub-iterations, the coupling interface information remains fixed, preserving the transformation data from the previous FVM-LBM exchange.Fig. 2. Temporal coupling process between FVM and LBM
In the present work, we implement a coupled computational framework combining the FVM through code_saturne [25] and the LBM via LUMA [26], integrated using the PLE coupling library [46]. On the LBM side, the LUMA solver employs a two-array streaming algorithm to advance the distribution functions. To show the coupling procedure, consider a minimal configuration with two MPI ranks, each assigned to each solver. From an algorithmic perspective, each coupling time step executes the following sequence of operations:
- update macroscopic variables for both LBM by Eqs. (13), (14) and (15) and FVM;
- compute the gradient fields on the LBM side using neighbouring points as indicated by the arrows in Fig. 1 through Eq. (34).
- predict macroscopic variables at coupling points using the second-order Taylor series on the LBM side;
- exchange macroscopic variables between solvers via the PLE coupling library;
- reconstruct unknown populations through the regularised scheme and the equilibrium distribution for velocity and temperature fields on the LBM coupling interface by Eqs. (35) and (12), respectively;
- advance one FVM and n LBM time steps, and return to action 1. Additionally, Algorithm 1 presents a pseudo code for the spatial FVM-LBM coupling interface treatment, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {F_{FVM}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {F_{LBM}}$$\end{document} denote the coupling interface points on the FVM and LBM sides, while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{\textrm{FVM}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega _{\textrm{LBM}}$$\end{document} represent the coupling interface points mapping to FVM and LBM subdomains.
Algorithm 1Pseudo code for the FVM-LBM coupling interface treatment.
Results
To assess the performance, robustness, and versatility of the proposed coupling framework, we validate it against a suite of seven benchmark problems designed to progressively challenge different physical and numerical aspects of the method.
- One-dimensional heat conduction: tests the performance of the coupling framework in purely diffusive regimes, including the sensitive analysis of grid resolution, overlapping size, number of LBM sub-iterations, and varying material properties across the FVM-LBM boundary.
- Two-dimensional natural convection in a cavity: examines the ability to resolve buoyancy-driven flows at Rayleigh numbers of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} , with particular attention to thermal stability and mass conservation across the coupling interfaces.
- Two-dimensional natural convection in a side-open cavity with porous media: evaluates heat and momentum transport in complex solid–fluid configurations, a common feature in energy systems such as solar receivers and cooling devices.
- Two-dimensional normal plate velocity: introduces forced convection with vertical inflow and a moving upper boundary. This setup not only provides an analytical solution for validation but also lays the groundwork for future extensions involving fully moving geometries.
- Two-dimensional thermal lid-driven cavity: combines shear-driven flow and thermal gradients, testing the framework’s capacity to capture secondary recirculation effects and mixed convection regimes.
- Three-dimensional natural convection in a cavity: extends the two-dimensional setup to full 3D and evaluates the framework’s capability to handle large-scale, fully coupled simulations in realistic geometries.
- Rayleigh-Bénard convection with a melting boundary: demonstrates the framework’s applicability to evolving boundary problems and multiphase thermofluid dynamics, employing the hybrid strengths of FVM and LBM for tracking moving solid–liquid interfaces. Together, these cases provide a comprehensive validation across a broad spectrum of physical regimes and modelling challenges, supporting the applicability of the framework to industrial and scientific problems.
One-dimensional heat conduction
Figure 3 presents a heat conduction problem designed to evaluate the numerical accuracy of the coupling model. The computational domain of length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L = 1$$\end{document} is divided into LBM (blue) and FVM (red) subregions, occupying 0.5L and 0.54L in the x-direction, respectively. A solid region of length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_s = 0.2L$$\end{document} is placed to the left of the LBM domain, while the remaining portion is occupied by a fluid. An overlapping region of length \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_o = 0.04L$$\end{document} is shared between the two solvers to ensure a smooth coupling interface. Since the temperature distribution is independent of the y-direction, only five grid points are used in that direction.Fig. 3. One-dimensional heat conduction: sketch of the problem setup
Constant high ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_h$$\end{document} ) and low ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c$$\end{document} ) temperatures are imposed at the left and right boundaries, respectively, while periodic boundary conditions are applied elsewhere. The initial temperature is set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_0 = (T_h + T_c)/2$$\end{document} . To investigate the coupling performance across varying material properties, the thermal conductivity ratio between solid and fluid, defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa = \lambda _s / \lambda _f$$\end{document} , is varied as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa = 0.1$$\end{document} , 1, 10, and 100.
According to Eq. (3), the steady-state solution for heat conduction without a source term satisfies:
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \boldsymbol{\nabla } \cdot \left[ \alpha _{(s,f)} \boldsymbol{\nabla } T \right] = 0, \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _s$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _f$$\end{document} denote the thermal diffusivities of the solid and fluid regions, respectively. For convenience, we introduce the dimensionless temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta $$\end{document} as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta = \frac{T-T_c}{T_h-T_c}. \end{aligned}$$\end{document}Eq. (36) admits analytical solution in the form
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \theta (x) = {\left\{ \begin{array}{ll} a_s x + 1, & \text {for} \quad 0 \le x \le L_s,\\[8pt] a_f x + b_f, & \text {for} \quad L_s < x \le L, \end{array}\right. } $$\end{document}where
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ a_s = -\frac{1}{L_s + (L - L_s)\kappa }, \quad a_f = \kappa a_s, \quad b_f = 1 + a_s L_s (1 - \kappa ). $$\end{document}For a more quantitative assessment of the accuracy of the proposed methodology, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l^2$$\end{document} -norm of the relative error between numerical ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\sigma }_{\textrm{cp}}$$\end{document} ) and analytical ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\sigma }_{\textrm{an}}$$\end{document} ) solutions is computed as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \epsilon = \frac{| \boldsymbol{\sigma }_{\textrm{an}} - \boldsymbol{\sigma }_{\textrm{cp}} |}{| \boldsymbol{\sigma }_{\textrm{an}} |} \times 100\% \end{aligned}$$\end{document}Figure 4 illustrates the variation of the relative error at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa =1$$\end{document} as a function of grid resolution ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1/50,\,1/100,\,1/200,\,1/400,\,1/500$$\end{document} and 1/600) for the pure FVM, pure LBM, and coupling models. The time steps are set to 0.002 for the LBM and 0.02 for the FVM in all simulations, with ten LBM sub-iterations performed within each FVM time step in the coupling framework. The relative error decreases monotonically with increasing grid resolution for all the considered models, with the FVM showing superior performance. In the coupled simulation, the relative error decreases to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.48~\%$$\end{document} at a resolution of 1/500, which is therefore adopted for the subsequent studies.Fig. 4. One-dimensional heat conduction: grid independence study
The influence of the overlapping domain is further investigated using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_o=4,\,8,\,12,\,16,\,20$$\end{document} and 30 grid points. Table 1 shows the relative error decreases continuously as the overlapping size increases. The behaviour arises because the equilibrium distribution functions are employed as the reconstruction operator on the LBM coupling interface, and they depend on the local velocity and temperature fields.Table 1. One-dimensional heat conduction: percentage relative discrepancy for different overlapping size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_o$$\end{document} (grid points)4812162030 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \, [\%]$$\end{document} 2.231.110.760.580.480.34
Additionally, the study of number of LBM sub-iteration (n) in Table 2 shows that the relative error gradually decreases as n increases. For a fixed grid resolution, increasing n reduces the LBM time step, which in turn decreases the thermal diffusivity in lattice unit ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} ) and increases the relaxation frequency, thereby further improving the numerical accuracy [34].Table 2. One-dimensional heat conduction: percentage relative discrepancy for different LBM sub-iteration (n)n1248101620 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \, [\%]$$\end{document} 4.002.131.110.590.480.320.26
Figure 5 compares the dimensionless temperature profiles obtained from the coupled model and the analytical solution for different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} . Our results show excellent agreement with the theoretical prediction across all cases. In the overlapping region, the FVM and LBM yield identical temperature fields, and a smooth transition is observed across the coupling interface.Fig. 5. One-dimensional heat conduction: dimensionless temperature distribution for different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} , i.e. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa = 0.1$$\end{document} (black), 1 (red), 10 (green), and 100 (blue). Present results and analytical predictions are represented by lines and circles, respectively
To quantitatively assess the conjugated heat transfer performance, Table 3 reports the relative discrepancy ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document} ) at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa = 0.1$$\end{document} , 1, 10, and 100. The very slight mismatch confirms the high accuracy and robustness of our thermal coupling model.Table 3. One-dimensional heat conduction: percentage relative discrepancy between present results and analytical predictions for different \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} 0.1110100 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \, [\%]$$\end{document} 1.230.480.370.35
Two-dimensional natural convection in a cavity
The present coupling framework is further validated against the two-dimensional natural convection problem. Figure 6 sketches the computational domain of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L\times L$$\end{document} , where the blue and red sub-regions of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.5L\times L$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.54L\times L$$\end{document} are simulated by the LBM and FVM, respectively, with an overlapping domain ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.04L\times L$$\end{document} ). The entire domain adopts the same uniform resolution, i.e., 1/500.Fig. 6. Two-dimensional natural convection in a cavity: sketch of the problem setup
Constant high ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _h=1$$\end{document} ) and low temperatures ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _c=0$$\end{document} ) are enforced on the left and right borders, while the top and bottom ones use the adiabatic condition. No-slip boundary conditions are enforced at each side. The fluid density is set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho =1.0$$\end{document} with the kinetic viscosity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu = 2.0\times 10^{-4}$$\end{document} . At the beginning of the simulation, a quiescent fluid of temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _0 = 0.5$$\end{document} is considered. For each FVM time step ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t = 0.02$$\end{document} ), the LBM performs 10 sub-iterations with time step of 0.002. Gravity is enforced along the y direction and set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g=-9.806$$\end{document} . Values of the Rayleigh number of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} are investigated in this study. The Prandtl number is set to Pr \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=0.71$$\end{document} .
Figure 7 shows the velocity field with streamlines at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} . A smooth flow field can be observed throughout the computational domain, especially near the coupling interfaces. The heated fluid ascends along the left wall due to the buoyancy force, impinges on the top face, travels toward the cold wall, and finally sinks near the cooling wall. Meanwhile, a reverse flow can be found from the cooling to the heating surface, leading to a steady clockwise rotational flow in the cavity. Finally, two independent vortices are generated in the central domain at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} , while two vortices move toward the heating and cooling walls and the third vortex occurs in the centre of domain at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document} . Using traditional pure FVM method, Chen et al. [47] and Luan et al. [21] numerically investigated two-dimensional natural convection at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} leveraging commercial software FLUENT. These findings are in excellent agreement with previous efforts [21, 47].Fig. 7. Two-dimensional natural convection in a cavity: streamline of the velocity field at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$= 10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document}
Figure 8 displays the temperature field at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} , demonstrating that the isotherms smoothly transition the coupling interfaces. Their distributions match well those shown in other studies [21, 47]. Notably, at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document} , Li et al. [10] observed visible discontinuities at the FVM–LBM interfaces, which can be attributed to inaccuracies introduced by the non-equilibrium extrapolation method used for reconstructing flow and thermal fields (see Figure 16 in Ref. [10]). In contrast, our coupling approach eliminates these discontinuities, ensuring smooth temperature profiles across the interface.Fig. 8. Two-dimensional natural convection in a cavity: temperature field with isotherms at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document}
To quantify the accuracy of coupling model, Fig. 9 shows the distributions of dimensionless velocity components along \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y/L=0.5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x/L=0.5$$\end{document} at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} . The velocity component ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_x$$\end{document} ) is normalised by the maximum value ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{x,max}$$\end{document} ) on the midline \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x/L=0.5$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_y$$\end{document} is normalised by the maximum value ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{y,max}$$\end{document} ) on the midline \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y/L=0.5$$\end{document} . Note that the midline of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y/L=0.5$$\end{document} crosses the LBM, FVM domains and the coupling interfaces while the whole line of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x/L=0.5$$\end{document} is in the FVM domain. Similarly, the velocity profiles obtained from the coupling framework are in excellent agreement with the references [21, 47].Fig. 9. Two-dimensional natural convection in a cavity: velocity distributions along the midlines \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x/L=0.5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y/L=0.5$$\end{document} at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document}
In Fig. 10, the dimensionless temperature distributions are quantitatively compared with the references [21, 47] on the midline of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y/L=0.5$$\end{document} . The temperature profiles change rapidly close to the heating and cooling surfaces, while they are relatively flat towards the centre of computational domain. They remain in excellent agreement with the reference data throughout. In addition, the temperature profiles are entirely continuous across the coupling interfaces, demonstrating the robustness of the coupling implementation.Fig. 10. Two-dimensional natural convection in a cavity: temperature distribution along the midline \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y/L=0.5$$\end{document} at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document}
To quantitatively assess heat transfer performance, Davis [48] employed second-order finite difference method to study two-dimensional natural convection and analysed the Nusselt numbers. The local ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{Nu}$$\end{document} ) and averaged ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\textrm{Nu}}$$\end{document} ) Nusselt numbers are calculated on the heating wall and they are defined as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \textrm{Nu} = - {L\over \Delta T} \left( {\partial T \over \partial x}\right) , \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} & \overline{\textrm{Nu}} = \int _0^L \textrm{Nu} \, dy. \end{aligned}$$\end{document}Tables 4 and 5 list the averaged, maximum and minimum values of the Nusselt number and the corresponding location on the heating surface at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} . The five parameters obtained from the coupling model match well with those in Refs. [21, 48].Table 4. Comparison of Nusselt number on the heating wall at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} PresentDavis [48]Luan et al. [21] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\textrm{Nu}}$$\end{document} 4.5054.5104.507 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Nu_{max}$$\end{document} 7.7427.7617.738 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y/L)_{max}$$\end{document} 0.0820.0850.083 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Nu_{min}$$\end{document} 0.7000.7360.746 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y/L)_{min}$$\end{document} 1.0001.0000.997Table 5Comparison of Nusselt number on the heating wall at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document} PresentDavis [48]Luan et al. [21] \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\overline{\textrm{Nu}}$$\end{document} 8.8038.9288.807 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Nu_{max}$$\end{document} 17.68618.07617.71 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y/L)_{max}$$\end{document} 0.03750.04560.0375 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Nu_{min}$$\end{document} 0.9961.0050.978 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(y/L)_{min}$$\end{document} 1.0001.0000.997
At Ra= \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} , mass conservation is further assessed at the coupling interfaces by monitoring the relative error of the mass flux ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho u_x$$\end{document} , see Eq. 38) as a function of time, as shown in Fig. 11. Referring to Fig. 6, the left coupling interface corresponds to data transfer from LBM to FVM, while the right interface represents transfer from FVM to LBM. The results show that the error converges to approximately \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.09~\%$$\end{document} at the left interface, while it becomes negligible at the right interface. The slightly higher error at the left interface is attributed to the fact that LBM data are stored at face locations, whereas the monitored quantities are evaluated at cell centres.Fig. 11. Two-dimensional natural convection in a cavity: change of relative error of mass flux with time at the coupling interfaces
After convergence is reached, the dimensionless mass flux is compared between the coupling and pure FVM models at the left and right interfaces, with the pure FVM data extracted at the corresponding interface locations. The results indicate good agreement between the two methods, with relative errors ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} ) of only \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.03\%$$\end{document} at both interfaces, demonstrating the accuracy of mass conservation in the coupled framework (Fig. 12).Fig. 12. Two-dimensional natural convection in a cavity: comparison of normalized mass flux at the left and right interfaces between the coupling and pure FVM models
Furthermore, the relative error of mass flux between the coupling and pure FVM models is examined for different numbers of LBM sub-iterations per FVM time step, with the results summarized in Table 6. For all five cases, the grid resolution over the entire domain is fixed at 1/500, and the FVM time step is maintained at 0.02. The results show that the errors decrease monotonically with increasing n at both interfaces, and similar error levels are observed at left and right interface of each case. This trend can be attributed to the weak compressibility inherent in the LBM, as increasing n effectively reduces the LBM time step and lattice velocity which is directly related to compressibility-induced errors.Table 6. Two-dimensional natural convection in a cavity: relative error of mass flux between the coupling and pure FVM models \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t_{FVM}/\Delta t_{LBM}$$\end{document} (n)Left interfaceRight interface49.25%8.98%82.82%2.79%102.03%2.03%161.19%1.20%201.02%1.03%
The benefit of introducing the CMs-LBM scheme is evaluated with reference to the steady-state convergence for both cases Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} . The velocity and temperature residuals are defined in Eqs. (41) and (42) [47] as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathrm {Velocity\ residual} = \frac{ \left\| \textbf{u}(t+\delta t) - \textbf{u}(t) \right\| _2 }{ \left\| \textbf{u}(t+\delta t) \right\| _2 }, \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathrm {Temperature\ residual} = \sum _{x,y} \left| \frac{ T(x,y,t+\delta t) - T(x,y,t) }{ T(x,y,t+\delta t) } \right| , \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert \cdot \Vert _2$$\end{document} is the discrete \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l^2$$\end{document} norm over the domain. Figure 13 compares the velocity residual variations with dimensionless physical time between the BGK and CMs coupling schemes at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} . At Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} , the same variation is observed at initial stage ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t < 1100$$\end{document} ). Then the residual remains at about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-9}$$\end{document} for BGK scheme while the CMs scheme captures further decreases. A similar phenomenon is observed at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document} . The velocity residual starts to diverge at about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=600$$\end{document} , after which the BGK scheme remains around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-8}$$\end{document} , whereas the CM schemes exhibit a continuous decrease.Fig. 13. Two-dimensional natural convection in a cavity: variations of velocity residuals with time steps for BGK- and CMs-LBM coupling schemes at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document}
Figure 14 depicts the time history of temperature residuals based on thermal BGK and CMs coupling schemes. Using the BGK scheme, the coupling model converges to approximately \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-9}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-10}$$\end{document} at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} , respectively, while the temperature residuals continue to decrease for the CMs scheme.Fig. 14. Two-dimensional natural convection in a cavity: variation of temperature residuals with time steps for BGK- and CMs-LBM coupling schemes at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document}
According to the analysis from the Figs. 13 and 14, the residuals from CMs scheme are lower than those of the BGK scheme in both flow and temperature fields, indicating the superior numerical convergence and stability of the CMs scheme in coupling model.
At Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document} , Fig. 15 shows the variation of the velocity residual with dimensionless physical time for the pure FVM, pure LBM and the coupled models. The discontinuity in the pure FVM result and the abrupt jump in the coupling model indicate that the velocity residual reaches zero for the FVM calculation. The excellent convergence performance of pure LBM ascribes to its smaller time step, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t_{LBM}=0.002$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t_{FVM}=0.02$$\end{document} . When t is less than approximately 160, the residual of coupling model is comparable to that of the pure FVM, but remains higher than that of pure LBM model. Once the residual from the pure FVM calculation drops to zero, the residual of coupled model approaches that of pure LBM. These observations demonstrate that the convergence behaviour of the coupling model lies between those of the pure FVM and LBM methods in the flow field.Fig. 15. Two-dimensional natural convection in a cavity: variations of velocity residuals with time for the pure FVM, pure LBM and coupling models at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document}
The comparison of the temperature residual is illustrated in Fig. 16 for the pure FVM, pure LBM and the coupling models at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document} . The results indicate that the coupling model exhibits a faster initial convergence speed than both pure FVM and pure LBM models. Subsequently, the residual of the coupling model falls between those of the pure FVM and pure LBM simulations, demonstrating that its convergence performance lies between the pure FVM and pure LBM approaches.Fig. 16. Two-dimensional natural convection in a cavity: variations of temperature residuals with time for the pure FVM, pure LBM and coupling models at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document}
Two-dimensional natural convection in a side open cavity with porous media
Natural convection in a side open cavi ty with porous media has ubiquitous engineering applications, such as solar receivers, ventilation, cooling of electronic [49], which involves complex flow and solid–fluid heat transfer processes. Referring to the two-dimensional study in Ref. [50], the proposed coupling model is validated through the thermo-hydraulic effects of discrete solid blocks on the natural convection inside a fluid-filled, horizontally heated cavity with the opening to one side.
The configuration of entire computational domain is shown in Fig. 17a. The square cavity ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L\times L$$\end{document} ) is simulated by LBM while the thermal reservoir in red domain is solved by FVM ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L\times 3L$$\end{document} ), with the shadow overlapping domain ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.04L\times L$$\end{document} ) solved by both methods. There are nine conducting, disconnected and uniformly distributed square blocks and the side length is set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D=0.2L$$\end{document} . Similarly to Ref. [50], the left boundary of the cavity is defined as a heating wall ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _h=1.0$$\end{document} ) while the remaining boundaries adopt an adiabatic boundary condition. In the thermal reservoir, the top and bottom borders enforce a constant lower temperature ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _c=0.0$$\end{document} ) with the zero-gradient velocity boundary while the right boundary is set to symmetry. Figure 17b illustrates the grid points distribution over the entire computational domain. The LBM portion is decomposed into uniform lattice points, with 500 points for the length of L. Each block is captured by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$100\times 100$$\end{document} lattice points, so the porosity equals to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi =0.64$$\end{document} in the cavity. A non-uniform grid distribution is applied to the FVM domain with the mesh size ranging from 0.002 to 0.02. The overlapping domain consists of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20\times 500$$\end{document} grids along the x and y directions. The temperature field is initialised as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _c=0.0$$\end{document} everywhere. The Rayleigh number is set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle Ra = g \beta \Delta T L^3/(\nu \alpha _f)=10^5$$\end{document} , where the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _f$$\end{document} being thermal diffusivity and kinematic viscosity of the fluids. The Prandtl number ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Pr=\nu /\alpha _f$$\end{document} ) is set to 1.0. Different ratios of solid–fluid thermal diffusivity ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa =\alpha _s/\alpha _f$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _s$$\end{document} being thermal diffusivity of solid) are investigated, i.e, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa = 0.1$$\end{document} , 1.0, 10 and 100, respectively.Fig. 17. Two-dimensional natural convection in a side open cavity with porous media: geometrical configuration a and grid points distribution b
Figure 18 exhibits the temperature field closed to the open cavity with the solid–fluid conductivity ratios ranging from 0.1, 1.0, 10 to 100. The temperature difference of two adjacent isotherms is set to 0.05 in Fig. 18. At \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa =0.1$$\end{document} , the dense isotherms occur in solid blocks due to the lower thermal diffusivity. As \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} increases, the effects of heat diffusion are gradually enhanced in blocks while the temperature distributions in the fluid are induced by the convection and diffusion effects. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} increases to 100, the isotherms concentrate on the fluid, meaning that the thermal resistance of the solid is smaller. Moreover, the isotherms smoothly cross the coupling interfaces which proves the correctness of coupling scheme.Fig. 18. Two-dimensional natural convection in a side open cavity with porous media in a porous medium: isotherm distributions closed to the cavity at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa =0.1$$\end{document} , 1, 10 and 100
Additionally, Fig. 19 shows the local velocity field depicted by streamlines at different \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} . The fluid enters the cavity through the bottom gaps between blocks or blocks and wall. Next to the heating wall, the heated fluid rises due to the buoyancy effect, and then it turns right before reaching the top wall. Eventually, the fluid flows through the top gaps between blocks or blocks and wall, and leave the cavity alongside the external vertical wall. The observations are in excellent agreement with the findings in Ref. [50]. In addition, it shows that the streamlines smoothly cross through the coupling interfaces.Fig. 19. Two-dimensional natural convection in a side open cavity with porous media: local streamline distributions
Table 7 quantitatively compares the averaged Nusselt numbers on the heated wall obtained by the present coupling model and Ref. [50] for different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} . A maximum relative error of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.95~\%$$\end{document} at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa =10$$\end{document} is obtained.Table 7. Two-dimensional natural convection in a side open cavity with porous media: comparisons of Nusselt number for different ratios of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} PresentRef. [50]Relative error0.11.6491.6681.14%11.9021.8731.55%102.5152.4432.95%1002.7352.6672.55%
Additionally, the dimensionless mass flow rate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_{\textrm{in}}$$\end{document} , entering the cavity through the side opening ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x/L=1$$\end{document} ) is used to evaluate the accuracy of the coupling scheme, and it is computed as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} m_{\textrm{in}} = \left\{ \begin{aligned}&-\int _0^L \frac{u \,L}{\alpha _f} \, \textrm{d}y, ~~~~~~~\text {when}~u<0, \\&~~0, ~~~~~~~~~~~~~~~~~~~~~~~~~\text {when}~u\ge 0. \end{aligned} \right. \end{aligned}$$\end{document}Table 8 lists the mass flow rates obtained from our model and Ref. [50], together with the relative errors. It shows the mass flow rate increases with the increase of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} , proving the improvement of heat transfer inside the cavity. However, our values are smaller, with a relative discrepancy around \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8.0~\%$$\end{document} . The mismatch is induced by the different boundary condition on the top, right and bottom boundaries in thermal reservoir since code_saturne does not support uncoupled velocity components on the boundary as adopted in Ref. [50].Table 8. Two-dimensional natural convection in a side open cavity with porous media: comparisons of mass flow for for different ratios of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} PresentRef. [50]Relative error0.13.8024.1578.54%14.3594.7498.21%105.2375.7008.12%1005.4835.9608.00%
Two-dimensional normal plate velocity
The normal plate velocity case is computed to evaluate the numerical performance of the coupling model under the forced convection heat transfer condition. The coupled computational domain is displayed in Fig. 20 where blue and red regions are solved by LBM and FVM, with height 0.55H and 0.5H, and an overlap of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_o=0.05H$$\end{document} . Height H is discretised with 200 grid points, while only five points are used in the x direction. Periodic boundary conditions are applied to the left and right sides. The upper edge moves with a constant uniform horizontal rightward velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_x$$\end{document} whilst a uniform vertical upward flow with velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_y$$\end{document} is injected through the bottom wall and withdrawn from the top one. The constant high ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_h$$\end{document} ) and low ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c$$\end{document} ) temperatures are set at the top and bottom borders with dimensionless units \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _h=1.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _c=0.0$$\end{document} . The fluid is initialised at rest with an averaged temperature \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _0 = 0.5$$\end{document} . The fluid kinematic viscosity is set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu =0.01$$\end{document} and the Reynolds number is defined by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{Re}={u_y H / \nu } = 10$$\end{document} . Three simulations are carried out by varying the Prandtl number as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{Pr}= 0.1, 1, 10$$\end{document} . The problem admits analytical solution in the form
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \theta = {T-T_c \over T_h-T_c} = {{\textrm{exp}( \mathrm {Re \, Pr} ~y/L) -1} \over {\textrm{exp}(\mathrm {Re \, Pr})-1}}. \end{aligned}$$\end{document}Fig. 20. Normal plate velocity: sketch of geometrical setup
With reference to Fig. 21, the present results are in excellent agreement with the analytical solution.Fig. 21. Two-dimensional normal plate velocity: comparisons of temperature distribution between the coupling and analytical solutions for different Pr numbers
The relative errors between the numerical and analytical solutions calculated by Eq. (38) are reported in Table 9.Table 9. Two-dimensional normal plate velocity: percentage relative discrepancy between present results and analytical predictionsPr0.1110 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \, [\%]$$\end{document} 0.412.003.65
Two-dimensional thermal lid-driven cavity flow
The thermal lid-driven cavity flow is another popular benchmark case to assess the performance of a numerical tool. Figure 22 describes the computational domain of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} . The blue square domain of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(0.66L)^2$$\end{document} is solved by the LBM while the domain except the red part is simulated by the FVM. A uniform resolution (1/500) is enforced to the entire computational domain, with 10 grid points in the overlapping region, perpendicular to the coupling interfaces. The lid moves to the right with a uniform velocity ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U_{lid}$$\end{document} ) while the remaining boundaries are set to no-slip.Fig. 22. Two-dimensional thermal lid-driven cavity flow: sketch of the computational domain
In addition, the low ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c$$\end{document} ) and high ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_h$$\end{document} ) temperatures are enforced on the top and bottom borders with dimensionless units \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _c=0.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _h=1.0$$\end{document} , whilst left and right boundaries adopt the adiabatic condition. The time steps are set equal to 0.02 and 0.002 on the FVM and LBM sides. The problem is governed by the Prantl number with Pr \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=\nu /\alpha =0.71$$\end{document} , the Grashof number with Gr \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle = {g\beta (T_h-T_c) L^3 / \nu ^2} = 10^6$$\end{document} , the Reynolds number with Re \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle = {U_{lid} L / \nu }$$\end{document} and the Richardson number Ri \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle = { \mathrm {Gr / Re^2}}$$\end{document} . Three scenarios are investigated by varying the Richardson number as Ri \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$= 0.1,~1,~10$$\end{document} . The kinematic viscosity is set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu =4.0\times 10 ^{-5}$$\end{document} and the thermal diffusivity is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\displaystyle \alpha ={\lambda / \rho c_p}=5.6\times 10 ^{-5}$$\end{document} .
Figure 23 compares the streamline distributions between the coupling framework and pure FVM model using code_saturne at Ri = 10, 1 and 0.1, and the results from the pure FVM simulation are selected as the baseline. At Ri \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$= 10$$\end{document} , two counter-rotating vortices are generated on the upper and lower parts thanks to the interaction of the buoyancy effect because of the unstable temperature gradient and the shear force due to the moving lid. Moreover, the lower vortex exhibits a larger size compared to the upper vortex, as the influence of buoyancy dominates over shear stress at Ri \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10$$\end{document} . The hot and cold fluids converge in the central region, where the densely packed streamlines are observed. A secondary vortex forms in the lower-left corner due to the dominant buoyancy effect. As the Ri number decreases to 1, the shear effect intensifies and balances the buoyancy effect, producing two counter-rotating vortices of comparable strength. However, the buoyancy-induced vortex structure becomes distorted due to the opposing motion of the lid. Notably, the secondary vortex persists in the same location as observed in the Ri \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10$$\end{document} case. When Ri further decreases to 0.1, the lid-driven effect dominates the entire cavity, forming a primary vortex at the centre and three small secondary vortices near the corners. In addition, the streamlines exhibit smooth continuity across the coupling interfaces, and the vortex structures in the coupled model agree well with those obtained from the pure FVM simulation. This consistency validates the accuracy of velocity coupling scheme.Fig. 23. Two-dimensional thermal lid-driven cavity flow: comparisons of streamline distribution between the coupled and pure FVM models
Figure 24 presents the dimensionless temperature distributions with isotherms between the coupled and pure FVM models at Ri \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10$$\end{document} , 1 and 0.1. The isotherms are plotted at constant intervals of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \theta = 0.1$$\end{document} , denoted by adjacent solid lines. Excellent agreement is observed between the coupled model and pure FVM results, with isotherms demonstrating smooth continuity across coupling interfaces, demonstrating the correctness of temperature coupling scheme. At Ri \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10$$\end{document} and 1, enhanced thermal mixing in the central region produces steeper temperature gradients due to the interaction between hot and cold fluids. When Ri decreases to 0.1, the lid-driven flow dominates the entire cavity, resulting in cold fluid occupying most of domain. The isotherm distributions reveal particularly steep temperature gradients along the bottom heated wall, attributable to combined shear effects and fluid impingement phenomena.Fig. 24. Two-dimensional thermal lid-driven cavity flow: comparisons of temperature field with isotherms between the coupled and pure FVM models
Table 10 shows the values of Nusselt number computed along the heating wall, comparing results from the pure FVM, coupled model, and Refs [51, 52]. The pure FVM simulations demonstrate closer agreement with reference data. Using the values from Ref. [51] as baseline, our coupled model shows relative errors of 5.60 %, 6.1 % and 3.69 % at Ri \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10$$\end{document} , 1 and 0.1, respectively. This discrepancy primarily stems from the local difference between constant FVM solver density and computed variable density in the LBM boundaries.Table 10. Two-dimensional thermal lid-driven cavity flow: comparisons of Nusselt number on the heating wall at Ri \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10$$\end{document} , 1 and 0.1RiRef. [51]Ref. [52]Pure FVMPresent104.8604.8484.8624.58815.7505.7395.7585.4020.112.16112.13812.58012.610
Three-dimensional natural convection in a cavity
The proposed framework is further validated against the three-dimensional natural convection flow in a cubic cavity. Figure 25 shows the computational domain of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^3$$\end{document} and a cross-section along the symmetry plane located at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=L/2$$\end{document} . The blue region of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.66L^3$$\end{document} is assigned to the LBM solver, while the rest of the domain except the red part is solved by the FVM. The no-slip condition is enforced for all sides. Constant temperatures corresponding to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _h=1.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _c=0.0$$\end{document} are enforced on the planes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x=0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x=L$$\end{document} , respectively, while the other boundaries adopt the adiabatic condition. Two values of the Rayleigh number are investigated, i.e., Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} .Fig. 25. Three-dimensional natural convection: a computational domain partition and b a cross-section on symmetrical plane
At Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} , a uniform grid density is applied to the FVM and LBM domains with the resolution of 1/100. In the overlapping region, the distance perpendicular to coupling interface is set to 0.1L. The kinematic viscosity and thermal diffusivity are set equal to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.667\times 10^{-4}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.347\times 10^{-4}$$\end{document} , respectively, so the Prandtl number is equal to 0.71. The time step on the FVM and LBM sides are 0.05 and 0.01, respectively, so each FVM time step corresponds to 5 sub-iterations on the LBM side.
Figure 26 presents the dimensionless temperature distributions along the horizontal and vertical centrelines on the plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=0.5L$$\end{document} , comparing the results obtained from the pure FVM, LBM and reference data [53], where a third-order accurate finite difference approach was adopted. The pure FVM and pure LBM cases are solved using code_saturne and LUMA, respectively. Horizontally, five temperature distributions are displayed at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y/L=0.1,~ 0.3,~ 0.5,~ 0.7$$\end{document} and 0.9. Note that two monitoring lines ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y/L=0.1$$\end{document} and 0.9) lie completely within the FVM domain. The temperature profiles are very close to the results from pure FVM, LBM and Ref. [53], and the smooth transitions of temperature distribution through the coupling interfaces manifest the correctness of thermal coupling scheme.Fig. 26. Three-dimensional natural convection: dimensionless temperature distributions along the horizontal (left) and vertical (right) centrelines on the plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=0.5L$$\end{document} at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document}
Following the above horizontal and vertical lines on the plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=0.5L$$\end{document} , Fig. 27 compares the velocity distributions, which are normalised by the reference velocity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{\textrm{ref}}=\sqrt{g\beta h \Delta T}$$\end{document} , between the pure FVM, pure LBM, and Ref. [53]. Our present results agree well with those of the pure FVM, LBM and reference. In the overlapping region, the same velocity distributions can be observed in the coupled FVM and LBM regions.Fig. 27. Three-dimensional natural convection: dimensionless velocity distributions along the horizontal (left) and vertical (right) centrelines on the plane of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=0.5L$$\end{document} at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document}
For a higher Rayleigh number of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} , the grid resolution is refined to 1/200 on both FVM and LBM sides. In the overlapping region, the distance perpendicular to the coupling interfaces decreases to 0.05L. The time steps on the FVM and LBM sides are set to 0.01 and 0.001, respectively. The kinematic viscosity and thermal diffusivity are set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4.17\times 10^{-4}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5.86\times 10^{-4}$$\end{document} . In this case, the results from the pure FVM computed by code_saturne are selected as baseline and used for comparison with those from the coupling model. The variation of temperature across the horizontal and vertical centrelines on the plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=0.5L$$\end{document} are shown in Fig. 28. Notably, the temperature distributions are identical between the pure FVM and coupling model.Fig. 28. Three-dimensional natural convection: dimensionless temperature distributions along the horizontal (left) and vertical (right) centrelines on the plane of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=0.5L$$\end{document} at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document}
Furthermore, Fig. 29 presents a quantitative comparison of the dimensionless velocity distributions along the horizontal and vertical centrelines for the pure FVM, the coupled model, and data in Ref. [53]. Interestingly, our results closely match those of the pure FVM and the reference. For the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_x/u_{\textrm{ref}}$$\end{document} profile, a slight deviation is observed near the top boundary between the numerical models and the reference. Nevertheless, the results from both the pure FVM and the coupled model remain physically consistent, as the system’s symmetry naturally leads to symmetrical velocity distributions.Fig. 29. Three-dimensional natural convection: dimensionless velocity distributions along the horizontal (left) and vertical (right) centrelines on the plane of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z=0.5L$$\end{document} at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^6$$\end{document}
Across the entire coupling domain, Fig. 30 illustrates the isothermal surface distributions at Ra = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} with a maintained dimensionless temperature difference of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \theta = 0.1$$\end{document} between adjacent surfaces. In both cases, the isothermal surfaces exhibit smooth continuity through the overlapping region, though minor divergences near the coupling interfaces due to the weak compressibility effects inherent in the LBM.Fig. 30. Three-dimensional natural convection: isothermal surfaces at Ra = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^5$$\end{document} (left) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} (right)
Rayleigh-Bénard convection with a melting boundary
The advantages of the proposed coupled approach are demonstrated in complex melting processes, where the complementary strengths of LBM and FVM are effectively combined. The LBM excels at resolving local phase-change dynamics and associated thermal flows, while the FVM is more suitable for large-scale single-phase thermal flow simulations. To exemplify this concept, Fig. 31 illustrates the two-dimensional computational domain of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L\times L/2$$\end{document} , with the LBM (blue) and FVM (red) domains occupying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L\times 0.33L$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L\times 0.21L$$\end{document} , respectively. The overlapping domain is set to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L\times 0.04L$$\end{document} . Additionally, a consistent grid resolution of L/500 is employed in both computational domains to resolve the phase-change dynamics.Fig. 31. Rayleigh Bénard convection with a melting boundary: sketch of the problem setup
The computational domain is initialised with the fluid and solid phases occupying the bottom and top regions, respectively, each with a height of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_f=H_s=0.25L$$\end{document} . The flow field begins from a quiescent state, with a uniform initial mass density of 1.0 throughout the domain. The temperature field is initially set to a constant value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _0 = 0.0$$\end{document} . Constant high and low temperatures ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_h$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c$$\end{document} ) are prescribed over the bottom and top boundaries, respectively, with dimensionless units of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _h=1.0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _c=0.0$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta _c$$\end{document} serves as the threshold temperature for solid melting (triggered when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta >\theta _c$$\end{document} ). Periodic boundary conditions are applied to the left and right boundaries, while no-slip velocity condition is enforced at the top and bottom boundaries. In order to properly account for the phase-change dynamics, additional source terms are incorporated into the momentum and energy Eqs. (1–3) and they read as follows [54, 55]
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \partial _t \boldsymbol{u} + \left( \boldsymbol{u} \cdot \boldsymbol{\nabla } \right) \boldsymbol{u}= & -\dfrac{1}{\rho _0} \boldsymbol{\nabla } p + \nu \boldsymbol{\nabla }^2 \boldsymbol{u} + \varepsilon \boldsymbol{g} \beta (T-T_0) - {\chi \over \rho _0} \boldsymbol{u}, \nonumber \\ {\partial _t T} + \boldsymbol{u} \cdot \boldsymbol{\nabla } T= & \boldsymbol{\nabla } \cdot (\alpha \boldsymbol{\nabla } T) - \displaystyle {\Lambda \over c_p }\partial _t \varepsilon , \end{aligned}$$\end{document}where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi =1-\varepsilon ^2$$\end{document} is a penalisation factor. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document} represents the liquid fraction ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =1$$\end{document} for liquid, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =0$$\end{document} for solid), and the value is initially set to 1 and 0 in the regions of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$y \in [0:0.25L]$$\end{document} and [0.25L : 0.5L], respectively.. The thermal properties include latent heat \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document} and the heat capacity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_p=1$$\end{document} . In addition to the Rayleigh and Prandtl number, the dynamics of this system is also governed by the Stefan number, St = \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_p \Delta T / \Lambda $$\end{document} , where the temperature difference is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta T = T_h - T_c$$\end{document} . In this study the dimensionless parameters are specified as St \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=1$$\end{document} , Pr \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=1$$\end{document} and Ra spans four orders of magnitude (i.e., Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^7$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^8$$\end{document} ). The characteristic time scale is defined as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0=\sqrt{L/(2 g \beta \Delta T)}$$\end{document} .
Figure 32 presents the temperature and liquid fraction at times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t/t_0=20$$\end{document} and 48 by considering Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} . The temperature field reveals the formation of two Bénard convection cells, driven by buoyancy forces under gravity, which induce a wavy melting front characteristic of symmetric convection regimes. These results align with the results shown by De Rosis & Giustini [54], where the LBM was adopted to perform the numerical simulations.Fig. 32. Rayleigh Bénard convection with a melting boundary: temperature and liquid fraction at dimensionless time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t/t_0=20$$\end{document} and 48 for Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} . Light blue and grey colours correspond to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =1$$\end{document} (liquid) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =0$$\end{document} (solid)
As shown in Fig. 33, the temperature field exhibits unsteady behaviour at the higher Ra number equal to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^8$$\end{document} , resulting in a significantly more irregular melting front morphology. Despite minor local deviations, the isotherms maintain consistent distribution across the overlapping region, demonstrating once again the robustness of the present methodology under unsteady flow conditions.Fig. 33. Rayleigh Bénard convection with a melting boundary: temperature and liquid fraction at dimensionless time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t/t_0=20$$\end{document} and 48 for Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^8$$\end{document} . Light blue and grey colours correspond to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =1$$\end{document} (liquid) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon =0$$\end{document} (solid)
To quantitatively evaluate the melting dynamics, Fig. 35 presents the temporal evolution of the domain-averaged liquid fraction, r, which is defined as
\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} r(t) = \frac{1}{L \times L/2}\sum _{\boldsymbol{x}} \varepsilon (\boldsymbol{x},t). \end{aligned}$$\end{document}Our findings reveal a distinct transition from an initially smooth stage ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t/t_0<10$$\end{document} ) to a sharp change in the melting behaviour. Notably, the rate of phase change, as characterised by the slope of the liquid fraction curve, exhibits an inverse relationship with the Ra number, that is, the higher Ra values correspond to reduced melting rates. This phenomenon can be attributed to the enhanced convective heat transfer at higher Ra (Fig. 34). Since the temperature acts as the primary factor governing the melting process, the variation of average temperature ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\theta }$$\end{document} ) in the fluid phase, representing the degree of mixing between the hot and cold fluids, is shown in Fig. 35 for different Ra numbers. Similar to the phase-change behaviour, consistent trends are observed in the average temperature evolution. Higher Ra conditions enhance fluid mixing and lead to a reduction in the the average temperature, thereby suppressing the overall phase-change process.Fig. 34. Rayleigh Bénard convection with a melting boundary: time evolution of the liquid fraction at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^7$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^8$$\end{document} Fig. 35. Rayleigh Bénard convection with a melting boundary: time evolution of the averaged temperature in fluid phase at Ra \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=10^5$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^6$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^7$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^8$$\end{document}
Conclusions and future work
We have developed a robust coupling framework between the finite volume method and the lattice Boltzmann method for thermal flow simulations, employing a central-moments-based collision operator and advanced reconstruction techniques. The framework is implemented through the code_saturne and LUMA solvers, integrated using the Parallel Location Exchange (PLE) coupling library for efficient communication between FVM and LBM solvers.
Quantitative validation against seven benchmark problems, including both diffusive and convective regimes, complex geometries, and phase change, demonstrates the accuracy and versatility of the proposed approach. Notably, the use of central moments in both velocity and temperature fields reduces interface discontinuities and improves convergence: for example, in natural convection at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {Ra} = 10^6$$\end{document} , temperature residuals dropped by two orders of magnitude compared to BGK-based implementations. The reconstructed temperature and velocity profiles show excellent agreement with analytical and reference data across all cases.
Overall, the method constitutes a significant advance in the coupling of LBM and FVM, enabling accurate and efficient multiscale thermal flow modelling. Future work will focus on exploring high-order temporal interpolation scheme and extending the domain-decomposition-based framework to moving boundaries and phase-change-driven flows characteristic of realistic industrial applications.