International Journal of Heat and Mass Transfer xxx (2016) xxx–xxx

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer

journal homepage: www.elsevier.com/locate/ijhmt

A simple lattice Boltzmann model for conjugate heat transfer research

S. Chen a,b,⇑, Y.Y. Yan b, W. Gong b
a
Institute for Modelling and Simulation in Fluodynamics, Nanoscience and Industrial Mathematics ‘‘Gregorio Millán Barbany”, Universidad Carlos III de Madrid, Leganes 28911, Spain
b
Faculty of Engineering, The University of Nottingham, University Park, Nottingham NG7 2RD, UK

a r t i c l e i n f o a b s t r a c t

Article history: In this paper a lattice Boltzmann (LB) model is proposed for conjugated heat transfer research. Through
Received 29 July 2016 taking the most advantages of the standard LB method, the present model can remedy the shortcomings
Received in revised form 19 October 2016 of the available related LB models via a simple way and meanwhile a number of intrinsic advantages of
Accepted 31 October 2016
the standard LB method are preserved. It does not require any specific treatment dependent on interface
Available online xxxx
topology and independent from the choice of lattice model. Moreover, it can be used for unsteady prob-
lems with complicated and time dependent interfaces. The accuracy and reliability of the present model
Keywords:
are validated by three nontrivial benchmark tests. The good agreements between the present numerical
Lattice Boltzmann method
Conjugate heat transfer
prediction and available open data demonstrate the applicability of the present model for complicated
Heterogeneous media conjugated heat transfer problems. Finally, the present model could be extended to some other important
areas straightforwardly, such as fluid–solid phase change modeling.
Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction complication will be dramatically enhanced for a conjugate prob-

lem with a complex interface [3].
To meet a diversity of purposes of practical applications, almost During the past three decades, the lattice Boltzmann (LB)
all realistic systems are constituted by heterogeneous/composite method has attracted increasing attention due to its some intrinsic
media. Due to the balance between cost and performance, in the advantages, such as relatively easy treatment of complicated
manufacturing industry, the popularly used heterogeneous/ geometry, high parallel computing efficiency and capturing inter-
composite media are made up of several material layers with action between different phases/components at a mesoscopic level
different thermodynamic properties. Consequently, complicated [4]. Especially, as it is a particle-based numerical solver, the LB
interfaces between different material layers are commonly found method can guarantee, automatically, the continuity of a certain
in realistic systems. In the area of heat transfer research, these macroscopic quantity and of its flux across an arbitrary interface
interfaces are classified as conjugate problems where the within the investigated domain, if the macroscopic quantity and
so-called conjugate boundary condition is applicable [1]. Namely, its flux can be recovered from the zeroth- and first-order moment
on the interfaces, the temperature profile and heat flux both should of the corresponding pseudo-particle distribution function, respec-
be continuous, according to the local thermodynamic equilibrium tively, in the LB framework. This feature is a potential great advan-
assumption and the law of energy conservation. From the view- tage for conjugate heat transfer research as in the LB framework
point of scientific computing, the challenge of a conjugate problem one need not explicitly treat the topology of the interface where
may come from solving a strongly ‘‘coupled” interface domain the conjugated boundary condition should be strictly satisfied.
which is subject to Dirichlet- and Neumann-like boundary restric- Wang et al. [5] perhaps are the pioneers to utilize this upside. They
tion at the same time (but please bear in mind that a conjugate constructed a LB model for conjugate heat transfer simulation
boundary condition is neither a standard Dirichlet/Neumann without any explicit treatment on interfaces and observed that
boundary condition nor a combination of them). No doubt, how the LB approach could solve conjugate heat transfer problems
to accurately and efficiently solve a conjugate problem has been more efficiently than traditional numerical techniques, benefitting
an important topic in heat transfer research [2,3]. Especially, to from such intrinsic advantage of the LB scheme. Unfortunately,
the popularly used numerical techniques, their algorithmic meanwhile their model suffers from three obvious shortcomings:
(1) in their model any interface has to be located exactly on the
middle point between two lattice grids; (2) the model can work
well only in steady scenarios; and (3) the specific heat capacity
⇑ Corresponding author at: Faculty of Engineering, University of Nottingham, UK.
should be homogenous within the investigated domain even the
E-mail address: [email protected] (S. Chen).

http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120
0017-9310/Ó 2016 Elsevier Ltd. All rights reserved.

Please cite this article in press as: S. Chen et al., A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat Mass Transfer (2016),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120
2 S. Chen et al. / International Journal of Heat and Mass Transfer xxx (2016) xxx–xxx

domain is filled with heterogenous media. In order to remedy the first-order temporal accuracy due to the introduction of the source
latter two defects, the same research group [6] proposed another term while the temporal accuracy of the standard LB method is
LB model. Although their new model reached their expectation second-order). Moreover, their strategy will become very compli-
successfully, the major superiority of the LB method for conjugate cated if it is extended for conjugate fluid–solid problems as where
problem modeling was lost. Namely, in Ref.[6] the interfaces one has to calculate spatial gradients of enthalpy flux emerging in
between heterogenous media should be handled explicitly. As a the source term. Numerical accuracy and stability both will be
result, their discussions were limited for the cases restricted by hampered seriously. In addition, the intrinsic advantage of local
straight-interface geometry [6]. To model conjugate problems with computing of the standard LB method will be damaged by the spa-
arbitrary interfaces, Li et al. [7] and Le et al. [8] designed an tial gradients in their source term.
improved scheme, respectively. However, as interpolation and/or In this work, we propose a simple LB scheme for conjugate heat
extrapolation are required, these schemes, especially the former, transfer research. It does not need any specific treatment on arbi-
are complicated and lose the intrinsic advantage of local comput- trary interfaces and meanwhile can keep the intrinsic advantages
ing of the standard LB method. Moreover, in many scenarios the of the standard LB method, such as local computing. Therefore,
normal heat flux across an interface can be hardly evaluated, which the present scheme can remedy the shortcomings of the available
limits the applicable range of these schemes. An immersed- relevant models for conjugate heat transfer simulation. Moreover,
boundary-like LB scheme was also proposed by Hu et al. [9] to deal it can be extended to model solid–liquid phase change problems
with conjugate heat transfer across curved interfaces. In their work straightforwardly. The rest of this paper is organized as follow:
curved interfaces were approximated by zigzags and one had to in Section 2 we present this simple LB model together with a com-
assume that the ratio of thermal conductivities and the ratio of parison analysis to previous efforts and the corresponding numer-
thermal diffusivities should be identical [9]. The former treatment ical validation is conducted in Section 3, followed by a conclusion
will change the hydrodynamic characteristics of the investigated on this work. Although in the present work we take the single-
object and the latter one is hardly satisfied by realistic systems. relaxation-time (SRT) LB model as an example to show how to
In addition, how to choose a suitable delta function is a crucial establish a simple LB approach for conjugate heat transfer simula-
issue for the immersed boundary method, which will influence tion, its multiple-relaxation-time (MRT) counterpart is also easy to
numerical accuracy critically. Unfortunately, it is a challenge to be constructed, following the MRT template used in Ref. [8].
determine an appropriate multidimensional delta function [10].
In order to fully utilize the intrinsic advantage of the LB method
2. LB model for conjugated heat transdfer
to model conjugate heat transfer with complicated interfaces,
which is not achieved by the aforementioned efforts [5–9], Karani
2.1. Macroscopic governing equation for conjugate heat transfer
and Huber [11] proposed a different strategy. They pointed out, for
the first time, the root why the standard LB method can not work
For conjugate heat trnasfer modeling, the best choice is the con-
for unsteady conjugate heat transfer simulation is the spatial vari-
servation form of energy equation which reads [11]:
ation of the heat capacitance (the product of the density and speci-
fic heat capacity) across interfaces can not be modelled @ t qC p T þ ra qC p Tua ¼ ra kra T ð1Þ
appropriately by the standard LB scheme. Accordingly, the authors
designed a source term, which was added into the standard LB where q; ua and T are the density, velocity and temperature of
evolving equation, to correct the errors generated by the standard working media, respectively. In addition, k and C p denote the ther-
LB method. Different from the schemes proposed in Refs. [6–9], mal conductivity and constant pressure specific heat capacity.
their strategy avoids any specific treatment depending on the As analyzed in detail by Refs. [11–13], however, the recovered
interface topology. However, it was observed that their strategy macroscopic energy governing equation by the standard LB
lacks of mathematical rigor in differentiating the piecewise capac- method reads:
itance constant function [12]. As the key to extend the standard LB
scheme to conjugate heat transfer research is to capture the spatial @ t qC p T þ ra qC p Tua ¼ ra jra qC p T ð2Þ
variation of the heat capacitance across interfaces accurately,
where j ¼ k=ðqC p Þ is the thermal diffusivity of working media. Con-
Huang and Wu [13] claimed their LB model for phase change mod-
sequently, across any arbitrary interface within the investigated
eling could be adopted, after a minor modification, for conjugate
domain, in the framework of the standard LB method, only the fol-
heat transfer simulation, without any explicit treatment on inter-
lowing equalities can be guaranteed:
faces. However, in their model [13] only the spatial variation of
the specific heat capacity, rather than the heat capacitance, can Tþ ¼ T ð3Þ
be modelled, so their claim is tenable only when the density is
identical over the whole investigated domain. It is a too strong na ½jra qC p Tþ ¼ na ½jra qC p T ð4Þ
restriction to simulate conjugate problems in realistic systems.
Recently, a enthalpy-based LB model for unsteady conjugate heat where na is normal to the interface, and ½ þ and ½  indicate the
conduction between solid heterogeneous media was published parameters at each side of the interface.
[12]. A normalized sensible enthalpy, rather than temperature Compared with the conjugate boundary restriction [2,3]:
used in previous studies [5–9,11], was adopted to construct their
evolving equation [12]. In order to offset the additional terms Tþ ¼ T ð5Þ
brought by the introduction of the normalized sensible enthalpy,
a source term including temporal difference was added accord- na ½kra Tþ ¼ na ½kra T ð6Þ
ingly. For simulation of conjugated heat conduction across solid–
solid arbitrary interfaces, no special treatment on interfaces was there is an obvious difference between Eqs. (4) and (6). Only when
required in their scheme[12]. Such strategy can remedy the short- ½qC p þ ¼ ½qC p  or at the steady status, Eq. (4) can approximate to
coming caused by the determination of the average heat capaci- Eq. (6) exactly. Consequently, for the LB approach, the key to accu-
tance (in the source term) on interfaces, which has to be rately model conjugate problems is to recover the diffusion term in
addressed by the LB scheme developed in Ref. [11], but with a cost the energy equation exactly (i.e. the term at the right hand of Eq.
of degrading the numerical accuracy (the LB model in Ref. [12] is (1)).

Please cite this article in press as: S. Chen et al., A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat Mass Transfer (2016),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120
 S. Chen et al. / International Journal of Heat and Mass Transfer xxx (2016) xxx–xxx 3

2.2. Derivation of the present LB model for conjugate heat transfer With the aid of Eq. (13), we can replace ra r and ra r1 in the
modeling source term in Eq. (8) by

As revealed by Karani and Huber [11], the key to recover the dif- 1 X
fusion term in the energy equation Eq. (1) exactly is to treat the ra rðxa Þ ¼ xj rðxa þ eja DtÞeja ð14Þ
c2s Dt j
spatial variation of the heat capacitance appropriately. In the pre-
sent work, inspired by Ref.[12], a sensible-enthalpy-like quantity

h ¼ ðqC p Þ0 T is introduced, where ðqC p Þ0 is the reference heat and

capacitance. With the aid of h , Eq. (1) can be transformed to
(please refer to the Appendix for the detail): 1 1 X 1
ra ðxa Þ ¼ xj ðxa þ eja DtÞeja ð15Þ
k k 1 h
 r c2s Dt j r

@ t h þ ra h ua ¼ ra

r h  ra h ra  ua ra r ð7Þ
qC p a ðqC p Þ0 r r Consequently, the source term S can be written as
qC p
where r¼ is the ratio of heat capacitance. If the heat capaci-
ðqC p Þ0
1 1 X X 1
tance is homogenous (namely, ra r1 ¼ ra r ¼ 0), the last two terms S¼ rð1  Þ ½g j ðxa ; tÞ  g j ðxa ; tÞeja xj ðxa þ eja DtÞeja
ðeqÞ
c2s Dt 2s j j
r
in Eq. (7) will automatically vanish.
 1 h
 X
Eq. (7) is a standard advection–diffusion equation respect to h ,  ua xj rðxa þ eja DtÞeja ð16Þ
with a source term S ¼  ðqC p Þ ra h ra r  r ua ra r. It can be solved
k  1 h c2s Dt r j
0

by the following LB scheme [4,11]:

One can observe that in the present work the calculation of S is
1 a completely local operation. Away from interfaces, the output of
ðxa ; tÞ þ xj DtS
ðeqÞ
g j ðxa þ eja Dt; t þ DtÞ  g j ðxa ; tÞ ¼  ½g j ðxa ; tÞ  g j
s Eq. (16) becomes zero since there is no spatial variation of the heat
ð8Þ capacitance in a homogenous material layer, consistent to the
above statement.
where xj and Dt represents the weight coefficients and dimension- As the reference heat capacitance ðqC p Þ0 is a constant over the
less time step, and ej and g j ðxa ; tÞ denotes the discrete velocity whole domain, divided by ðqC p Þ0 , Eq. (7) can be transformed to:
direction and the corresponding pseudo-particle distribution func-
tion at the lattice grid xa . As discussed in previous research [8,14], k k 1 T
@ t T þ ra Tua ¼ ra ra T  ra T ra  ua ra r ð17Þ
for an advection–diffusion equation, a ‘‘simpler” lattice model with qC p ðqC p Þ0 r r
fewer discrete velocity directions can be used to improve computa-
tional efficiency, for example a D2Q5 lattice model for two- It is obvious that for homogeneous media Eq. (17) will reduce to
dimensional problem [8]. However, please bear in mind that the the popularly used macroscopic governing equation for tempera-
present scheme is independent from the choice of lattice model. ture field. The additional terms are the correction when the inves-
In Eq. (8) s is the dimensionless relaxation time and it satisfy tigated domain consists of heterogeneous media.
the relationship: In the present work, Eq. (7) is referred to as an ‘‘enthalpic” for-
mulation and Eq. (17) is its ‘‘temperature” counterpart. The advan-
k
j¼ ¼ c2s ðs  0:5ÞDt ð9Þ tage of the ‘‘enthalpic” form Eq. (7) is it can be extended
qC p straightforwardly to model solid–liquid phase change problems
as in such research area the enthalpy-base method has received
where cs is the speed of sound in the LB framework [4,14].
increasing attention due to its simplicity [9,16]. It will be discussed
The equilibrium distribution in Eq. (8) reads [14]
in detail in our future work as it is beyond the scope of the present


 eja ua work.
¼ xj h
ðeqÞ
gj 1þ 2 ð10Þ
cs Compared with previous efforts [6–9], the present model does
not need special treatments on the topology of the interfaces
and the local temperature can be obtained through the zeroth-order between heterogeneous media. Compared with Ref. [5], the pre-
moment of g j ðxa ; tÞ: sent model can model unsteady conjugate heat transfer with com-
plicated interfaces via a simple way. Compared with Karani and
1 X
T¼ g ð11Þ Huber’s model [11], the present work avoids the puzzle of deter-
ðqC p Þ0 j j
mining the average heat capacitance on interfaces. Compared with
Ref. [13], the present work breaks the limitation of identical den-
In succession, how to evaluate the source term S appropriately
sity over the whole investigated domain. Compared with the
becomes a critical step to guarantee the advantage of local com-
model designed in Ref. [12], the present model keeps the intrinsic
puting of the standard LB method. According to Eq. (13) in Ref.
advantage of local computing of the standard LB method, besides
[11], one can get
numerical stability and accuracy. Therefore, the present model
k k realizes to model conjugate problems through taking the most
 ra h ¼ r r h
ðqC p Þ0 qC p a intrinsic advantages of the LB method. Meanwhile, the present


1 X model possesses a number of attracting advantages, such as its
¼r 1
ðeqÞ
½g ðxa ; tÞ  g j ðxa ; tÞeja ð12Þ capablity to model conjugate heat transfer in media with evolving
2s j j
microstructures where interfaces and thermophysical properties
In addition, in the LB framework there is an approximation [4]: change over time (please see the Appendix for the detail).
Finally, if there are internal heat resources in the investigated
1 X domain, another source term that represents the internal heat
ra /ðxa Þ ¼ xj /ðxa þ eja DtÞeja þ OðDt2 Þ ð13Þ
c2s Dt j
resources should be added into Eq. (8) accordlingly, as shown in
Refs. [11,14]. It will be illustrated in detail through the third vali-
where /ðxa Þ is a certain scalar quantity at the lattice grid xa . dation benchmark case below.

Please cite this article in press as: S. Chen et al., A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat Mass Transfer (2016),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120
4 S. Chen et al. / International Journal of Heat and Mass Transfer xxx (2016) xxx–xxx

3. Numerical validation

In order to validate the present model, three non-trivial bench-

mark tests are adopted. The first one is transient heat conduction
in three-layered stratified media [12] which can validate the
present model’s accuracy for unsteady conjugate problems. The
second one is conjugate mixed convection heat transfer in a
lid-driven enclosure with thick bottom wall [18], which can check
the capability of the present model to conserve conductive and
advective heat flux simultaneously on interfaces as in which ther-
mal fluid flow is not always parallel to the interface. The last
benchmark test is magnetohydrodynamic (MHD) mixed convec-
tion with Joule heating in a lid-driven cavity which contains a heat
conducting horizontal circular cylinder [19]. It can demonstrate
the reliability of the present model for curved interfaces embeded
in complicated fluid and heat flow. The D2Q9 lattice is used for
solving flow field and the D2Q5 lattice is employed to compute Fig. 2. Temperature profiles along the vertical central line at various instants: solid
energy field, similar with our previous work [14]. The LB scheme lines-present results, dots-solutions in Ref. [12].
for fluid flow simulation and boundary treatment used in our pre-
vious work [14,15] is adopted here again. In our simulation, ðqC p Þ0
is set as an average value of various heat capacitance in the inves-
tigated domain, which is helpful to improve numerical stability.

3.1. Transient heat conduction in three-layered stratified media

The transient heat conduction in three-layered stratified media

investigated by Ref. [12] is a typical issue in heat transfer research.
The schematic configuration is depicted by Fig. 1. At time t ¼ 0, the
temperature over the whole domain is identical, as T c ¼ 0. Since
t > 0, the temperature on the downside is elevated and fixed at
T h ¼ 1:0. The vertical sides of the investigated domain both are adi-
abatic. Firstly, the thermophysical properties of each layer are set
as: k1 ¼ k3 ¼ 1; k2 ¼ 0:1; ðqC p Þ1 ¼ ðqC p Þ3 ¼ 1:0 and ðqC p Þ2 ¼ 2:0.
The same as Ref. [12], a grid resolution 100  300 is employed.
Fig. 2 illustrates the temperature profiles along the vertical cen-
tral line at different instants until the steady status obtained by the
present model, compared with the solutions given in Ref. [12] Fig. 3. Temperature profiles along the vertical central line at various instants: dots-
(please refer to Fig. 3 therein). The present results agree well with present results, dashed lines-analytic solutions in Ref. [17].
the data in Ref. [12], which demonstrates the present model’s
applicability for unsteady conjugate problems.
In order to check the applicability of the present model for high
heat capacitance ratio, the thermophysical properties proposed in
Ref. [17] are also adopted: k1 ¼ k3 ¼ 1; k2 ¼ 0:1; ðqC p Þ1 ¼
ðqC p Þ3 ¼ 1:0 and ðqC p Þ2 ¼ 0:1=3. The analytic solutions for this
case were available in Ref. [17]. Fig. 3 plots the comparison
between the present numerical results and the analytic solutions
in Ref. [17]. The agreement demonstrates the present model can

Fig. 4. Convergence behavior of the present model for transient heat conduction in
three-layered stratified media.

work well for high heat capacitance ratio. Fig. 4 depicts the conver-
gence behavior of the present model, where the normalized error E
is defined as
RX jT  T  j
E¼ ð18Þ
Fig. 1. Schematic configuration of three-layered stratified media. RX jT  j

Please cite this article in press as: S. Chen et al., A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat Mass Transfer (2016),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120
 S. Chen et al. / International Journal of Heat and Mass Transfer xxx (2016) xxx–xxx 5

Tc =0 , u=1 , v= 0 3.2. Conjugate mixed convection heat transfer in a lid-driven enclosure

with thick bottom wall

Fig. 5 shows the schematic configuration of lid-driven enclosure

T’
x =0 T’
x =0
with thick bottom wall investigated in Ref. [18]. The square cavity

H g is H  H filled with fluid. The vertical walls both are adiabatic. The
solid bottom wall with a thickness h is heated by a fixed high tem-
u= 0 u= 0
perature T h ¼ 1 while the moving top lid is cooled by a constant
temperature T c ¼ 0. The gravity g is downward. The boundary con-
v= 0 v= 0 ditions are illustrated by Fig. 5. The Richardson number Ri, which
describes the ratio between buoyant convection term and force
y
h convection term, is unity. In the present work h=H ¼ 0:1 and the
ratio of thermal conductivity between fluid and solid thick bottom
x H wall is k ¼ kf =ks ¼ 0:1; 1 and 10. The detailed description, such as
Th=1, u= 0 , v= 0 the governing equations, please refer to Ref. [18]. A grid resolution
100  110 is adopted.
Fig. 5. Schematic configuration of lid-driven enclosure with thick bottom wall.
Figs. 6–8 plots the streamlines and isotherms at various k. With
a low k (e.g. k ¼ 0:1), temperature gradient in fluid is large due to
and Dx denotes the grid interval. In Eq. (18), T  represents the ana- its low thermal conductivity. Against k increasing, fluid becomes
lytic solutions given by Ref. [17] and the subscript X means the colder and temperature is distributed inside the solid bottom wall
whole investigated domain. The slope of the error line is about owing to the insulation-like effect of the thick bottom wall. The
1:92, which indicates the accuracy of the present model approaches fluid at the top-right corner of the cavity is always the coldest over
to second order. the whole domain. The observation is in good agreement with that
reported by Ref. [18] (please refer to Figs. 3–5 therein).

0.95
0.9
0.85
0.8
0.75
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05

(a) (b)
Fig. 6. (a) streamlines and (b) isotherms for h=H ¼ 0:1; k ¼ 0:1 and Ri ¼ 1.

0.95
0.9
0.85
0.8
0.75
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05

(a) (b)
Fig. 7. (a) streamlines and (b) isotherms for h=H ¼ 0:1; k ¼ 1 and Ri ¼ 1.

Please cite this article in press as: S. Chen et al., A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat Mass Transfer (2016),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120
6 S. Chen et al. / International Journal of Heat and Mass Transfer xxx (2016) xxx–xxx

0.95
0.9
0.85
0.8
0.75
0.7
0.65
0.6
0.55
0.5
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05

(a) (b)
Fig. 8. (a) streamlines and (b) isotherms for h=H ¼ 0:1; k ¼ 10 and Ri ¼ 1.

The profile of y-component velocity along the horizontal central

line of the square cavity is depicted by Fig. 9. One can observe that
a higher k will prevent the fluid motion. It agrees with the conclu-
sion in Ref. [18] that flow strength was decreased against k
increasing.
Table 1 lists the comparison of the average Nusselt number Nu
on the interface between fluid and solid bottom wall, which is cal-
RH
culated by Nu ¼  H1 0 @T @y
dx. There are only slight differences
between the present numerical prediction and the published data
[18].
Fig. 10 illustrates the profiles of temperature and x-component
velocity along the vertical central line of the investigated domain,
with a comparison to the data obtained by the commercial CFD
software FLUENTÒ. The comparison shows the present model can
capture the fine structures of flow and temperature field
accurately.

Fig. 9. y-component velocity along the horizontal central line of lid-driven cavity.
3.3. MHD mixed convection with Joule heating in a lid-driven cavity
which contains a heat conducting horizontal circular cylinder
Table 1
Average Nusselt number on the interface with various k. Finally, the MHD mixed convection in a lid-driven cavity inves-
tigated in Ref. [19] is adopted to validate the present model’s accu-
k ¼ 0:1 k¼1 k ¼ 10
racy for modeling conjugate heat transfer with complicated
Ref. [18] 2.899 2.468 1.468 interfaces. As shown by Fig. 11, a heat conducting horizontal circu-
Present results 2.9696 2.5303 1.4635
lar cylinder, whose diameter is d, is located in the center of the cav-

Fig. 10. Profiles of (a) temperature and (b) x-component velocity along the vertical central line of the investigated domain for h=H ¼ 0:1; k ¼ 1 and Ri ¼ 1: dots-CFD, lines-
present results.

Please cite this article in press as: S. Chen et al., A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat Mass Transfer (2016),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120
 S. Chen et al. / International Journal of Heat and Mass Transfer xxx (2016) xxx–xxx 7

reflect the effect of Joule heating an additional term J e v 2 should

be added into the right hand of Eq. (1), where J e is the Joule heating
parameter and v is the y-component velocity. Therefore, a corre-
sponding additional source term xj DtJ e v 2 =r should be added into
the LB evolving equation Eq. (8). The dimensionless parameters
used in the present simulation are the Joule heating parameter
J e ¼ 1, the Richardson number Ri ¼ 1, the Hartmann number
Ha ¼ 10 and the Reynolds number Re ¼ 100. The ratio of thermal
conductivity between solid cylinder and MHD fluid K ¼ ks =kf varies
from 1 to 10. The diameter of the cylinder varies in the range
d=H ¼ 0:2 and 0:6. The detailed description about the investigated
domain and macroscopic governing equations please refer to Ref.
[19]. A grid resolution 100  100 is adopted.
Figs. 12 and 13 plot the streamlines and isotherms when the
size of the cylinder is changing while the solid–fluid thermal con-
ductivity ratio is fixed at K ¼ 5. The influence of cylinder size on
Fig. 11. Schematic configuration of MHD mixed convection with joule heating in
flow and temperature field is significant. With a small cylinder
lid-driven cavity containing a heat conducting horizontal circular cylinder.
(d=H ¼ 0:2), the fluid flow is featured by a pair of counter rotating
rolls. While for a large cylinder (d=H ¼ 0:6), the fluid flow is dom-
inated by a clockwise vortex generated by the movement of the left
ity. The horizontal walls of the cavity are adiabatic. The tempera- vertical wall, together with twin counterclockwise small vortices
ture on the right vertical wall is T h ¼ 1 while that on the left lateral caused by the buoyant force. The dramatic change of the flow field
wall is T c ¼ 0. The left lateral wall is moving with velocity v ¼ 1 results from that a large cylinder can reduce the available space for
while the other walls are stationary. The gravity g goes down- the buoyancy-induced recirculation. The convective distortion of
wards. In addition, the working MHD fluid is subject to uniform the isotherms near the right top corner of the enclosure decreases
external magnetic field B0 . Accordingly, the Joule heating should against d increasing. The observation agrees well with that
be considered, which acts like an internal heat source. In order to reported by Ref. [19].

0.05

0.15 0.95
0.85

0.25 0.75

0.35 0.65
0.45
0.55

(a) (b)
Fig. 12. (a) streamlines and (b) isotherms for d=H ¼ 0:2 and K ¼ 5.

0.05

0.15
0.25

0.45
0.55

0.35

0.65 0.85
0.95
0.75

(a) (b)
Fig. 13. (a) streamlines and (b) isotherms for d=H ¼ 0:6 and K ¼ 5.

Please cite this article in press as: S. Chen et al., A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat Mass Transfer (2016),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120
8 S. Chen et al. / International Journal of Heat and Mass Transfer xxx (2016) xxx–xxx

0.05 0.05
0.95
0.95
0.15 0.15
0.85
0.85

0.25

0.75

0.75
0.25

0.35 0.35 0.65

0.45 0.65
0.55
0.55 0.45

(a) (b)
Fig. 14. Isotherms of (a) K ¼ 10 and (b) K ¼ 10 at d=H ¼ 0:2.

the present model. Moreover, it can be used for unsteady problems

Table 2 with complicated and time dependent interfaces and changeable
Average Nusselt number on the left wall with various d and K ¼ 5. thermophysical properties. The accuracy and reliability of the pre-
d=H ¼ 0:2 d=H ¼ 0:6 sent model are validated by three nontrivial benchmark tests. The
Ref. [19] 1.022651 1.547871 good agreements between the the present numerical prediction
Present results 1.082642 1.620414 and available open data demonstrate the applicability of the pre-
sent model for complicated conjugated heat transfer problems.
Especially, the present model could be extended to some other
important areas, such as fluid–solid phase change modeling.
Table 3
Average Nusselt number on the left wall with various K and d=H ¼ 0:2.
Finally, although in the present study we only take a single-
relaxation-time LB model as an example to show how to address
K¼1 K ¼ 10 the variation of thermophysical properties of working media, the
Ref. [19] 0.981962 1.033299 extension to its multiple-relaxation-time counterpart is straight-
Present results 1.011489 1.105778 forward [8]. It will be considered in our future work.

Acknowledgments

Fig. 14 illustrates the effect of K on the distribution of isotherms. This work has received funding from the Universidad Carlos III
The corresponding streamlines are not plotted as the impact of K de Madrid, the European Unions Seventh Framework Programme
on flow field is relatively slight and they are very similar with for research, technological development and demonstration under
Fig. 12(a). According to this figure, one can observe that the ratio grant agreement No. 600371, el Ministerio de Economa y Compet-
of thermal conductivity between solid cylinder and fluid critically itividad (COFUND2014-51509), el Ministerio de Educacin, cultura y
affects the distribution of isotherms. With a high K (e.g. K ¼ 10), Deporte (CEI-15-17) and Banco Santander. S. Chen would also
the isotherms in the vicinity of the solid cylinder will be distorted acknowledge the support from the British Newton Alumni Fellow-
significantly. The above phenomena were also reported by Ref. ship Scheme.
[19].
Tables 2 and 3 list the average Nusselt number Nu on the left Appendix A. Derivation of Eq. (7)
wall at different configurations, compared with the open data in
Ref. [19]. It can be found that Nu increases with d and K, which is 
Replacing T in Eq. (1) by h , one can obtain the following
consistent with that reported by Ref. [19]. equation:

h
@ t rh þ ra rh ua ¼ ra kra
 
ðA:1Þ
4. Conclusion ðqC p Þ0
The same as all previous research [5–9,11–13], at current stage
Conjugated heat transfer is a popular problem in industry. For
we only consider the scenario where the heat capacitance varies
traditional numerical tools, how to solve conjugated problems
spatially but keep constant over time. Consequently, the left hand
with complicated geometry is still a challenge. Recently, some
side of Eq. (A.1) can be written as:
scholars have conducted a number of efforts to develop LB-based
approaches to address this challenge. Unfortunately, the shortcom-
ings of these LB-based approaches are obvious, too. In the present
work, we firstly analyze the advantages of the LB method for con-
jugate problems and reveal the root which restrict the standard LB
method to model conjugate heat transfer. Based on the analyses,
we propose a new LB model which can remedy the shortcomings With the aid of Eq. (A.2), Eq. (A.1) can be transformed to:

of the available relevant LB models via a simple way. Some impor- h
rð@ t h þ ra h ua Þ ¼ ra kra  h ua ra r

ðA:3Þ
tant advantages of the standard LB method are preserved. Any ðqC p Þ0
specific treatment dependent on interface topology is avoided by

Please cite this article in press as: S. Chen et al., A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat Mass Transfer (2016),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120
 S. Chen et al. / International Journal of Heat and Mass Transfer xxx (2016) xxx–xxx 9

Dividing Eq. (A.3) by r, there is: [5] J. Wang, M. Wang, Z. Li, A lattice Boltzmann algorithm for fluidCsolid conjugate
heat transfer, Int. J. Therm. Sci. 46 (2007) 228–234.
 
1 h h [6] F. Meng, M. Wang, Z. Li, Lattice Boltzmann simulations of conjugate heat
 ua ra r
 
@ t h þ ra h ua ¼ ra kra ðA:4Þ
r ðqC p Þ0 r
transfer in high-frequency oscillating flows, Int. J. Heat Fluid Flow 29 (2008)
1203–1210.
[7] L. Li, C. Chen, R. Mei, J.F. Klausner, Conjugate heat and mass transfer in the
As lattice Boltzmann equation method, Phys. Rev. E 89 (2014) 043308/1–043308/
 21.
1 h k k 1
ra kra ¼ ra r h  ra h ra ðA:5Þ [8] G. Le, O. Oulaid, J.F. Zhang, Counter-extrapolation method for conjugate
r ðqC p Þ0 qC p a ðqC p Þ0 r interfaces in computational heat and mass transfer, Phys. Rev. E 91 (2015)
033306/1–033306/11.
With the aid of Eq. (A.5), Eq. (A.4) can be written as [9] Y. Hu, D. Li, S. Shu, X. Niu, Full Eulerian lattice Boltzmann model for conjugate
heat transfer, Phys. Rev. E 92 (2015) 063305/1–063305/12.

k k 1 h [10] B. Yang, S. Chen, C. Cao, Z. Liu, C. Zheng, Lattice Boltzmann simulation of two
 
@ t h þ ra h ua ¼ ra ra h  ra h ra  ua ra r
qC p ðqC p Þ0 r r cold particles settling in Newtonian fluid with thermal convection, Int. J. Heat
Mass Transfer 93 (2016) 477–490.
ðA:6Þ [11] H. Karani, C. Huber, Lattice Boltzmann formulation for conjugate heat transfer
in heterogeneous media, Phys. Rev. E 91 (2015) 023304/1–023304/10.
namely Eq. (7). [12] H. Rihab, N. Moudhaffar, B.N. Sassi, P. Patrick, Enthalpic lattice Boltzmann
formulation for unsteady heat conduction in heterogeneous media, Int. J. Heat
If the heat capacitance varies over time, namely the last term of
Mass Transfer 100 (2016) 728–736.
Eq. (A.2) h @ t r – 0, it is also can be treated straightforwardly by

[13] R. Huang, H. Wu, Total enthalpy-based lattice Boltzmann method with
the present scheme. What should be done is just to add h @ t r into
 adaptive mesh refinement for solid-liquid phase change, J. Comput. Phys.
315 (2016) 65–83.
the source term S in Eq. (8). For the LB method, there is no obvious [14] S. Chen, Z. Liu, C. Zhang, Z. He, Z. Tian, B. Shi, C. Zheng, A novel coupled lattice
extra computational cost to evaluate h @ t r if the forward Euler

Boltzmann model for low Mach number combustion simulation, Appl. Math.
time discretization scheme is adopted. Comput. 193 (2007) 266–284.
[15] S. Chen, K. Luo, C. Zheng, A simple enthalpy-based lattice Boltzmann scheme
for complicated thermal systems, J. Comput. Phys. 231 (2012) 8278–8294.
References [16] D. Chatterjee, S. Chakraborty, An enthalpy-based lattice Boltzmann model for
diffusion dominated solid–liquid phase transformation, Phys. Lett. A 341
[1] T.L. Perelman, On conjugated problems of heat transfer, Int. J. Heat Mass (2005) 320–330.
Transfer 3 (1961) 293–303. [17] Y. Sun, I.S. Wichman, On transient heat conduction in a one-dimensional
[2] X. Chen, P. Han, A note on the solution of conjugate heat transfer problems composite slab, Int. J. Heat Mass Transfer 47 (2004) 1555–1559.
using SIMPLE-like algorithms, Int. J. Heat Fluid Flow 21 (2000) 463–467. [18] H.F. Oztop, C. Sun, B. Yu, Conjugate-mixed convection heat transfer in a lid-
[3] N. Sato, S. Takeuchi, T. Kajishima, M. Inagaki, N. Horinouchi, A consistent direct driven enclosure with thick bottom wall, Int. Commun. Heat Mass Transfer 35
discretization scheme on Cartesian grids for convective and conjugate heat (2008) 779–785.
transfer, J. Comput. Phys. 321 (2016) 76–104. [19] M.M. Rahman, M.A.H. Mamun, R. Saidur, S. Nagata, Effect of a heat conducting
[4] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, horizontal circular cylinder on MHD mixed convection in a lid-driven cavity
Oxford University Press, Oxford, 2001. along with joule heating, Int. J. Mech. Mater. Eng. 4 (2009) 256–265.

Please cite this article in press as: S. Chen et al., A simple lattice Boltzmann model for conjugate heat transfer research, Int. J. Heat Mass Transfer (2016),
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.10.120