The temperature decomposition method for periodic thermal flows with conjugate heat transfer

International Journal of Heat and Mass Transfer 150 (2020) 119288

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/hmt

The temperature decomposition method for periodic thermal flows with conjugate heat transfer Mayssaa Jbeili a, Junfeng Zhang a,∗ Bharti School of Engineering, Laurentian University, 935 Ramsey Lake Road, Sudbury, Ontario, P3E 2C6, Canada

a r t i c l e

i n f o

Article history: Received 25 September 2019 Revised 4 December 2019 Accepted 28 December 2019

Keywords: Heat transfer Periodic flow Temperature decomposition method Conjugate interface Computational fluid dynamics

a b s t r a c t Periodic structures are often found in many industrial applications. Simulations for such systems can be performed over a one-module domain when neglecting variation in fluid properties, and the results can be applied to a multiple-module flow passage. This approach is widely found in previous studies; however, it is limited to systems with relatively simple boundary situations: either the temperature or heat flux can be specified over the wall surfaces. Recently Li and Zhang (J. Heat Transf. 140:112002, 2018; Numer. Heat Tr. B Fund. 74:559, 2018) have proposed a temperature decomposition method, which makes the one-unit simulations possible for the above situations. However, the conjugate heat transfer, for example, between the flowing fluid and the structure materials (solid structure), has not been considered. In this study, we extend the temperature decomposition method to periodic thermal flows with conjugate heat transfer between the fluid and structure. The physical temperature is split into two parts, namely the transient part and the equilibrium part. The governing equations and conjugate relations at the interface for each temperature parts are established, and the corresponding inlet/outlet boundary conditions are proposed. It can be mathematically shown that the temperature obtained from this method satisfies the original energy equation in the fluid and solid domains, the thermal boundary conditions on the system walls, and the conjugate relations at the solid-fluid interfaces. Physical interpretations and numerical implementation steps are also provided. Several simulations are conducted to verify and demonstrate the validity, accuracy, and potential usefulness of our method, and the results and discussions clearly show that the TDM method correctly reflects the underlying relationship in temperature, which otherwise is not available, for conjugate periodic thermal systems. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction The advantage of periodic structures in heat transfer systems lies in their ability to enhance heat transfer with the increase of contact area between the fluid and the wall surface. These periodic structures are encountered in many industrial applications, such as heat exchangers and heat sinks [1–3], and can also be found in recent microfluidic and nanofluid applications [4,5]. When the fluid property changes in a periodic structure can be neglected, identical fluid flow and similar temperature distribution can be observed in periodic units after a certain distance from the flow entrance. Flows in this situation are called fully developed periodic thermal flows [6–8]. For such systems, the flow and temperature fields can be calculated by considering only one periodic unit with appropriate inlet/outlet boundary conditions, whereas the results can be applied to a flow passage of multiple unit length [6,7,9–

∗

Corresponding author. E-mail address: [email protected] (J. Zhang).

https://doi.org/10.1016/j.ijheatmasstransfer.2019.119288 0017-9310/© 2020 Elsevier Ltd. All rights reserved.

12]. Clearly this approach is favorable for using less computational resources; however it is limited to systems with a uniform temperature or with explicitly specified heat flux over all surfaces. This requirement cannot be met in many heat transfer systems, such as microchannel heat sinks where fluid flows between surfaces with different wall temperatures [13]. To address this limit, a temperature decomposition method (TDM) was recently developed by Li and Zhang [14,15] to simulate periodic thermal flows with general wall conditions, such as non-uniform wall temperatures, mixed temperature and heat flux conditions, and convective boundary conditions. In this method, the physical temperature is split into two parts: the transient part that reflects the temperature variation among periodic units, and the equilibrium part representing the thermal equilibrium state far downstream from the entrance. Similar to previous methods, TDM carries the one-unit computational advantage and the calculation results are readily applicable to multiple units along the periodic flow passage. These one-unit methods for periodic thermal flows [6,10,15,16] have considered the temperature fields only in the

2

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

fluid domain and thus require explicit information of the thermal conditions over solid surfaces. In practical heat transfer devices, such information is usually not available. For instance, periodic thermal flows through wavy and corrugated pipes have been studied extensively in the literature [6,7,9,10,17]. However, the thermal condition on the inner surface, which is required by existing one-unit methods, is complex and difficult to be described even when the outer surface condition is relatively simple [18]. For a less ambiguous example, we can consider a pin-fin heat sink for microelectronic devices [19]. The heat generated by microelectronic elements (chips, batteries, etc.) is first passed to the heat sink from the base surface, and then it is removed by the cooling fluid flowing through the void space among the fins. The exact thermal condition on the fluid-solid contact area (fin surfaces) varies from location to location, and it cannot be described in advance. Actually all system parameters can affect the thermal condition on the fin surface, including the geometry and configuration of fins, the thermal properties of the heat sink material, the heat generation amount and distribution at the base, and the fluid flow condition (fluid properties and flow field). To simulate such systems with conjugate heat transfer, the energy equation needs to be solved simultaneously in both solid and fluid domains, with appropriate conjugate conditions considered at the solid-fluid interface. Kelkar and Choudhury [20] calculated the flow and thermal fields with natural convection due to periodic heating blocks in a vertical adiabatic channel, and the temperature and heat flux continuity was considered at the solid-fluid interface. With the heat generation rate in each block available, the temperature field can be readily decomposed into a linear variation part to account for the heat addition from blocks and a periodic part identical in all modules. However, the adiabatic boundary condition is often invalid and the method in Ref. [20] cannot be utilized in general periodic thermal flow systems with conjugate interfaces. In this paper, we will extend the work by Li and Zhang [15] to periodic thermal systems with the solid-fluid conjugate heat transfer considered. Following the treatments in Refs. [14,15], the physical temperature is split into the transient and equilibrium parts in both solid and fluid domains. In addition to the governing energy equations for the two temperature parts, respective thermal conditions are developed for system boundaries, the inlet and outlet, and the conjugate interface. The transient and equilibrium temperature parts can then be solved independently for one periodic unit of the structure, and the simulation results can then be used to construct the temperature fields along the flow passage with the bulk inlet temperature provided. It can be shown mathematically that the temperature obtained from our approach satisfies the governing equations for thermal energy and the system boundary/interface conditions, and therefore, is indeed the true temperature solution of the system. We also discuss the physical meanings of the transient and equilibrium parts. At last, several demonstration examples are presented to illustrate potential applications of our method in periodic heat transfer systems. The buoyancy effect in natural convection systems could also be considered by following similar treatments in Refs. [12,20,21]. 2. Theory and method description The two-dimensional (2D) periodic system shown in Fig. 1 is used as a general example to describe the fully developed periodic thermal flows with conjugate interface. A pressure difference between inlet and outlet is necessary to drive the fluid flow. For the fluid domain f , both the flow and thermal fields should be solved; on the other hand, for the solid part s , only temperature is involved. When solving the temperature fields in the fluid and solid domains, in addition to the explicit thermal conditions applied on the system boundaries  f and  s , appropriate relations

of temperature and its gradient at the conjugate interface  i also need to be satisfied. Next we describe the mathematical formulations and numerical strategies for solving the flow and thermal distributions in this system using only one periodic unit. 2.1. Periodic fluid flow The flow field in the fluid domain f is governed by the following incompressible continuity (Eq. (1)) and momentum (Eq. (2)) equations:

∇ · u = 0,

(1)

∂u ∇P + ∇ · (uu ) = − + ν∇ 2 u, ∂t ρ

(2)

where u is the flow velocity, P is the pressure, ρ is the fluid density, ν is the kinematic viscosity, and t is time. Here we assume all fluid properties constant and therefore a fully developed periodic flow can be established after some distance from the entrance. At this state, the following relationships can be written for the velocity u and pressure P at locations of the same relative position in different periodic modules [6,10,11,16]:

u(x ± mL, y ) = u(x, y ),

(3)

P (x ± mL, y ) = P (x, y ) ∓ mPL ;

(4)

where L is the streamwise length of the periodic unit, PL is the pressure drop over each periodic unit, and m is a natural number. Assuming the flow and pressure field is given by u(x, y) and P(x, y) in one particular module (0 ≤ x ≤ L), it is easy to show that u(x ± mL, y) and P(x ± mL, y) from Eqs. (3) and (4) will satisfy the continuity and momentum equations Eqs. (1) and (2) exactly in the module ± mL ≤ x ≤ (1 ± m)L; since the identical velocity and shifted-by-a-constant pressure profiles cause no change to any terms in the governing differential equations. The no-slip boundary condition, if satisfied in the original unit (0 ≤ x ≤ L), will be satisfied as well in the module ± mL ≤ x ≤ (1 ± m)L. Therefore, once available in one periodic module, the velocity and pressure fields in other modules can be readily obtained via Eqs. (3) and (4). To solve the flow field in one periodic unit as shown in Fig. 1b, additional inlet/outlet boundary conditions are necessary. Substituting x = 0 and m = 1 in Eqs. (3) and (4) yields

u(L, y ) = u(0, y );

(5)

P (L, y ) = P (0, y ) − PL ;

(6)

which provide inlet-outlet relations for the flow velocity and pressure. Patankar et al. [6] suggested to split the fluid pressure P into two components: a global pressure term −PL x/L and a local re: duced pressure P

P (x, y ) = −

PL L

(x, y ). x+P

(7)

, the momentum equation Eq. (2) is Using this reduced pressure P rewritten to

∂u ∇ P PL + ∇ · (uu ) = − + ν∇ 2 u + , ∂t ρ ρL

(8)

and Eq. (6) now changes to a perfect periodic boundary condition:

(L, y ) = P (0, y ). P

(9)

Eqs. (8) and (1) now can be solved for the one-unit domain as shown in Fig. 1b by an appropriate numerical method, such as the finite difference or finite volume methods, with the no-slip boundary condition on the fluid boundaries  f and  i and the periodic

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

3

Fig. 1. Two-dimensional schematic illustrations of (a) a periodic flow passage with conjugate heat transfer between solid (s , grey area) and fluid (f ) domains, and (b) one periodic unit of length L of the system with inlet (x = 0) and outlet (x = L). The vertical dashed lines in (a) divide the flow passage into individual periodic modules. The coordinate system is set with the x direction along the flow direction and the y direction in the transverse direction. Several types of boundaries and interface are involved in this system: the system boundaries  s and  f separating the solid and fluid domains, respectively, from the outside environment; the solid-fluid interface  i ; and the inlet and outlet (indicated by superscripts) boundaries of the solid and fluid (indicated by subscripts) domains for the isolated unit in (b).

boundary conditions Eqs. (5) and (9) at the inlet and outlet. The last term in Eq. (8) is often incorporated as a body force when solving the flow field [6,9,11,16,22]. 2.2. Temperature decomposition method with conjugate interface The energy equations for the temperature distribution T in the fluid and solid domains are written respectively as follows:

∂ Tf + u · ∇ Tf = α f ∇ 2 Tf , ∂t ∂ Ts = αs ∇ 2 Ts , ∂t

f ,

(10)

s .

(11)

Here α is the thermal diffusivity, and the subscripts f and s denote the fluid and solid domains respectively. On the system boundaries  f and  s , we consider the general Robin condition as:



 ∂ Tf a f Tf + b f = cf , ∂nf  f   ∂ Ts as Ts + bs = cs . ∂ ns 

(12)

(13)

s

Here the parameters af , bf , cf , as , bs , and cs are all timeindependent constants; however, their values may vary at different locations. The local boundary normal directions are denoted by nf on  f and ns on  s , pointing inward into the fluid (f ) or solid (s ) domain (Fig. 1b). With proper values assigned to these parameters, all typical thermal boundary conditions, including the temperature condition, the heat flux condition and the convective condition, can be incorporated. For example, an adiabatic condition can be obtained on the fluid wall by setting a f = c f = 0 and b f = 1; and an isothermal wall at Tw on the solid wall can be achieved with as = 1, bs = 0, and cs = Tw . Across the solid-fluid interface  i , the temperature and energy flux are usually assumed to be continuous. The energy flux continuity is a consequence of the energy conservation principle; and it is usually true except for situations with chemical reactions or the electromagnetic effects at interfaces. For solid-solid interfaces, possible thermal resistance due to the surface roughness could cause a temperature jump at the interface [1,23]. Moreover, studies suggested thermal resistance between solid-liquid and even liquidliquid interfaces may have a significant effect in microscopic transport systems [24,25]. To incorporate the possible temperature discontinuity at the solid-fluid interface  i , we consider the following conjugate condition [26,27]:





Ci T f − Ts = κs

∂ Tf ∂ Ts = κf , ∂ ni ∂ ni

i

(14)

Here ni is the normal direction at the interface  i pointing into the fluid side (Fig. 1b), Ci is the thermal interface conductance, and κ f and κ s are the conductivities of the fluid and solid, respectively. A large Ci value corresponds to the situation with no temperature jump across the interface. To solve the energy equations Eqs. (10) and (11) for periodic thermal flows in one unit as shown in Fig. 1b, we need to add boundary conditions for Tf and Ts at the inlet x = 0 and outlet x = L, which however are not available. To make the solutions of Tf and Ts in one periodic unit possible, we follow the temperature decomposition method (TDM) by Li and Zhang [14,15], and split the regular temperature T into two parts: the transient part θ and the equilibrium part γ :

T f (x, y ) = aθ f (x, y ) + γ f (x, y ), Ts (x, y ) = aθs (x, y ) + γs (x, y ),

f , s .

(15) (16)

Please note that here the word transient is used for the temperature component that decays along the flow passage in space; and this is different from the temporal temperature variations in unsteady flows. The parameter a is constant for one specific periodic module and its magnitude decays exponentially from one module to the next one downstream (See Eq. (37)). More details on the numerical calculation and physical meaning of this coefficient will be presented later in this section. The governing equations for these four temperature components are the same as those for the regular temperatures Tf and Ts . For boundary conditions, we drop the non-homogeneous terms cf and cs for the transient temperature θ , however, make no change for the equilibrium temperature γ . The conjugate conditions at the interface  i for θ and γ also remain the same as Eq. (14) for temperature T. The complete equation and boundary condition systems are

∂θ f + u · ∇ θ f = α f ∇ 2θ f , ∂t

f

(17)

∂θs = αs ∇ 2 θs , s ∂t   ∂θ f af θf + bf = 0, ∂nf 

(18)

(19)

f

 a s θs + b s

Ci



∂θs ∂ ns

 = 0,

(20)

s

 ∂θ ∂θ θ f − θs = κs s = κ f f , ∂ ni ∂ ni

i

(21)

4

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

for the transient temperature θ , and

∂γ f + u · ∇ γ f = α f ∇ 2γ f , ∂t

f

(22)

∂γs = αs ∇ 2 γs , s ∂t   ∂γ f af γf + bf = cf , ∂nf  f   ∂γs as γs + bs = cs , ∂ ns 

(23)

2.3. Determinations of decay rate λL and transient factor a

(24)

(25)

s

Ci



 ∂γ ∂γ γ f − γs = κs s = κ f f , i ∂ ni ∂ ni

(26)

for the equilibrium temperature γ . While it is difficult to establish a relationship for the temperature T between periodic units, we are able to do so for its transient part θ and equilibrium part γ . The θ system is similar to those for the normalized temperature in previous studies of periodic thermal flows with surface heat flux conditions [6,7,9,11]; we therefore propose an exponential decay in θ along the periodic units as in Refs. [14,15]:

θ f (x ± mL, y ) = e∓mλL L θ f (x, y ),

f ,

(27)

θs (x ± mL, y ) = e∓mλL L θs (x, y ),

s .

(28)

The parameter λL is the decaying rate that describes the overall θ variation in the streamwise direction [16]. On the other hand, we assume the exact periodic relationship for the equilibrium part γ among periodic units following Refs. [14,15]:

γ f (x ± mL, y ) = γ f (x, y ),

f ,

(29)

γs (x ± mL, y ) = γs (x, y ),

s .

(30)

Substituting x = 0 and m = 1 into the above equations provides the necessary inlet-outlet relations for solving θ and γ in one periodic unit shown in Fig. 1b:

θ f (L, y ) = e−λL L θ f (0, y ) in  f

θs (L, y ) = e−λL L θs (0, y ) in s ; (31)

γ f (L, y ) = γ f (0, y ) in  f

γs (L, y ) = γs (0, y ) in s . (32)

At last, we should note that the mathematical system for θ is underdetermined: both governing equations Eqs. (17) and (18) are linear; and all the boundary constraints Eqs. (19)–(21), and (31) are homogeneous. As a consequence, no particular solution for θ can be obtained even with all parameters in the wall conditions Eqs. (19) and (20) been specified. First θ = 0 is always a solution, which is useless for our present study. Furthermore, if any nonzero solution θ ∗ exists, βθ ∗ (β is an arbitrary factor) will also be a solution to the θ equations and boundary constraints. To fix this problem, we follow the practices in previous studies [6,7,9,11], and set the bulk average value of θ at inlet x = 0 to 1:



θ¯0 =

 inf

requirements for the thermal field, including the governing equations Eqs. (10) and (11), the general thermal boundary conditions Eqs. (12) and (25), and the conjugate condition Eq. (14). This means the temperature from our TDM method is indeed the solution to the original periodic thermal flow system.

u(0, y )θ f (0, y )dy  = 1,  in u (0, y )dy

(33)

f

which makes the mathematical system for θ determined and solvable with appropriate numerical schemes [6,7,11]. Eqs. (15)–(33) set up a complete system for the transient temperature θ and the equilibrium part γ . Based on these equations, one can mathematically verify that the temperature T obtained from θ and γ via the Eqs. (15) and (16) satisfies all the

When solving the transient temperature θ from the above equations, the decay coefficient λL is involved; however, its value cannot be obtained directly from the system parameters. Following previous practices [6,7,11,14,15], we apply the energy conservation law to the one-unit control volume in Fig. 1b and the following expression for λL is derived (See the online Supplementary Document):

⎡ 1 L

λL = − ln ⎣

ρ c p, f

κs 



in f

 s

 ∂θs 



∂ ns d  + κ f  f ∂θ

uθ f − α f ∂ xf



x=0



∂θ f ∂nf



⎤

d

  s dy − κs  in ∂θ ∂ x x=0 dy s

⎦.

(34) This equation is used in the θ calculation process to iteratively update the λL value until a convergence criterion is reached [6,11,16]. For the steady calculations in this paper, we start with 10−2 for the initial λL value, and the system can gradually converge to a steady λL value. Detailed treatments for unsteady flows will be presented in Section 3.3. After obtaining the θ and γ distributions for one periodic unit, next the value of the transient coefficient a needs to be found so that the original temperature T can be calculated via Eqs. (15) and (16). Applying the bulk average operation similar to that in Eq. (33) to Eq. (15), we then have

a0 = T¯0 − γ¯ 0 ,

(35)

where the inlet bulk average values T¯0 for the fluid temperature and γ¯ 0 for the equilibrium temperature are given as



T¯0 =

 inf

u(0, y )T f (0, y )dy   in u (0, y )dy f



γ¯ 0 =

 inf

u(0, y )γ f (0, y )dy  .  in u (0, y )dy

(36)

f

To construct the temperature field in a particular unit, the inlet bulk temperature T¯0 needs to be provided, and the value of γ¯ 0 can be easily calculated from the γ f distribution field from the oneunit simulation. The transient coefficient for this unit, a0 , can then be obtained from Eq. (35). In the m-th periodic unit downstream, the equilibrium part γ remains the same Eqs. 29 and (30), and the transient portion am θ decays to a0 e−mλL L θ (Eqs. (27) and (28)). In other words, we can say that the transient factor am for the m-th unit downstream is related to the factor a0 of the 0-th unit as

am = e−mλL L a0 .

(37)

This relationship allows us to construct the physical temperature fields over a multiple-unit length under different inlet conditions (i.e., different T¯0 values), based on one single computer simulation over the one-unit domain. Clearly as m increases, the magnitude of the coefficient am gets smaller and the temperature value T approaches γ . At the state where am → 0 and therefore T → γ , the thermal equilibrium is established (i.e., the temperature distribution in each module becomes identical), and this is the reason for calling γ the equilibrium temperature component. Moreover, the transient factor a0 can have negative or positive values depending on the inlet bulk temperature T¯0 at the entrance. When T¯0 > γ¯ 0 , then a0 > 0 and this results in a cooling process: temperature decreases with flow as approaching the thermal equilibrium state. On the other hand, if a0 < 0, heating occurs since T¯0 < γ¯ 0 . It should

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

5

Fig. 2. Geometry illustrations of the straight channel (a) and one period of the wavy channel (b) used in this study. For both channels, the boundaries  s and  f correspond respectively to the outer solid and fluid walls, and the solid-fluid interface  i separates the solid (grey area) and fluid (white area) domains. Geometry parameters used to describe the channel structures in the text are also labeled in these figures.

be mentioned that the TDM method is a pure mathematical treatment and it is not related to any particular numerical schemes; and therefore it can be adopted with any computational methods for flow and heat transfer, such as the finite difference method, the finite volume method, the finite element method, and the recent lattice Boltzmann method [28,29]. The general process for utilizing TDM in periodic thermal flow simulations is outlined in the online Supplementary Document.

In addition, we impose uniform velocity and temperature profiles are at the channel inlet:

3. Method validation and numerical simulations

Here we set the coordinate system as the x axis horizontal along the channel and the y axis vertical, with the origin at the left end of the solid-fluid interface (Fig. 2a). The fluid viscosity ν and thermal diffusivity α are adjusted to set the Reynolds number Re = U0 H/ν = 40 and the Prandtl number P r = ν /α = 0.7 [14,15,32]. Clearly, this is a heating process since a cold (T = 0) fluid current enters the system with one sidewall hot (T = 1 at bottom) and another adiabatic (∂ T /∂ y = 0 at top). We would like to emphasis that this long channel system is simulated directly using LBM, and the decomposition method is not involved at all. The simulated results for the streamwise velocity u and temperature T are displayed in Fig. 3(a, b). Here we only show results in the 20H length from the inlet for a clear representation of the flow and temperature field variations, although the simulation is conducted for a channel of total length L = 50H = 2500. To view the overall temperature variation along the entire channel length, we plot the bulk temperature T¯ along the channel in Fig. 3c. It can be seen that, after a length of about the channel width H from the entrance, the flow quickly evolves into a stable state (Fig. 3a). On the other hand, the temperature field is still developing (Fig. 3b), and no apparent patterns or relationships can be observed. The bulk temperature in Fig. 3c exhibits an exponential trend; however, no stable state is established yet even after a long distance of 50H. The flat channel system simulated here can be actually viewed as a very simple periodic system: any position along the channel can be taken as the inlet of a periodic module of any finite length, and the channel is made up by repeating such identical modules along the flow direction. Therefore, the TDM concept and formulations should be applicable to this system as well. To verify this speculation, a polynomial series solution is derived for θ of this system (See the online Supplementary Document). The equilibrium part γ can be relatively easily solved for this simple system as follows. For the equilibrium part γ , the governing equations Eqs. (22) and (23) can be simplified to

In this section, several demonstration simulations are presented to validate our TDM method and illustrate its usefulness and potential feasibility. For simplicity, only two-dimensional (2D) laminar flows are considered. The lattice Boltzmann method (LBM) [14,29] was used in these numerical simulations; however, other computational methods can also be employed with TDM with no difficulty. Detailed descriptions about LBM for flow and heat transfer simulations can be found in the literature [11,29–31]. Due to the absence of publications on periodic thermal flows with conjugate heat transfer, we are not able to validate our method by comparing TDM results to literature data as usual. 3.1. Thermal flow through flat channel: TDM vs. direct simulation In the previous section, we mathematically show that the temperature from TDM is in fact the correct solution for the periodic thermal system; and the physical meanings of TDM elements (θ , γ , λL , and a) are discussed in details. Here, we design a simple, however, thoughtful system to further validate our TDM method numerically. Here, flow and temperature simulations are performed for a 2D flat channel consisting of solid and fluid domains and an interface along the channel (Fig. 2a). The length of the channel is L = 2500, while the interface separates the fluid flow section of width H = 50 from the solid part of width H/2 = 25. The conductivity of the solid part is considered three times larger than that of the fluid, κs /κ f = 3; and the interface conductance is taken as Ci = 0.02, resulting in an appreciable temperature jump at the interface, as we will see in the results. All these properties are nondimensional in lattice units. The thermal conditions at the system boundaries are: the Dirichlet condition with an constant temperature T = 1 at the bottom solid surface, and the Neumann condition with an adiabatic top fluid wall:



∂T ∂y



=0 ;

Ty=−H/2 = 1.

(38)

y=H

For the fluid section, we apply the no-slip condition for velocity

uy=0 = uy=H = 0.

(39)

ux=0 = [U0 , 0] ;

Tx=0 = 0

(40)

and no-gradient conditions at the outlet:

 ∂ u  =0 ; ∂ x x=L

∂ 2γ f =0 ; ∂ y2

 ∂ T  = 0. ∂ x x=L

∂ 2 γs =0 ∂ y2

(41)

(42)

considering the periodicity in Eqs. (29) and (30) for this steady system with the classical Poiseuille flow in the channel. The boundary

6

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

Fig. 3. Velocity (a) and temperature (b) distributions from the direct long-channel simulation, and the bulk temperature T¯ (c, dark solid line) along the fluid part of the channel. Also, in (c) four prediction curves (dashed lines), starting from x∗ /H = 6, 18, 30, and 42, are plotted using the bulk temperature at the starting position and the analytical decaying rate λL = 8.969381 × 10−4 . These curves at the channel end are enlarged in the inset of (c) to show the prediction accuracy.

conditions for γ are the same as those for the original temperature T in Eq. (38):



∂γ f ∂y



=0 ;

γs(y=−H/2) = 1.

(43)

y=H

Since the general solutions to Eq. (42) are linear functions with constant gradients in the y direction, the adiabatic condition at the top boundary y = H requires γ f to be a constant-value function with no gradient in the fluid domain. Accordingly, γ s in the solid domain also needs to have a constant value to satisfy the heat flux requirement in the conjugate condition Eq. (26). Meanwhile, the zero γ gradients at the interface suggest γs = γ f across the interface. Finally, the boundary condition γs = 1 at y = −H/2 fixes the constant value for γ at 1, and therefore we have

γ f = γs = 1.

(44)

for this system. Obviously the bulk equilibrium value γ¯ 0 is also 1, and this should be the utmost value for the bulk temperature T¯ in Fig. 3c, if the channel is long enough till the thermal equilibrium state is established. At a particular position x∗ along the channel, if we take this position as the inlet of a periodic unit of length l = x − x∗ (x is a location downstream from x∗ ) and apply the TDM treatment here, the corresponding transient factor at this position a(x∗ ) can be obtained from the bulk temperature T¯ (x∗ ) as

a(x∗ ) = T¯ (x∗ ) − 1;

(45)

and the transient factor a(x) at the exit position x should be





∗ a(x ) = e−λL (x−x ) T¯ (x∗ ) − 1 ,

(46)

according to Eq. (37). The bulk temperature at position x can then be expressed as





∗ T¯ (x ) = a(x ) + 1 = e−λL (x−x ) T¯ (x∗ ) − 1 + 1, x ≥ x∗ .

(47)

The decaying rate in this system is λL = 8.969381 × 10−4 from our series solution. This relation can be used to predict the bulk temperature T¯ (x ) downstream from that at the current position x∗ , if the TDM theory reflects the variation features along the channel.

To test this statement, the predicted bulk temperature variations from Eq. (47) are also plotted in Fig. 3c as dashed lines with the starting locations x∗ = 30 0, 90 0, 150 0, and 210 0 shown as filled circles. Clearly Eq. (47) predicts the downstream bulk temperature fairly accurate. For a more quantitative comparison, the simulated and predicted T¯ profiles at the channel outlet are enlarged in the inset of Fig. 3c. The maximum difference is found between the prediction from leftmost starting position x∗ = 6H = 300 (bottom dashed line) and the direct simulation (top solid line), with the bulk temperature values at the exit x = 50L as 0.8983 and 0.9056, respectively. The relative difference is only 0.81%. This is quite impressive, considering that the prediction is made from a relatively close location to the inlet (x∗ /H = 6) and throughout a long distance of 44H. In addition, with series solution of θ available, we can easily construct the temperature profile from TMD, using the transient factor a(x) from the local bulk temperature. Fig. 4 compares the temperature profiles from our direct LBM simulation (blue dots) and those constructed with TDM (red dashed lines) at different locations along the channel. At the conjugate interface (horizontal dashed lines), the solid temperature profiles appear steeper than the fluid profiles, and this observation is consistent to the three times larger solid conductivity in our system. The temperature discontinuity is also noticeable at the interface: With the bottom solid hotter and top fluid cooler, the heat flux flows from solid to fluid across the interface, and this results in a temperature drop from Ts to Tf because of the low interface conductance Ci = 0.02 here. At the later stage of the process (take Fig. 4f as an example), the temperature field has evolved close to the equilibrium state γ = 1, and the energy flux across the conjugate interface is small. Accordingly, the temperature difference across the interface also reduces and becomes less significant. The relative error between the simulated and predicted profiles are also calculated as



E2 =

(TLBM − TT DM )2 dy 

TT2DM dy

.

(48)

Here, TLBM represents the temperature from our direct LBM simulation, and TTDM is the prediction from our TDM method based on the local bulk temperature; and the two integrals cover both

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

7

Fig. 4. Comparison of temperature profiles from our direct long-channel simulation (blue dotted lines) and the predictions according to our decomposition method Eqs. (15) and (16) (red dashed lines) at several axial locations. The black dashed lines indicate the conjugate interface between lower solid and upper fluid domains. Also displayed in each panel is the relative difference E2 between the simulated and predicted profiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the solid and fluid domains. The calculated error values are also displayed in panel titles Fig. 4. It is clear from both the plotted profiles and the relative errors that our TDM method can describe the temperature variation in this channel accurately, except the entrance part x < 4H. This is reasonable since the application of TDM requires both the flow and thermal fields to be well developed from the entrance as in typical periodic thermal flow studies [6,7,9–11]. These analysis and comparison demonstrate that our TDM method can indeed describe the temperature variation in periodic flows with conjugate interfaces. 3.2. Wavy channel flow: one- vs. three-module TDM simulations Wavy channels are widely found in heat transfer applications and have been extensively investigated in previous work, most of which are limited to systems with uniform wall temperature [8– 10,16,33]. Recently, the temperature decomposition method was used to study similar systems under variable thermal wall conditions [14,15]. However, wavy channels with conjugate solid-fluid interface have not yet been investigated. Fig. 2b shows the geometry of the wavy channel used in our simulation. It consists of a fluid flow passage with varying width described by the following equation:

H (x ) = Havg − 2A cos

2π x  L

,

(49)

and a solid part at the bottom with constant thickness Hs , to incorporate the conjugate interface. Here, H(x) is the local channel width for the flow passage, Havg is the average width over a periodic unit, and A is the wavy amplitude. Now, the maximum channel width Hmax = Havg + 2A occurs at the middle x = L/2, and the

minimum width Hmin = Havg − 2A is found at the inlet x = 0 and outlet x = L. The solid wall thickness Hs is set as Hs = L/6. Following previous studies, we use Havg /L = 13/28 and A/L = 1/8. The Reynolds number is Re = U0 Havg /ν = 100 and the Prandtl number is P r = ν /α = 0.7 [8,10,11]. In our calculation, the flow is driven by a body force term (the last term in Eq. (8)), which is constant on all fluid nodes and its magnitude is adjusted during the computation to achieve the required Reynolds number. A temperature discontinuity at the solid-fluid interface in introduced by using a relatively low conductance Ci = 0.005, and the solid-fluid conductivity ratio is set as κs /κ f = 3. For this proposed system, we first consider a relatively simple situation with uniform temperature and heat flux on the walls:



Tbottom = 1 ;

∂T ∂n



= −0.02.

(50)

top

In this case, a constant heat flux is being transferred to the fluid since the top wall is under negative normal gradient. The bottom solid wall is kept at constant temperature. With the underlying periodic relationships correctly considered in a numerical scheme, in principle one can conduct the periodic thermal flow simulation with any positive numbers of periodic modules. Obviously the one-module would be most desirable and efficient for steady laminar flows; however, by comparing the results from simulations with different numbers of modules, we can numerically verify that the periodic thermal relations described above are consistent to the governing equations and boundary condition constraint. This comparison approach has been adopted in previous studies [10,11,15]. Two simulations are conducted separately: one includ-

8

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

Fig. 5. Comparisons of the calculated flow streamlines (a), transient temperature θ (b), and equilibrium temperature γ (c) from the one- and three-module simulations. The one-module results are displayed in red in (a) and as contour lines in (b) and (c), while the three-module results are shown in blue in (a) and as background color patches in (b) and (c). The vertical dashed lines are displayed to indicate individual periodic modules. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

ing only one module, and the other including three modules. The length of one period is L = 120; and therefore the domain lengths in the flow direction are L = 120 for the one-module simulation and 3L = 360 for the three-module simulation, respectively. All these dimension values are in lattice units. Results of the streamlines and the two temperature components, the transient temperature θ , and the equilibrium temperature γ , are shown in Fig. 5. For the fluid part, the streamlines from the one-module simulation are displayed in red in Fig. 5a, while those from the threemodule simulation are in blue. The gray patch represents the solid part and there is no flow there. The flow patterns appear identical between the two simulations and among the three units, suggesting the periodic flow condition has been correctly and accurately incorporated in our calculations. Fig. 5b displays the contour plots for the transient temperature θ in the fluid and solid domains, with the color patches from the three-module simulation and the

contour lines from the one-module simulation. The contour lines in the second and third periods were obtained by multiplying the θ results from the one-module calculation with coefficients e−λL L and e−2λL L according to Eqs. (27) and (28), with the decaying rate λL = 8.192918 × 10−4 . The λL value from the three-module calculation is identical up the last digit. Similarly, Fig. 5c displays the distributions and variations for the equilibrium temperature γ from both calculations. The contour lines in the first module (0 ~ L) are directly from the one-module simulation, and those in the second (L ~ 2L) and third (2L ~ 3L) modules are simply duplication of those in the first module, considering the exact periodic inlet-outlet relationship for γ f and γ s in Eqs. (29) and (30). Excellent agreement between the contour lines and the color patches, in both the solid and fluid domains and for both θ and γ , is noticed in Figs. 5b and c. This observation, as well as the nearly identical λL values from the two separate simulations, strongly suggest that our TDM

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

method can reflect the underlying variation features of temperature along periodic thermal flows. Next we perform a close examination on the temperature variations in Fig. 5b and c. With a higher thermal conductivity in the solid domain (κs /κ f = 3), the temperature variation across the channel is more abrupt in the fluid part and more gentle in the solid part, as indicated by both the color changes and the gaps between contour curves. The isothermal condition Tbottom = 1 in Eq. (50) implies θ = 0 and γ = 1 on the bottom solid surface; these conditions, in combination with the good solid conductance, lead to an approximately constant temperature field with minimal variations in the solid part. On the other hand, the uniform heat flux condition (∂ T /∂ n )top = −0.02 in Eq. (50) introduces an adiabatic condition for θ ((∂ θ /∂ n )top = 0) and a uniform influx for γ ((∂ γ /∂ n )top = 1). As a result, we see θ contours perpendicular to the top wall in Fig. 5b, but the γ contours intersect with the top wall at similar angles in Fig. 5c. The transient component θ represents the temperature portion which decays exponentially; and this is consistent to the similar spatial distributions, however with decreasing magnitude, in the three periodic units. As for the equilibrium component γ , it becomes identical in each periodic unit. The hot spot occurs around the center of the top wall, since here the fluid channel width is the maximum and the heat transfer path from the top to the bottom walls is the longest. For both the θ and γ distributions, due to the relatively low conductance Ci = 0.005 assigned to the solid-fluid interface, a temperature discontinuity can be observed, shown as the mismatch in color patches and contour lines across the interface. Now with the values of θ and γ obtained in both fluid and solid domains, the physical temperature T along the wavy channel can be calculated as per Eqs. (15) and (16), provided that the inlet bulk fluid temperature T¯0 is available and the transient coefficient a can then determined. In our next demonstration, we set T¯0 = 4, which results in a cooling process along the channel since the equilibrium temperature γ f varies in a much lower range between 1 and 2.8 (see Fig. 5c). The bulk values at the fluid inlet γ¯ 0 calculated from the one- and three-module simulations are identical up to the seventh digit as 1.957092. The transient coefficient a0 at inlet x = 0, as per Eq. (35), thus equals T¯0 − γ¯ 0 =2.042908; and those at the inlets to the second (x = L) and third (x = 2L) units are, respectively, a1 = a0 e−λL L =1.85161692 and a2 = a0 e−2λL L =1.67823769. The original temperature T in this 3L long passage can then be calculated from two approaches: (1) using the coefficients a0 , a1 , and a2 with the θ and γ results from the one-module calculation, or (2) using the coefficient a0 only and the θ and γ data from the three-module calculation. The calculated temperature fields T(x, y) from these two approaches are shown in Fig. 6a, similarly with the color patches from the three-module simulation and the contours from the one-module calculation. More quantitatively, the temperature profiles along three horizontal and three vertical positions are compared in Figs. 6b and c, using symbols for results from the one-module simulation and lines for the threemodule results. Again excellent agreement is observed between the one- and three-module simulation results, suggesting that the TDM method has reasonably represented the relationships of temperature variations among periodic structures with conjugate heat transfer. Several similar features noticed in Fig. 5b and c can be found in Fig. 6a as well. Minimal temperature variation is observed in the solid domain due to its high conductivity and small thickness. The isothermal contours appear parallel to the solid-fluid interface and form similar angles near the top wall; and this is consistent to the boundary conditions specified in Eq. (50). The obvious color difference between the two sides of the conjugate interface is associated to the low interface conductance (i.e., large thermal resistance), leading to

9

a temperature jump at that interface. Fig. 6b displays the temperature variation profiles along the flow at three transverse locations: y = L/24 (black), y = 0 (blue), and y = −L/24 (red). Results from the three-module simulation (curves) and the one-module simulation (symbols) match each other impressively well, and both show a decreasing trend which confirms the cooling process as previously mentioned. The temperature discontinuity can be clearly identified in Fig. 6c, where the transverse profiles at the maximum width locations x/L = 0.5, 1.5, and 2.5 are plotted. Across the temperature jump at the interface, the profile slope appears steeper on the fluid side and more gentle on the solid side, which is due to the higher thermal conductivity in solid according to Eq. (14). These transverse profiles are also helpful in verifying that the thermal conditions Eq. (50) at the top (constant gradient) and bottom (constant value) walls are correctly satisfied by our TDM method. Both the contours in Fig. 6a and the profiles in Fig. 6(b, c) clearly show the decreasing temperature trend along the flow direction, and this further confirms the cooling process in this system. 3.3. Unsteady flow with moderate Reynolds number As the Reynolds number increases, the flow in a wavy channel becomes unsteady with periodic, quasiperiodic, and even chaotic fluctuations in fluid velocity [34]. Accordingly, the thermal fields in such unsteady systems also exhibit temporal variations [9,10,16,35]. There are two numerical schemes for generating moderate Reynolds number unsteady flows in the literature, namely the constant body force (or equivalently, constant pressure gradient) approach [9,16,22] and the constant flow rate approach [10,34]. Greiner et al. [35] compared these two methods and found that the results are similar. Here, we conduct a similar calculation to test the performance of our TDM method in the unsteady flow regime, and we adopt the constant body force approach, in which a presumed body force is applied to all fluid nodes. The Reynolds number is obtained from the time-averaged mean velocity after the flow reaches a statistical steady state [16,35]. All the system parameters and wall conditions are kept the same as those in Section 3.2. For a better spatial resolution of the flow and thermal fields, we increase the length of the computational domain to L = 360. To solve the transient component θ , we also follow the treatment in previous simulations of unsteady periodic thermal flows [16,35], and the decay factor λL is calculated from Eq. (34), however, with the integral terms replaced with their respective time-averaged values via the backward moving average [16]. The initial guess value 10−2 is used at the beginning of the calculation before the backward moving average can be conducted. Fig. 7 a displays the temporal variations of flow velocity and temperature components for the center of the periodic module after the statistical steady state is established. Here the velocity has been normalized with the time-averaged mean velocity U¯ 0 and the time has been normalized with Havg /U¯ 0 . The variation patterns appear quasiperiodic as described by Guzman and Amon [34]. The Reynolds number calculated from the resulting mean velocity U¯ 0 in our simulation is Re = 846, and the variation period is estimated as 3.933. Furthermore, instantaneous flow streamlines and thermal fields at four evenly-spaced instants in one variation period are plotted in Fig. 7b. The streamline patterns are similar to those reported in previous publications [10,34]; however, the thermal fields here are not directly comparable to other publication results due to the different thermal conditions on the boundary walls, the presence of the solid domain and the conjugate interface in our system, and moreover that our results are for the transient and equilibrium temperature components and there are no counterparts in previous publications. The features associated with the θ and γ fields discussed in Section 3.2 are still valid in general, such as the isothermal tem-

10

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

Fig. 6. The calculated temperature field via Eqs. (15) and (16) with T¯0 = 4 at the inlet over a 3L long segment of the wavy channel (a), the streamwise variations at y=L/24 (black), 0 (blue), and -L/24 (red) (b), and the variations across the channel width at the middle of each period (c). The background color patches in (a) and the lines in (b, c) are from the three-module simulation, and the isothermal contour lines in (a) and the symbols in (b, c) are from the one-module results. The vertical dashed lines in (a, b) are displayed to indicate individual periodic modules; and those in (c) represent the solid-fluid interface. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

perature at the bottom solid wall, the small gradient in solid due to the larger solid conductivity, the discontinuity at the solid-fluid interface, and the thermal conditions on the top wall (adiabatic for θ and constant flux for γ ). Unlike the smooth curves in Fig. 5, clearly the isothermal contours in Fig. 7b are highly disturbed by the disordered flow field. The hot patch in the transient part θ near the inlet is stretched longer due to the stronger convection current alone the centerline, and the difference between the inlet and exit becomes smaller. This is similar to previous observations for the normalized temperature in periodic thermal flows with isothermal walls [10,22]. The equilibrium temperature on the top heating wall is higher in the unsteady system (~ 3.7) than that in the steady case (~ 2.8). This might also be caused by the strong current along the centerline, which creates a pair of larger circulation vortexes (separation bubbles) in the wide section. To understand the effect of the separation bubbles on the γ distribution, first please recall that the equilibrium component γ represents the thermal equilibrium state of the system far downstream, at which the constant heat injection from the top wall is balanced by the heat release via the isothermal bottom wall. The flow and thermal conditions at the inlet and outlet are exactly identical; and we will use the narrow and wide sections in our next analysis to refer streamwise loca-

tions along the channel. For the steady flow with a low Reynolds number (see the streamlines at Re = 25 in Refs. [10,11] as examples), after leaving the narrow section, the fluid follows the wavy wall geometry, expands in the wide section, and merges back as re-entering the narrow section. There is no separation bubbles, and the heat injection from the top wall in the wide section can be carried out directly by the streamwise current to the narrow section. Similarly, in the lower portion of the system, the fluid current picks up the heat from the upper stream at the narrow section and it them brings it to the cold low boundary (the solid-fluid interface in our system). However, as the Reynolds number increases, a pair of separation bubbles are developed near the wavy boundaries at the wide section [9,10]. These separation bubbles prevent the streamwise current from touching the hot boundary on top and the cold boundary at bottom; instead, the heat transfer across the channel takes place via the following route: [hot wall at top] → [upper separation bubble] → [streamwise current at centerline] → [lower separation bubble] → [cold wall at bottom]. Clearly this less direct pathway reduces the heat transfer efficiency. The negative impact of circulation vortexes on heat transfer is also evident in both the instantaneous and time-averaged plots in Fig. 7b: In the circulation areas, the isothermal contours are less

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

11

Fig. 7. (a) Temporal variations of flow velocity and temperature components at the center of the computation domain (x/L = y/L = 0), and (b) instantaneous streamlines (top row) and isothermal contours for the transient (middle row) and equilibrium (bottom row) temperature components at Reynolds number Re=846. The flow and temperature fields are taken at four evenly-spaced time instants (labeled on top of each column) over one variation period, and the time-averaged streamlines and temperature distributions are displayed in the fifth column. The four time instants are also indicated as vertical dashed lines in (a). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

crowded, indicating gentle temperature variations and poor heat conduction in those regions. For the two Reynolds numbers considered in this work, the unsteady system at Re = 846 has larger separation bubbles than the steady Re = 100 case, and therefore the heat transfer across the channel is less efficient. With a constant temperature at the bottom solid wall and a constant heat injection flux at the top fluid wall, the less efficient heat transfer for Re = 846 results in a higher temperature at the top wall. Please be aware that at the thermal equilibrium state, the fluid is saturated, and it only serves as the transport medium but does not carry any heat away. This is different from typical heat transfer systems, where the fluid is responsible for taking heat away from the structure surfaces; and in those situations, the strong streamwise current at a high Re is preferable since the large flow rate can remove heat faster from the system. 3.4. Effects of solid conductivity and interface conductance At last, as a further example to demonstrate the potential applications of our method, we further investigate the effects of the interface conductance Ci and the solid-fluid conductivity ratio κ s /κ f on the temperature variation. The same geometry shown in Fig. 2b is adopted; however, the thermal conditions on the outer walls are modified to represent more general and realistic situations that may be encountered in heat transfer devices: Along the top fluid wall, the heat injection flux magnitude varies from zero at the narrowest section to the maximum at the widest section:





∂T ∂n





= −0.01 1 − cos

2π x 

top

L

;

(51)

On the bottom solid wall, we consider a convective condition

∂T ∂n



= −0.1(1 − T ),

(52)

bottom

meaning that the bottom solid wall is exposed to a surrounding environment with a reference temperature of T∞ = 1. Four different scenarios are considered in this study; and each case is indi-

cated by a two-letter label: ’k’ for a same solid conductivity as fluid (κs /κ f = 1) or ’K’ for a better solid conductivity (κs /κ f = 3); and ’c’ for a low interface conductance Ci = 0.005 or ’C’ for a high interface conductance Ci = 10 to efficiently remove the temperature discontinuity. Accordingly, case KC corresponds to the system with best conduction and case kc as the worst. The θ and γ results for the four cases kC, kc, KC, and Kc are displayed in Fig. 8. All calculations are conducted using only one periodic unit. The isothermal contour lines for θ in Fig. 8a1-4 appear perpendicular to the top wall; this conforms with the adiabatic boundary condition ∂ θ /∂ n = 0 at the top fluid wall according to Eqs. (51) and (19). On the other hand, in the same figures, the relation between the θ distributions near the bottom solid wall and the applied boundary condition ∂ θ /∂ n = −0.1θ is less apparent, however, apprehensible: Due to the proportionality between θ and its gradient above, along the bottom wall, the light blue color (higher θ value) is associated with dense contour lines (larger θ gradient) at the narrow section, and on the contrary the dark blue color (low θ value) and sparse contours are found at the wide section in the middle. As for the γ f results in Fig. 8b1-3, we see contours form similar angles to the top wall since the same variable heat flux condition ∂ γ /∂ n = −0.01[1 − cos(2π x/L )] is applied for all the four cases as per Eqs. (51) and (24). The convective condition for γ , ∂γ /∂ n = −0.1(1 − γ ), at the bottom solid wall results a similar relation between the γ value and its gradient as we describe above for θ . The hot spot of γ is found near the center of the top wall, and this is caused by the maximum heat injection from the top wall and also the long transverse distance to the bottom wall at the wide section. The solid conductivity effect could be clearly identified by examining the gap distance between contour lines near the interface for both θ and γ . Similar line gaps can be observed in the fluid and solid domains in cases kc and kC, where κs = κ f , indicating similar gradients on the two sides of the interface. When the solid conductivity increases in cases Kc and KC (κs = 3κ f ), smaller gradients on the solid side are expected thus to ensure the continuity of heat flux across the interface (see Eq. (14)). This is evident

12

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

Fig. 8. The transient temperature θ (a1-4) and the equilibrium temperature γ (b1-4) for the four cases with different solid conductivity and interface conductance. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table 1 Characteristic properties for the four cases of wavy channel flows with different solid conductivity and interface conductance.

κs = κ f (k) Parameters decay rate: λL (10−4 ) maximum γ : γmax inlet bulk average: γ¯ 0 Transient factor: a0 90% θ reduction length: n90%

Ci = 10 (kC) 6.624468 2.251705 1.563889 -0.7638889 29.0

in Fig. 8a3-4 for θ and Fig. 8b3-4 for γ , with wider line gaps in solid and narrower gaps in fluid (approximately three times different). The effect of the interface conductance Ci can be observed by comparing the left vs. right panels of the θ and γ plots in Fig. 8. In cases kC and KC, the relatively high conductance value Ci = 10 allows the heat flux to be transferred smoothly across the interface with negligible thermal resistance. Smoother color variations and continuous contour lines are obtained near the interface for both θ and γ . On the other hand, with a low conductance value Ci =0.005 in cases kc and Kc, a temperature jump occurs at the interface due to high resistance, as indicated by the mismatch in color and discontinuity of contour lines on the two sides of interface. For a more quantitative comparison, we list several representative values for the four cases in Table 1: the decaying rate λL for θ , the maximum γ value in the entire periodic unit γmax , and the bulk average of γ at the inlet γ¯ 0 . The relatively small λL values suggest that the process toward the thermal equilibrium state will be slow for all cases, as to be confirmed next in Fig. 9. Among the four cases, the maximum λL = 1.098707 × 10−3 is found in case KC, the best situation with good solid conductivity and interface conductance, and the minimum λL = 5.253905−4 in case kc with low solid conductivity and interface conductance. This relation is reasonable, since case KC is most favorable for heat transfer and thus the equilibrium state can be established quicker; and vise versa for case kc. Cases kC and kC fall in between, with a slightly larger λL for case Kc, meaning the heat transfer efficiency would be slightly better in case Kc than kC for the particular parameters used here. When the thermal equilibrium is reached, the heat

Ci = 0.005 (kc) 5.253905 2.390650 1.701583 -0.9015830 36.5

κs = 3κ f (K) Ci = 10 (KC) 10.98707 2.040396 1.353855 -0.5538551 17.5

Ci = 0.005 (Kc) 7.528128 2.189161 1.501660 -0.7016597 25.5

flux injected into the system from the top wall counterbalanced by the flux exiting the system through the bottom solid wall. Therefore, in a system with better thermal conduction and low interface resistance, heat can be transferred from the top to bottom walls relatively easily, resulting in a lower overall temperature field and a smaller transverse temperature gradient. For this reason, the lowest values among the four cases, both for the maximum γmax and inlet bulk average γ¯ 0 , are found in case KC with best thermal conduction. Again case kc is the worst (with the highest γmax and γ¯ 0 ), and case Kc is marginally better than case kC. These observations and discussions further demonstrate that the conjugate conductance has been correctly considered in our TDM method for periodic thermal flows. As previously mentioned, the advantage and attractiveness of the one-module simulation technologies, such as the TDM method here, are the ability to effectively predict the temperature variation along a flow passage of multiple modules. In this part, we explicitly demonstrate this fact using the θ and γ results in Fig. 8 with different solid domain conductivity and interface conductance combinations. The inlet bulk temperature is set at T¯0 = 0.8; correspondingly the flow process would be a heating process since T¯0 < γ¯ 0 for all cases in Table 1. The transient coefficients a0 calculated from Eq. (35) are also listed in Table 1. The negative values indicate that heating occurs in all these processes. The temperature fields in the first eight periods are also shown in Fig. 9. Please note the four systems here have the same thermal conditions applied on the solid and fluid walls Eqs. (51) and (52), and the same bulk temperature T¯0 = 0.8 is imposed at the

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

13

Fig. 9. The calculated temperature fields in the first eight periods of a wavy channel from one-module calculations with different solid conductivity and interface conductance. The red lines correspond to the isothermal contours with T = T∞ = 1, which separates the heating and cooling regimes on the bottom solid wall. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

inlet. The only difference lies in the different solid conductivity κ s and the interface conductance Ci as described above. The features associated with these different parameters discussed above in Fig. 8 can also be seen here for the temperature T. With the heat injection flux on top (Eq. 51), hot spots appear near the top wall at the wide section locations, where the heat flux is at its maximum. The temperature in these hot spot regions also increases in the flow direction as more heat is added to the fluid. However, the role of the bottom solid wall is less straightforward, since the convection condition Eq. (52) there implies that the solid wall is exposed to a medium with an ambient temperature T∞ = 1. It is interesting to find that the bottom solid wall acts as a heat source to the fluid in the first 4 ~ 5 units, since the fluid is so cold that heat is added to the system from the T∞ = 1 outside medium via the solid wall. After that, the fluid temperature has increased enough and then the heat transfer direction across the bottom solid wall reverses: The solid wall behaves as a heat sink and the fluid releases heat to the medium across the solid wall. To better illustrate the transition between these two regimes, the isothermal line at T = 1 is plotted and is shown in red color in Fig. 9 for all four cases. This line reaches the bottom boundary perpendicularly, since there is no heat flux at that point when the wall temperature is the same as the ambient temperature. Again the earliest source-sink transition occurs at 4.22L in case KC with the best thermal conduction, and the latest 5.04L is found in case kc with a low solid conductivity and a high interface resistance, consistent to our previous observations and discussions.

Clearly the thermal equilibrium states are not established yet in the first eight periodic units in Fig. 9. This is mainly due to the low inlet temperature value T¯0 = 0.8 compared to the equilibrium value γ¯ 0 (1.35 and 1.70 in Table 1), and the small decaying factor λL . To quantify the channel length required to approach thermal equilibrium, we calculate the number of periodic units required from the inlet to have the transient part reduce to 10% of its initial state as

e−n90% λL L = 0.1

→

n90% = −

ln 0.1 , λL L

(53)

and the n90% values are given in Table 1 as well. Consistent to our analysis above, the longest distance of 36.5L is found in the worst conduction case kc and the shortest distance of 17.5L for the best conduction system KC. We would like to have an overall view of how these equilibrium states can be approached; however, simply plotting the temperature fields over this many modules in a figure will hinder the detailed representation of temperature variations. Therefore, for all four cases, we display the centerline fluid temperature over a 60L length in Fig. 10. The equilibrium temperature γ along the centerline of one period is also plotted on the right edge of main panel and in the inset (3) of Fig. 10. The γ curves provide a comparison standard for identifying the equilibrium state. The flows start with an inlet centerline temperature slightly low than the imposed bulk value T¯0 = 0.8 (see the temperature contours at the left channel end in Fig. 9), and the temperature increases monotonically in the first ~ 6 periods, although the increase slope may vary (Fig. 10a). This increasing trend can be better understood by checking Fig. 9 again. In the first 4 ~ 5 periods,

14

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

Fig. 10. Fluid temperature profiles along the centerline y = 0 over a 60-period length of the wavy channel with different solid conductivity and interface conductance. The profiles near the inlet (x/L < 6) are enlarged in inset (1) to show the monotonic increase in the early part of the flow, and inset (2) shows the detailed temperature variation behaviors in each periodic unit. The equilibrium temperature profiles for the four cases considered here are also displayed in inset (3) and the right edge of the main panel for the convenience of our comparison and analysis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

heat fluxes flow into the fluid stream from both the top boundary and the solid interface, causing the fluid temperature to continuously increase. Later on as the temperature on the bottom solid wall exceeds the ambient temperature T∞ = 1, the outside media starts to act as a thermal sink, and such cooling effect is more profound to the fluid in the narrow section where the flow width is small. Furthermore, the heat injection from the top boundary is more intense at the wide section and it drops to zero at the narrow section. This spacial distribution of these top-heating and bottom-cooling effects results in the particular zig-zag variation for x > 6L in Fig. 10 (see inset Fig. 2 for a detailed view). Due to the flow convection, the local maximum and minimum in fluid temperature occur downstream of the widest and narrowest locations. The longer delay in maximum temperature might be attributed to the relatively farther transverse distance from the heat source center on top wall to the centerline, while the cooling solid interface is very close to the fluid centerline. Despite these local zig-zag variations in each period, overall the fluid temperature increases along the flow, and in general it exhibits an exponential fashion and approaches the equilibrium state gradually. As expected, the equilibrium state at relatively low temperature is established earliest in the best conduction case KC, at x =∼ 35L. After that point, the temperature field in case KC becomes identical in each period and also the same as the equilibrium temperature γ . On the contrary, the temperature in the worst conduction case kc is still developing even after a length of 60L. For the two intermediate systems, case Kc outpaces case kC marginally: the equilibrium state is approximately reached at x = 60L in case Kc, but case kC still needs some

extra distance. These temperature profile features are consistent to the quantitative values in Table 1. 4. Summary and concluding remarks In this paper, we have extended the temperature decomposition method (TDM) to periodic thermal flows with conjugate heat transfer involved. The physical temperature is split into two components: the transient part and the equilibrium part. Governing equations and appropriate thermal constraints on the system boundaries, conjugate interface, and the period inlet and outlet are proposed, such that the temperature field along the periodic structures can be obtained from a one-module calculation. In addition to the detailed mathematical descriptions, physical interpretations and numerical implementations are also provided. Several carefully designed example simulations are conducted, and the results are analyzed and discussed in details and depth. The results and discussions clearly show that the TDM method correctly reflects the underlying relationship in temperature, which otherwise is not available, for conjugate periodic thermal systems. These simulations sufficiently verify and demonstrate the validity, accuracy, and potential usefulness of our method. We would like to mention that, although the observations and discussions for the example simulations may appear intuitive and straightforward from fluid mechanics and heat transfer fundamentals, it is not trivial to obtain such information from a numerical approach. The aim of this study is to extend the original TDM method [14,15] to periodic thermal flows with conjugate solid-fluid inter-

M. Jbeili and J. Zhang / International Journal of Heat and Mass Transfer 150 (2020) 119288

faces. For simplicity and clarity, the method descriptions and simulations are limited to 2D laminar systems. The lattice Boltzmann method is used in our simulation, simply to benefit from our existing programs. With these being said, the decomposition method is not related to any particular numerical methods, and it can be implemented in other numerical schemes like finite-element or finite-volume methods as well. For heat transfer systems with turbulent flows and natural convection, numerical techniques available in the literature [7,12,20,21,36] can be utilized in combination with the TDM method in such situations. CRediT authorship contribution statement Mayssaa Jbeili: Data curation, Formal analysis, Writing - original draft. Junfeng Zhang: Conceptualization, Methodology, Supervision, Writing - review & editing. Acknowledgments We thank the anonymous reviewers for carefully reviewing the manuscript and providing critical comments and constructive suggestions. This work was supported by the Natural Science and Engineering Research Council of Canada (NSERC). The calculations have been enabled by the use of computing resources provided by WestGrid (westgrid.ca), SHARCNet (sharcnet.ca), and Compute/Calcul Canada (computecanada.org). Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ijheatmasstransfer.2019. 119288 References [1] F.P. Incropera, D.P. DeWitt, T.L. Bergman, A.S. Lavine, Fundamentals of Heat and Mass Transfer, John Wiley & Sons, New York, 2006. [2] H. Soumerai, Practical Thermodynamics Tools for Heat Exchanger Design Engineers, John Wiley & Sons, New York, 1978. [3] M.A. Ahmed, M.M. Yaseen, M.Z. Yusoff, Numerical Study of Convective Heat Transfer from Tube Bank in Cross Flow Using Nanofluid, Case Studies in Thermal Engineering 10 (2017) 560–569, doi:10.1016/j.csite.2017.11.002. [4] R. Wang, W. Wang, J. Wang, Z. Zhu, Analysis and Optimization of Trapezoidal Grooved Microchannel Heat Sink Using Nanofluids in a Micro Solar Cell, Entropy 20 (2018) 9. [5] M.A. Ahmed, M.Z. Yusoff, K.C. Ng, N.H. Shuaib, Numerical Investigations on the Turbulent Forced Convection of Nanofluids Flow in a Triangular-Corrugated Channel, Case Studies in Thermal Engineering 6 (2015) 212–225. [6] S.V. Patankar, C.H. Liu, E.M. Sparrow, Fully Developed Flow and Heat Transfer in Ducts Having Streamwise-Periodic Variations of Cross-Sectional Area, ASME J. Heat Transf. 99 (1977) 180. [7] E. Stalio, Direct Numerical Simulation of Heat Transfer Enhancing Surfaces, University of Bologna, Bologna, Italy, 2002 Ph.D. thesis. [8] H.M.S. Bahaidarah, N.K. Anand, H.C. Chen, Numerical Study of Heat and Momentum Transfer in Channels with Wavy Walls, Num. Heat Transf. A 47 (2005) 417. [9] G. Wang, S.P. Vanka, Convective heat transfer in periodic wavy passages, International Journal of Heat and Mass Transfer 38 (1995) 3219.

15

[10] A.G. Ramgadia, A.K. Saha, Numerical study of fully developed flow and heat transfer in a wavy passage, International Journal of Thermal Science 67 (2013) 152. [11] Z. Wang, H. Shang, J. Zhang, Lattice boltzmann simulations of heat transfer in fully developed periodic incompressible flows, Physical Review E 95 (2017) 63309. [12] D. Angeli, E. Stalio, M.A. Corticelli, G.S. Barozzi, A fast algorithm for Direct Numerical Simulation of natural convection flows in arbitrarily-shaped periodic domains, Journal of Physics: Conference Series 655 (2015) 12054. [13] A.M. Adham, N. Mohd-Ghazali, R. Ahmad, Thermal and hydrodynamic analysis of microchannel heat sinks: a review, Renewable and Sustainable Energy Reviews 21 (2013) 614–622. [14] P. Li, J. Zhang, Simulating heat transfer through periodic structures with different wall temperatures: a temperature decomposition method, Journal of Heat Transf. 140 (2018) 112002. [15] P. Li, J. Zhang, The Temperature Decomposition Method for Periodic Thermal Flows with General Wall Conditions, Numerical Heat Transfer B 74 (3) (2018) 559–577. [16] E. Stalio, M. Piller, Direct numerical simulation of heat transfer in converging— diverging wavy channels, Journal of Heat Transfer 129 (2007) 769. [17] G. Russ, H. Beer, Heat transfer and flow field in a pipe with sinusoidal wavy surface-II. experimental investigation, International Journal of Heat and Mass Transfer 40 (1997) 1071–1081. [18] Y. Yang, P. Chen, Numerical Optimization of Turbulent Flow and Heat Transfer Characteristics in a Ribbed Channel, Heat Transf. Eng. 36 (2015) 290–302. [19] Y. Peles, A. Kosar, C. Mishra, C. Kuo, B. Schneider, Forced Convective Heat Transfer across a Pin Fin Micro Heat Sink, Int. J. Heat Mass Transf. 48 (2005) 3615–3627. [20] K.M. Kelkar, D. Choudhury, Numerical Prediction of Periodically Fully Developed Natural Convection in a Vertical Channel with Surface Mounted Heat Generating Blocks, Int. J. Heat Mass Transf. 36 (1993) 1133–1145. [21] M. Piller, S. Polidoro, E. Stalio, Multiplicity of Solutions for Laminar, Fully-Developed Natural Convection in Inclined, Parallel-Plate Channels, Int. J. Heat Mass Transf. 79 (2014) 1014–1026. [22] B. Niceno, E. Nobibe, Numerical Analysis of Fluid Flow and Heat Transfer in Periodic Wavy Channels, Int. J. Heat Fluid Flow 22 (2001) 156–167. [23] J.P. Holman, Heat Transfer, McGraw-Hill, New York, 1968. [24] S. Maruyama, T. Kimura, A Study on Thermal Resistance over a Solid-Liquid Interface by the Molecular Dynamic Method, Therm. Sci. Eng. 7 (1999) 63–68. [25] H.A. Patel, S. Garde, P. Keblinski, Thermal resistance of nanoscopic liquid-liquid interfaces: dependance on chemistry and molecular architecture, Nano Lett. 5 (2005) 2225–2231. [26] G. Le, O. Oulaid, J. Zhang, Counter-Extrapolation Method for Conjugate Interfaces in Computational Heat and Mass Transfer, Physical Review E 91 (2015) 33306. [27] K. Guo, L. L, G. Xao, N. AuYeung, R. Mei, Lattice Boltzmann method for conjugate heat and mass transfer with interfacial jump conditions, Int. J. Heat Mass Transf. 88 (2015) 306–322. [28] J. Wendt, Computational Fluid Dynamics: An Introduction, Springer, Berlin, Germany, 2009. [29] J. Zhang, Lattice Boltzmann Method for Microfluidics: Models and Applications, Microfluidics and Nanofluidics 10 (2011) 1–28. [30] Z. Guo, C. Shu, Lattice Boltzmann Method and Its Applications in Engineering, World Scientific Publishing, Singapore, 2013. [31] Q. Chen, X. Zhang, J. Zhang, Improved treatments for general boundary conditions in lattice boltzmann method for convection-diffusion and heat transfer processes, Physical Review E 88 (2013) 33304. [32] G.M. Brown, Heat or Mass Transfer in a Fluid in Laminar Flow in a Circular or Flat Conduit, AIChE Journal 6 (1960) 179. [33] C. Nonino, S. Savino, S.D. Giudice, FEM for the 3-D analysis of conjugate conduction-convection heat transfer in cross-flow microheat exchangers, Int. J. of Num. Method Heat Fluid Flow 25 (2015) 1322–1339. [34] A.M. Guzman, C.H. Amon, Transition to chaos in converging-diverging channel flows: ruelle-Takens-Newhouse scenario, Phys. Fluid 6 (1994) 1994–2002. [35] M. Greiner, R.J. Faulkner, V.T. Van, H.M. Tufo, P.F. Fischer, Simulations of ThreeDimensional Flow and Augmented Heat Transfer in a Symmetrically Grooved Channel, ASME J. Heat Transf. 122 (20 0 0) 653. [36] A.K. Saha, S. Acharya, Unsteady Simulation of Turbulent Flow and Heat Transfer in a Channel with Periodic Array of Cubic Pin-Fins, Numerical Heat Taansfer A 46 (2004) 731–763.