# Numerical study of heat transfer over banks of rods in small Reynolds number cross flow

Gabriel Gamrat, Michel Favre-Marinet, Stéphane Le Person

► To cite this version:

Gabriel Gamrat, Michel Favre-Marinet, Stéphane Le Person. Numerical study of heat transfer over banks of rods in small Reynolds number cross flow. International Journal of Heat and Mass Transfer, 2007, pp.853-423. 10.1016/j.ijheatmasstransfer.2007.04.038. hal-00265920

HAL Id: hal-00265920

https://hal.science/hal-00265920v1

Submitted on 4 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Numerical study of heat transfer over banks of rods in small
Reynolds number cross-flow

Gabriel Gamrat, Michel Favre-Marinet *, Stéphane Le Person

Laboratoire des Ecoulements Ge'ophysiques et Industriels, CNRS-UJF-INPG, 1025 rue de la Piscine, BP 53 X, 38041 Grenoble Cedex, France

This work presents numerical computations of heat transfer over banks of square rods in aligned and staggered arrangements with porosity in the range 0.44–0.98. It is focused on low Reynolds number flows (0.05–40). Two thermal boundary conditions were investigated, namely constant wall temperature and constant volumetric heat source. The effects of bank arrangements and porosity as well as the effects of Prandtl and Reynolds numbers on the Nusselt number are examined. In the case of constant volumetric heat source, the results are approximated with a power equation adapted for the case of low *Re* number flows. This study shows that the thermal boundary condition on the solid surface influences heat transfer when thermal equilibrium is reached in the bank of rods.

Keywords: Bank of rods; Laminar flow; Thermal boundary condition; Thermal equilibrium

**1. Introduction**

Much work has been done in the past on convective heat transfer in banks of tubes or rods in cross-flow. One of the most extensive reviews in the field of cross-flow heat exchanger is that of Zukauskas [1], who proposed correlations between the Nusselt, Reynolds and Prandtl numbers for various arrangements of cylindrical tube banks. These correlations are available for moderate to high values of the Reynolds number (1 < *Re* < 2 × 10⁶).

On the other hand, banks of rods have very often been
used as a geometrical model for low Reynolds number
flows through porous media. Spatially periodic models
have been considered to compute the permeability of the
medium as a function of porosity and Reynolds number
[2–10]. These models are very attractive for numerical sim-
ulations, since the computations may be restricted to a sim-

ple cell extracted from the periodic pattern. The role of
finite Reynolds number flow and the deviation due to
non-linearities from the original Darcy's law have been
extensively discussed in the literature. There are much less
numerical works on heat transfer over banks of rods in low
Reynolds number cross-flow [11–16]. In the context of por-
ous media, one of the issues is that of local thermal equilib-
rium of the fluid and the solid matrix constituting the
porous medium. This problem is much more complex than
the isothermal one, since heat transfer not only depends on
the porosity and the Reynolds number, but also on the Pra-
ndtl number and on the thermal conditions on the solid
surfaces. The range of parameters and boundary condi-
tions found in [11–16] are shown in Table 1.

The objective of our work was to establish a database
for the heat exchange coefficient in banks of squared rods
with the thermal condition of uniform volume source heat-
ing and for low Reynolds number flows. The motivations
were twofold. Firstly, the thermal condition in heat
exchangers is often neither uniform flux nor uniform tem-
perature heating. The influence of this condition on the
heat transfer coefficient is negligible in turbulent flows,
but may be significant for low Reynolds number flows,

<footnote>* Corresponding author. Tel.: +33 (0) 4 76 82 50 49; fax: +33 (0) 4 76 82 52 71.

*E-mail addresses:* Gabriel.Gamrat@hmg.inpg.fr (G. Gamrat), Michel. Favre-Marinet@hmg.inpg.fr (M. Favre-Marinet), Stephane.Leper- son@hmg.inpg.fr (S. Le Person).</footnote>
Nomenclature

<table><tr><td>a<sub>sf</sub></td><td>interfacial surface area per unit length, m</td></tr><tr><td>C<sub>d</sub></td><td>drag coefficient</td></tr><tr><td>C<sub>F</sub></td><td>Forchheimer coefficient</td></tr><tr><td>c<sub>p</sub></td><td>Specific heat at constant pressure, J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>C<sub>r</sub></td><td>resistance factor</td></tr><tr><td>Da</td><td>Darcy number</td></tr><tr><td>d</td><td>solid element size, m</td></tr><tr><td>e</td><td>channel width, m</td></tr><tr><td>h</td><td>convective heat transfer coefficient, W m<sup>-2</sup> K<sup>-1</sup></td></tr><tr><td>K</td><td>permeability, m<sup>2</sup></td></tr><tr><td>K<sub>app</sub></td><td>apparent permeability, m<sup>2</sup></td></tr><tr><td>k</td><td>thermal conductivity, W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>L</td><td>longitudinal and transversal pitch, m</td></tr><tr><td>M&#775;</td><td>mass flow rate per unit length, kg s<sup>-1</sup> m<sup>-1</sup></td></tr><tr><td>Nu<sub>2e</sub></td><td>Nusselt number based on the mean fluid temperature (Eq. (17))</td></tr><tr><td>Nu<sub>2e,b</sub></td><td>Nusselt number based on the bulk fluid temperature</td></tr><tr><td>Nu<sub>2e,x</sub></td><td>local Nusselt number based on the bulk fluid temperature</td></tr><tr><td>Nu<sub>d</sub></td><td>Nusselt number (Eq. (11))</td></tr><tr><td>Pe</td><td>Peclet number (=Re<sub>2e</sub>Pr)</td></tr><tr><td>Po</td><td>Poiseuille number</td></tr><tr><td>Pr</td><td>Prandtl number</td></tr><tr><td>q<sub>v</sub></td><td>volumetric heat source, W/m<sup>3</sup></td></tr></table>

<table><tr><td>Re<sub>2e</sub></td><td>Reynolds number in a channel (Eq. (18))</td></tr><tr><td>Re<sub>d</sub></td><td>Reynolds number (Eq. (9))</td></tr><tr><td>Re<sub>D</sub></td><td>Darcian Reynolds number (Eq. (10))</td></tr><tr><td>T</td><td>temperature, K</td></tr><tr><td>u<sub>D</sub></td><td>Darcy velocity, m s<sup>-1</sup></td></tr><tr><td>x*</td><td>dimensionless distance (= x/2e P<sub>e</sub>)</td></tr></table>

Greek symbols

<table><tr><td>ε</td><td>porosity</td></tr><tr><td>ρ</td><td>density, kg m<sup>-3</sup></td></tr><tr><td>μ</td><td>dynamic viscosity, kg m<sup>-1</sup> s<sup>-1</sup></td></tr><tr><td>φ</td><td>heat flux density, W m<sup>-2</sup></td></tr><tr><td>σ<sub>0</sub></td><td>dimensionless temperature</td></tr><tr><td>θ</td><td>dimensionless temperature</td></tr></table>

Subscripts

<table><tr><td>b</td><td>bulk</td></tr><tr><td>D</td><td>Darcy</td></tr><tr><td>f</td><td>fluid</td></tr><tr><td>max</td><td>maximal</td></tr><tr><td>min</td><td>minimal</td></tr><tr><td>p</td><td>pressure</td></tr><tr><td>s</td><td>solid</td></tr><tr><td>w</td><td>wall</td></tr></table>

tions for the drag coefficient and the heat exchange coefficient of a cylinder in two-dimensional cross-flow. Following this approach, Bavière et al. [20] considered a rough-wall consisting in periodically distributed parallel-epipeds of square cross-section. They estimated the drag coefficient by using the formula for the drag force on very slender prolate spheroids in creeping flows. Their work was restricted to isothermal flows and is currently being extended to improve the determination of the drag coefficient and to take into account heat transfer in the micro-channel. The present paper is therefore devoted to numerical computations of the flow and heat exchange in banks of rods of square cross-section heated by volume

which are predominant in the field of microheat transfer. It is then important to test the sensitivity of the heat transfer coefficient to the thermal condition for the design of micro-heat exchangers. Additionally, the situation of uniform volume source heating is encountered in experimental works on arrays of cylinders with cross-flow convection where the cylinders are electrically heated at uniformly distributed rate [17].

Secondly, we are developing a numerical model for roughness effects on microchannel flows using a discrete-element method initially proposed by Taylor et al. [18,19] for predicting the rough-wall skin friction and heat transfer coefficient in turbulent flows. This method needs correla-

<table><caption>Table 1
Conditions from the literature for heat transfer computations in periodic arrays of rods</caption><thead><tr><th>References</th><th>Rod cross-section</th><th>Geometrical arrangement</th><th>Re</th><th>Pr</th><th>Porosity</th><th>Thermal conditions</th></tr></thead><tbody><tr><td>Martin et al. [11]</td><td>Cylindrical</td><td>Aligned in squared or triangular arrays</td><td>3–160</td><td>0.72</td><td>0.8–0.99</td><td>Constant wall temperature or heat flux</td></tr><tr><td>Kuwahara et al. [12]</td><td>Squared</td><td>Staggered</td><td>2 × 10<sup>-3</sup>–10<sup>3</sup></td><td>10<sup>-2</sup>–10<sup>2</sup></td><td>0.36–0.91</td><td>Constant wall temperature</td></tr><tr><td>Ghosh Roychowdhury et al. [13]</td><td>Cylindrical</td><td>Aligned or staggered</td><td>40–1000<sup>a</sup></td><td>Not given</td><td>0.5–0.8</td><td>Constant wall temperature</td></tr><tr><td>Mandhani et al. [14]</td><td>Cylindrical</td><td>Staggered</td><td>1–500</td><td>0.1–10</td><td>0.4–0.99</td><td>Constant wall temperature or heat flux</td></tr><tr><td>Nakayama et al. [15]</td><td>Squared</td><td>Aligned or cylinders in yaw</td><td>10<sup>-2</sup>–6×10<sup>3</sup></td><td>1</td><td>0.25–0.875</td><td>Constant wall temperature</td></tr><tr><td>Saito and de Lemos [16]</td><td>Squared</td><td>Staggered</td><td>4–400</td><td>1</td><td>0.44–0.9</td><td>Constant wall temperature</td></tr></tbody></table>

<footnote>a Re is based on the mean velocity at the minimum cross-sectional area in [13].</footnote>
sources uniformly distributed inside the rods. Computations were also performed with the condition of uniform temperature heating for comparison with published data.

## 2. Numerical model

### 2.1. Computation domain, physical equations and boundary conditions

The geometrical pattern considered in this study consists of infinitely long rods of square cross-section periodically distributed either in the aligned or in the staggered arrangement (Fig. 1). This configuration is characterized by the side *d* of the square solid elements and by equal transverse and longitudinal pitch *L*. The porosity of both arrangements is expressed by $\epsilon = 1 - \frac{d^2}{L^2}$. The flow is considered as two-dimensional with the mean direction along the x-axis. The x- and y-directions are called longitudinal and transverse hereafter.

The maximum value of the Reynolds number based on the Darcy velocity and the size of the solid elements $Re_d$ is equal to 40 in the present study. It is well known that the two-dimensional flow around a single circular cylinder is stable in this low-range of $Re_d$ [21]. It is most likely that the case of an array of cylinders is more stable than an isolated cylinder. Dybbs and Edwards [22] cited by Kaviany [23] investigated the flow through packing of spheres and for complex arrangements of cylinders. They observed the onset of instabilities for $Re_p > 150$, where the Reynolds number $Re_p$ is based on the average pore velocity and an average characteristic length scale for the pores. This critical value corresponds to $Re_d$ much larger than 40 for the low-range of porosity. It was then assumed that instabilities did not occur for the present low values of $Re_d$ and that the flow around the solid elements was symmetrical with respect to the x-direction. Considering the above assumption, the computational domain as depicted in Fig. 1, therefore, consists of one wavelength $\lambda_x$ in the main flow direction along the x-axis and of only one half wavelength $\lambda_y$ in the y-direction ($\lambda_x = L$ or $2L$ for the aligned or the staggered arrangements, respectively). The flow was assumed to be laminar and incompressible with constant physical properties. Viscous dissipation was neglected.

With the above simplifications, the governing equations are:

$$ \text{Continuity equation } \nabla \cdot \vec{U} = 0, \qquad (1) $$

$$ \text{Momentum equation } \rho(\vec{U} \cdot \nabla \vec{U}) = -\nabla p + \mu \nabla^2 \vec{U}, \qquad (2) $$

$$ \text{Energy equation for the fluid phase } \rho c_p (\vec{U} \cdot \nabla T) = k_f \nabla^2 T. \qquad (3) $$

Most computations were performed with the condition of uniform volumetric heat source within the rods. In this case, the energy equation for the solid phase was

$$ \nabla^2 T + \frac{q_v}{k_s} = 0, \qquad (4) $$

where $q_v$ is the volumetric heat source.

The flow was supposed to be fully developed at the scale of the rods array, which allows assuming periodical conditions for the velocity field in the stream direction. The condition of uniform volumetric heat source also allows assuming periodical conditions for the heat transfer problem. The periodicity conditions are

$$ \vec{\Theta}(x, y) = \vec{\Theta}(x + \lambda_x, y) \qquad (5) $$

for any flux (velocity, heat flux). Since this condition could not be accounted for within the solid by the software used in this study, the control domain was defined with the inlet and outlet boundaries within the fluid (domain “a” in Fig. 1).

Due to periodicity, state variables as pressure and temperature can be written as the sum of a mean linear gradient and a periodic component

$$ \phi(x, y) = \frac{d\phi}{dx} x + \tilde{\phi}(x, y) \qquad (6) $$

$$ \text{with } \tilde{\phi}(x, y) = \tilde{\phi}(x + \lambda_x, y). \qquad (7) $$

The term $dp/dx$ is deduced from our computations while the term $dT/dx$ is determined by the energy conservation in the fluid domain. It was then assumed that $dT/dx$ in the solid and fluid phases is equal to $q_s d^2/2M$, where $M$ is the mass flow rate in the control domain per unit length in the spanwise direction.

Symmetry boundary conditions were assumed at the surfaces of the computational domain parallel to the main flow direction. They are written in form

<figure><img src="image_2.png" /><figcaption>Fig. 1. Computational domain.</figcaption></figure>
$$ \nabla \vec{\Theta} \cdot \vec{n} = 0, \quad \vec{\Theta} \cdot \vec{n} = 0, \quad \nabla \phi \cdot \vec{n} = 0, \qquad (8) $$

where $\vec{n}$ is the vector normal to the symmetry surface. The no-slip velocity condition and the continuity of temperature were assumed at all fluid-solid interfaces.

For the case of uniform temperature heating, $q_v$ was set to zero in Eq. (4) and a uniform temperature was prescribed at all fluid-solid interfaces. The thermal boundary condition in the fluid in x-direction was modified as explained later.

The flow may be defined by the Darcy velocity $u_D$ obtained by averaging the velocity over the total surface $L^2$ of an elementary cell. The Reynolds number $Re_d$ is based on $u_D$ and the size of the rods

$$ Re_d = \frac{u_D d}{\nu}. \qquad (9) $$

Additionally, the Darcian Reynolds number is defined by

$$ Re_D = \frac{u_D K^{0.5}}{\nu}, \qquad (10) $$

where $K$ is the permeability in the limit of creeping flow. Note that $K$ is a function of $d$ and $\epsilon$.

The Nusselt number is defined by

$$ Nu_d = \frac{\varphi d}{(T_s - T_f)k_f}, \qquad (11) $$

where $\varphi$, $T_s$ and $T_f$ are the heat flux density averaged over the fluid/solid interface, the temperature averaged over the solid surface and the fluid temperature averaged over the surface open to flow ($L^2 - d^2$), respectively.

## 2.2. Numerical scheme, meshing and numerical accuracy

The set of equations was solved with Fluent 6.1.22 by means of a second order upwind finite volume scheme. The SIMPLEC algorithm was used in order to improve convergence with regard to pressure-velocity coupling. The double precision solver was used and the convergence of results was assumed when the average pressure gradient and the average heat flux reached a constant value. The typical level of scaled residuals decreased below $10^{-8}$ for the continuity equation and $10^{-11}$ for the energy equation. An orthogonal grid generated by Gambit 2.1.2, was used with the size of the mesh cells equal to $L/200$. The grid then contained 400 × 100 cells in x, y directions, respectively, for the staggered arrangement. This mesh size was deduced from tests conducted for three different grids A, B, C, namely 200 × 50, 400 × 100, 800 × 200. The three meshes were tested with the porosity equal to 0.985 and the $Re_d$ number equal to 20. The pressure gradient and the $Nu$ number converged to their asymptotic values when the mesh size was decreased. The difference in the pressure gradient was equal to 0.5% between grids B and C and increased up to 0.7% between B and A. The difference in the Nusselt number was equal to 0.5%

between grids B and C and increased up to 2% between grids B and A.

## 2.3. Model validation

Computations were carried out for various cases of isothermal flows across banks of cylinders. The permeability was found in perfect agreement with the results of Martin et al. [11] for aligned cylinders of circular cross-section and with the results of Nakayama et al. [15] for aligned cylinders of square cross-section and zero yaw angle.

The numerical model was checked for laminar thermally and hydraulically fully-developed flow along a two-dimensional channel of height $e$ with symmetrical uniform temperature surfaces. For this thermal condition, the periodic condition for energy equation could not be kept in the form of Eqs. (1) and (2). According to the procedure recommended in [12], it was replaced by the condition of identical profiles of dimensionless temperature at inlet and outlet of the computation domain

$$ \theta(y) = \frac{T - T_w}{T_b - T_w}|_{x=0} = \frac{T - T_w}{T_b - T_w}|_{x=2L}. \qquad (12) $$

Iterative computations consisted in re-injecting at the channel inlet the temperature profile shape found at the outlet until convergence was obtained. The Nusselt number $Nu_{2e,b}$ was normalized with the hydraulic diameter $2e$, the wall and fluid bulk temperatures. The results are compared with those of Ash quoted by Shah and London [24] and those of Kuwahara et al. [12] in Fig. 2. The asymptotic trends, as given by Pahor and Strand [25] and Grosjean et al. [26] also quoted by Shah and London [24] are plotted in the same figure. The general trend of the variation $Nu_{2e,b} = f(Pe)$ is well recovered by the present computations. The agreement with the constant value of $Nu_{2e,b}$ (=7.54) observed for high $Pe$ is excellent. Axial conduction affects the convective heat transfer in the channel when the Peclet number is decreased. $Nu_{2e,b}$ departs from 7.54 for $Pe \approx 10-40$, depending on the authors and increases up

<figure><img src="image_3.png" /><figcaption>Fig. 2. Axial conduction effect on the fully-developed Nusselt number in a two-dimensional channel.</figcaption></figure>
to about 8.1 for the low Pe-range. The present results are slightly higher than those of the literature. However, the discrepancy is only 1.2% for Pe = 0.14. The largest difference with the results of Ash (+3.5%) is observed for Pe = 7. Our result is however in good agreement with the formulas of Pahor and Strand [25] and Grosjean et al. [26].

### 3. Results

#### 3.1. Hydrodynamics

For 2D cross-flow through an array of rods, the momentum equation is reduced to

$$0 = \frac{dp}{dx} + \frac{\mu}{K} u_D, \qquad (13)$$

when inertia is neglected. In the limit of creeping flow, the permeability does not depend on the flow velocity and may be related to $\epsilon$ by the Carman-Kozeny equation

$$K = \frac{d^2\epsilon^3}{C(1 - \epsilon)^2}. \qquad (14)$$

Bejan [27] reported that the Kozeny constant $C = 150$ for the cross-flow through the staggered arrangement of cylinders. The present computations relate the Darcy velocity $u_D$ to the pressure gradient for given geometric properties of the array. An “apparent” permeability $K_{app}$ and the corresponding value of $1/C$ are deduced from these results by using successively Eqs. (15) and (14)

$$K_{\text{app}} = \frac{\mu u_D_{\text{comp.}}}{-\frac{dp}{dx_{\text{comp.}}}}. \qquad (15)$$

The dimensionless parameter $1/C$ is plotted in Fig. 3 for the aligned and staggered arrangements. It is seen that $1/C$ is independent of $Re_D$, regardless the porosity, so that the linear Darcy’s law is satisfied for the low $Re_D$ number flows, as expected. For the low range of $Re_D$, the constant $C$ is equal to 130 for the staggered arrangement when the porosity $\epsilon = 0.44$. However, $1/C$ significantly decreases when $\epsilon$ is increased beyond $\epsilon = 0.8$. The similar trend is observed for the aligned arrangement, although the permeability is generally higher for this situation. However, the difference between the two sets of results only changes from about 2% for $\epsilon = 0.98$ up to 24% for $\epsilon = 0.44$ with respect to the staggered arrangement. Fig. 3 shows that inertia influences the flow when $Re_D$ is increased. With the choice of $K^{1/2}$ as length scale for the Reynolds number, this effect appears at a value of $Re_D$ ($\approx 1-10$) roughly independent of $\epsilon$. This limit value of $Re_D$ is slightly smaller for the staggered arrangement, as expected.

For the high $Re$ number flows, the Forchheimer modification of Darcy’s law has to be used

$$0 = \frac{dp}{dx} + \frac{\mu}{K} u_D + \frac{C_F}{K^{0.5}} \rho u_D^2, \qquad (16)$$

where $C_F$ is the Forchheimer coefficient. Fig. 4 presents the variations of the Forchheimer coefficient as a function of

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Normalized permeability as a function of $Re_D$ number for different porosities; (a) aligned arrangement and (b) staggered arrangement.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Variations of the Forchheimer coefficient as a function of porosity.</figcaption></figure>

porosity for the aligned and staggered arrangements. For both cases, the Forchheimer coefficient slightly increases when the porosity is decreased. Fig. 4 shows that the inertia effects are more pronounced for the staggered arrangement. This is obviously due to the many changes of direction of the stream in this configuration. Nakayama et al.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Dimensionless pressure gradient across the bank of rods. Comparison with Nakayama and Zukauskas results. Aligned arrangement. ε = 0.875.</figcaption></figure>

[15] investigated the effects of yaw on the pressure drop across a bank of cylinders of square cross-section and compared their results with those of Zukauskas [1]. Their results are plotted in Fig. 5 for the case of zero yaw angle and aligned arrangement together with the data of Zukauskas [1] presented in their paper. The model of the present study was used with the same value of porosity as that of Nakayama et al. [15]. Fig. 5 shows that the present results are in good agreement with the data obtained by these authors in the same conditions.

## 3.2. Uniform temperature heating

Computations were carried out with the constant wall temperature boundary condition and the staggered arrangement in order to compare the present results with previously published data. The calculation procedure described by Kuwahara et al. [12] was adopted in the current simulations. The computation domain was therefore shifted by L/2 in the x-direction as depicted in Fig. 1 (computational domain “b”). The periodic condition for the thermal field was taken into account as described in the previous section.

In the case of small porosity, the surface open to flow can be regarded as a series of narrow two-dimensional channels of width e (=L - d) following one another. The fluid domain may then be modelled by a succession of such channels, forming a “S-Z”-shaped channel about 3d in length in the computational domain (Fig. 1. Note that only one half of the longitudinal channels is included in the computational domain). This suggests to introduce the Nusselt number $Nu_{2e}$ and the Reynolds number $Re_{2e}$ based on the hydraulic diameter 2e of a two-dimensional plane channel and to compare the results with those of a single plane channel of dimensionless length d/e. $Re_{2e}$ is defined with the bulk velocity in a channel. The new dimensionless numbers are related to $Nu_d$ and $Re_d$ by

$$Nu_{2e} = Nu_d \frac{2e}{d}, \qquad (17)$$

$$Re_{2e} = \frac{2Re_d}{\sqrt{1 - \varepsilon}}. \qquad (18)$$

The results of Kuwahara et al. [12] were processed in this way for $Pr = 1$ and are drawn in Fig. 6 together with the present ones. $Nu_{2e}$ is normalized by $T_s - T_f$, like in Eq. (11) in order to do the comparison with [12].

Although the same trend is observed for the variation of $Nu_{2e}$ versus $Re_{2e}$, the current results are significantly higher (by about 40%) than those of Kuwahara et al. [12]. The small value of $\varepsilon$ used for this comparison (=0.36) corresponds to $d/e = 4$ for each channel located between the solid elements, which is not high enough to prevent an entrance effect at each channel inlet. Hence we carried out computations for a single channel of dimensionless length $d/e = 4$ in the two extreme cases when the flow is hydraulically and thermally fully-developed from the channel inlet (Case A) or oppositely when the velocity and temperature profiles are flat at the channel inlet (Case B). The results of these computations (Fig. 6) show that Case B gives rise to an enhanced axial conduction effect for low values of $Pe$. The data of the fully-developed situation are slightly higher than in Fig. 2 because the reference temperature is now the mean temperature $T_f$ instead of the bulk temperature $T_b$. The actual situation is obviously intermediate between Cases A and B because the distribution of velocity and temperature is more or less rearranged in the space between two consecutive channels of the control domain. The results of the computations carried out with the actual situation of a staggered arrangement are in good agreement with Case A for the lowest values of $Pe$ and follow the trend of Case B when $Pe > 100$. On the contrary, the results of Kuwahara et al. [12] are significantly lower than those of the two limiting cases. It is therefore concluded that they underestimate the Nusselt number for the low values of the porosity.

The computations were not possible for very low values of $Pe$ owing to results accuracy. For fully-developed flow

<figure><img src="image_8.png" /><figcaption>Fig. 6. Constant temperature heating. Comparison of the numerical results with two extreme cases of heat transfer in channel flow.</figcaption></figure>
through constant temperature parallel plates, it is well known (see for example, Bejan [27]) that the difference between the wall and fluid bulk temperatures decreases exponentially along the flow direction

$$ \sigma_0 = \frac{T_w - T_b(x)}{T_w - T_b|_{x=0}} = \exp(-4x^* Nu_{0-x}), \quad (19) $$

where $x^*$ is the dimensionless distance along the channel defined by $x^* = \frac{x}{2\epsilon} \frac{1}{Pe}$ and $Nu_{0-x} = 7.54$. The dimensionless temperature $\sigma_0$ at the exit of the computation domain may be estimated by Eq. (19) with $x = 3d$. The parameter $\sigma_0$ therefore decreases very rapidly and becomes vanishingly small when Pe is decreased. As a consequence, it is very difficult to carry out the computation for low values of Pe because the boundary condition Eq. (12) cannot be verified with sufficient accuracy. If we assume that the computations keep an accuracy at the level of $\alpha$, we can compute from Eq. (19) the value of Pe giving $\sigma_0 = \alpha$. We can then estimate the order of magnitude of the minimal Pe for which the computations are possible

$$ Pe_{\min} = \frac{-6 \times 7.54}{\ln \alpha} \frac{\sqrt{1 - \epsilon}}{1 - \sqrt{1 - \epsilon}}. \quad (20) $$

The result is not very sensitive to $\alpha$. The present numerical experiments suggested to take $\alpha \approx 10^{-5}$. For this value, Eq. (20) gives $Pe_{\min} \approx 15$ for $\epsilon = 0.36$.

This analysis shows that the computations are limited to a range of Pe larger than about 10, especially for low values of $\epsilon$. Hence it is very surprising that Kuwahara et al. [12] were able to obtain results for values of $Re_{2e}$ as low as $10^{-2}$, since Eq. (19) already gives $\sigma_0 = 2.6 \times 10^{-79}$ for $x = 3d$, $Re_{2e} = 1$, $\epsilon = 0.36$ and $Pr = 1$.

Fig. 7 shows the distribution of the local Nusselt number normalized with the fluid bulk temperature along the walls of the rods for the first half of the control domain. As explained above, the control domain inlet (point A) corresponds to the middle of a channel. For the transverse channel, $Nu_{2e,x}$ was obtained with the averaged heat flux on the opposite sides BD and CE of this channel. Fig. 7

<figure><img src="image_9.png" /><figcaption>Fig. 7. Distribution of local Nusselt number along the fluid/solid interface (constant wall temperature); $Pe = 14$, $Pr = 7$, $\epsilon = 0.36$.</figcaption></figure>

confirms that the entrance effect is weak in consecutive channels for a moderate value of $Pe$ (=14), since the Nusselt number rapidly decreases from the high values observed near the corners B and E to the constant 7.54 of the fully-developed regime.

### 3.3. Uniform volumetric source heating

The heating method by volumetric heat source in the solid elements avoids the computation difficulties encountered in the previous case. Instead of a continuously decreasing difference between the fluid and the wall temperature along the stream, the heat balance applied to the fluid contained in the control volume now implies a positive linear temperature gradient in the x-direction. The temperature field is then composed of this mean gradient superposed to local variations in the control domain. The periodicity condition implies the same longitudinal heat flux in the fluid and in the solid. However, the temperature is nearly uniform in each rod of the control domain if the conductivity of the solid is much higher than that of the fluid as in the present work ($k_s/k_f = 195$). The bank of rods then consists of successive nearly isothermal elements with increasing temperature from one element to the following one. As a result, the temperature jump between two adjoining elements gives rise to a significant heat transfer rate between the blocks and the fluid in the transverse channels. Contrary to the previous case, this region significantly contributes to the total heat exchange.

#### 3.3.1. Influence of the geometrical arrangement

For sake of clarity, the same definition of Reynolds number $Re_D$ as in the previous section was used to plot the heat transfer results. Fig. 8 shows that the $Nu_d$ variations are weak for low values of $Re_D$ and intensify when $Re_D$ is increased for both arrangements. This increase of $Nu_d$ is obviously due to an enhanced convective effect. It is observed that the Nusselt number $Nu_d$, similarly to the permeability with regard to inertia effects, is more influenced by this convective effect for the staggered arrangement (Fig. 8b).

#### 3.3.2. Influence of Prandtl number

The influence of Prandtl number is very important due to the wide range of this parameter used for flows through tube banks. The experimental and numerical results are commonly approximated by the power equation

$$ Nu_d = bRe_d^m Pr^n. \quad (21) $$

For low values of $Re_d$ the exponent of the Prandtl number varies around the level of $n = 1/3$ which is the theoretical value for a laminar boundary layer on a flat plate. Zukauskas [1] proposed $n = 0.36$ as sufficiently accurate for tube banks in various arrangements.

For low Reynolds number flows, the power formula Eq. (21) must be slightly modified by adding a constant value, which obviously accounts for the case of creeping flow
<figure><img src="image_10.png" /><figcaption>Fig. 8. Uniform volumetric source heating. Variation of Nusselt number as a function of $Re_D$ number; $Pr = 7$, (a) aligned arrangement and (b) staggered arrangement.</figcaption></figure>

$$Nu_d = a + bRe_d^m Pr^n. \qquad (22)$$

Computations were carried out for $Pr$ in the range 1–100. Fig. 9 shows the variations of $Nu_d$ as a function of $Re_d^m Pr^n$ for aligned and staggered arrangements, respectively. For both arrangements the exponent of $Re_d$ was found to fit well the results for $m = 0.5$. On the other hand, the best agreement with the computed data was found with different values of the exponent of $Pr$ for the aligned ($n = 0.2$) and the staggered ($n = 0.3$) arrangements. For the aligned arrangement, thermal boundary layers are forming only on the longitudinal walls of a block while the heat transfer on the front and rear walls is partially suppressed by an effect of “shading”, which could explain the lower exponent of the Prandtl number for the aligned arrangement. It is worth noting that the “a” and “b” constants of Eq. (22) depend on the porosity of the banks as shown in Fig. 10. In the limit of creeping flow, the difference between the values of $Nu_d$ for the two arrangements decays, since the constant “a” is nearly the same for both arrangements. The results were approximated by the least square method and are given by the following empirical expressions

<figure><img src="image_11.png" /><figcaption>Fig. 9. Effect of $Re_d$ and $Pr$ on Nusselt number; (a) aligned arrangement and (b) staggered arrangement.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Constants of Eq. (22) as a function of porosity.</figcaption></figure>

$$Nu_d = 3.02(1 - \varepsilon)^{0.278} \exp(2.54(1 - \varepsilon)) + ((1 - \varepsilon)0.44 + 0.092)Re_d^{0.5}Pr^{0.2}, \qquad (23)$$

$$Nu_d = 3.02(1 - \varepsilon)^{0.278} \exp(2.54(1 - \varepsilon)) + ((1 - \varepsilon)1.093 + 0.357)Re_d^{0.5}Pr^{0.3} \qquad (24)$$

for the aligned and staggered arrangements, respectively.
In order to assess the heat performance of banks it is interesting to relate the Nu<sub>d</sub> number to the hydraulic resistance. The Colburn analogy between the wall friction and the heat transfer coefficient may be written: $\frac{2St}{C_d} = Pr^{-2/3}$, where St = $\frac{Nu_d}{Re_d Pr}$ is the Stanton number and C<sub>d</sub> is the dimensionless resistance force exerting on the cylinders inside the control domain. Writing the permeability K<sub>app</sub> = $\frac{2L^2}{C_d Re_d}$, one obtains the expression

$$ \frac{Nu_d K_{app}}{L^2 Pr^{1/3}} = 1. \qquad (25) $$

When regarding the flow on a flat plate, the ratio $\frac{2St}{C_d}$ is related to the heat transfer rate obtainable per unit of pumping power and the Colburn analogy is valid both for laminar and turbulent boundary layer in absence of pressure forces. The flow through banks of tubes differs from the flow on a flat plate because the main part of hydraulic resistance is due to the pressure forces. This could explain why the current results are lower than the theoretical value 1 for the left-hand term of Eq. (25) (Fig. 11).

As can be seen in Fig. 11, the heat transfer performance curves exhibit a maximum. It occurs in the range of Re<sub>D</sub>, where the inertial force is negligible (Fig. 3) but the heat

transfer is already affected by the convection effect. For Pr = 7 the hydrodynamic entrance region formed at the inlet of each channel between neighbouring solid elements is considerably shorter than the thermal one, which explains the increase of thermal performance with Re<sub>D</sub> number. The maximum heat transfer performance is observed when the thermal boundary layers meet at the end of a channel as in the works of Bejan [27,28] on heat transfer optimization. For Pr = 1 the maximum is less pronounced than in the other case. Fig. 12 presents the comparison between aligned and staggered arrangements for Pr = 7. The heat transfer performance is slightly higher for the staggered arrangement. This effect is more pronounced for the small values of the porosity.

### 3.4. Influence of the thermal boundary condition

In most industrial situations, convection heat transfer in banks of tubes in cross-flow at moderate and high Reynolds numbers can be effectively modelled with a boundary condition of constant temperature surface. In fact, the temperature distribution in the cross-section of the tubes is almost uniform and the large heat capacity of the outer cooling fluid maintains its temperature at a constant level in the flow direction. This situation changes however when the outer flow is characterized by a low Reynolds number. In such case the low thermal capacity of the outer flow is responsible of a strong temperature gradient in the flow direction. In the conditions of the present study, the residence time of the fluid in a channel is long enough to ensure that the fluid reaches the temperature of the solid. In other words, the fluid reaches thermal equilibrium with the solid. This situation is illustrated by Fig. 13, which shows that the fluid rapidly reaches the solid temperature in each longitudinal channel. As a result the average temperature gradient (as in Eq. (6)) is the same for the solid and fluid phases. However, the local temperature gradient varies

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Thermal performance of tube banks in staggered arrangement; (a) Pr = 1 and (b) Pr = 7.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Thermal performance. Comparison between aligned and staggered arrangements for Pr = 7.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 13. Dimensionless temperature distribution $(T - T_{min})/(T_{max} - T_{min})$; $T_{min}$ and $T_{max}$ are the temperatures at inlet and outlet of the visualization window, respectively. Uniform volumetric source heating. $Pr = 7$, $Re_d = 0.5$, $\epsilon = 0.44$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15. Influence of the thermal boundary condition on heat transfer, $Re_d = 5$.</figcaption></figure>

substantially along the flow direction for rods of high thermal conductivity. As noted early, fluid regions of high temperature gradient then separate nearly constant temperature solid elements. On the other hand, writing the convection term of Eq. (3) in the form

$$||\vec{U}||\|\nabla T\| \cos(\vec{U}, \nabla T), \qquad (26)$$

one can conclude that the angle between the two vectors is close to 90° so that the convection term is almost negligible in a transverse channel despite the strong temperature gradient $||\nabla T||$. In the absence of convection, an elementary calculation shows that the local Nusselt number $Nu_{2e,x}$ should be equal to 4. This is confirmed by Fig. 14, which shows that $Nu_{2e,x}$ is close to 4 in the transverse channel BDEC. Two determinations of $Nu_{2e,x}$ were defined in the transverse channel by using the local wall temperature and heat flux on each side of the channel and the fluid bulk temperature. Like in the case of uniform temperature heating previously presented in Fig. 7, $Nu_{2e,x}$ rapidly approaches a constant value in the longitudinal channels. Since the boundary condition at the channel walls is not exactly constant temperature, the $Nu_{2e,x}$ asymptotic value is slightly above 7.54.

<figure><img src="image_18.png" /><figcaption>Fig. 14. Distribution of local Nusselt number along the fluid/solid interface (constant volumetric heat source).</figcaption></figure>

Fig. 15 presents the global Nusselt number $Nu_{2e}$ obtained for the two thermal conditions considered in this study and for $Re_d = 5$. The difference between the two cases is essentially due to the different contribution of the lateral walls BD and CE to the global heat transfer, as it is clearly shown in Figs. 7 and 14. For the case of constant temperature, the local $Nu_{2e,x}$ rapidly tends to the asymptotic value of 7.54, corresponding to fully-developed flow and heat transfer in a transverse channel with the walls BD and CE at the same temperature. In the case of constant volumetric heat flux, the heat transfer in a transverse channel is dominated by the conduction flux between the two walls BD and CE, as remarked before (see also Fig. 13) so that $Nu_{2e,x}$ tends to 4. This explains why the global $Nu_{2e}$ is slightly higher for uniform temperature heating than in the other case. This difference due to the thermal boundary condition disappears for the high range of Reynolds number where the outer fluid is always colder than the adjacent solid phase. On the contrary, low Re number flows are characterized by thermal equilibrium of the fluid and solid phases in the successive channels. As demonstrated by Fig. 13 for constant volumetric heat flux, the fluid temperature is thus equal to that of the adjacent solid at the outlet of a longitudinal channel and exhibits a continuous distribution intermediate between the temperatures of two successive solid elements within a transverse channel. It is concluded that the difference in heat transfer for the two boundary conditions is directly related to the condition of thermal equilibrium of the fluid and solid phases.

Since a compact solid matrix in a porous medium favours thermal equilibrium of the fluid and solid phases, it is obvious that the difference in $Nu_{2e}$ for the two boundary conditions is more pronounced for the small values of porosity. Fig. 15 also shows that the difference between the boundary conditions is higher when $Pr$ is decreased. Kim and Jang [29] proposed the criteria of local thermal equilibrium in the form
$$ \frac{\dot{M}c_p}{ha_{sf}} \ll 1, \qquad (27) $$

where h and $a_{sf}$ are the convective heat transfer coefficient and interfacial surface area per unit length in the spanwise direction, respectively. Eq. (27) can be equivalently written in the following dimensionless form

$$ \frac{Re_d Pr}{Nu_{2e}} \frac{1}{2\sqrt{1-\varepsilon}} \left( \frac{1}{\sqrt{1-\varepsilon}} - 1 \right) \ll 1, \qquad (28) $$

which shows that local thermal equilibrium is favoured by low values of $Re_d$, $Pr$, $\varepsilon$ and high values of $Nu_{2e}$. In fact, the present results (Fig. 15) show that the Nusselt number obtained with the two boundary conditions starts to depart when the left-hand side of Eq. (28) is about 1.5.

Finally, it may be remarked that the global Nusselt number $Nu_d$ (or $Nu_{2e}$) is normalized with the solid average temperature in the case of constant volumetric heat flux. It thus includes the thermal resistance due to internal conduction inside the solid elements. We can therefore expect that $Nu_d$ decreases with the solid conductivity. The decreased solid conductivity could also change the interfacial heat transfer coefficient but we have verified that it has a minor effect on the Nusselt number when compared with that of increased thermal resistance in the solid elements.

## 4. Conclusions

The present work is devoted to numerical simulations of the flow and associated heat transfer through banks of rods with different arrangements. The work focuses on the hydrodynamic resistance and more specifically on the heat transfer coefficient for low Reynolds number flows. This situation is of special interest because little data is available on this problem in the open literature. This study is also related to the problem of thermal equilibrium in a porous medium. The hydrodynamic resistance presented in term of permeability was found to agree well with published results. The present results show that the bank of rods can be modelled as a succession of narrow channels for estimating the Nusselt number in the case of tightly packed rods or equivalently for small values of porosity. Two different boundary conditions were investigated, namely constant wall temperature and constant volumetric heat source inside the solid elements. The convective heat transfer coefficient obtained with the former one was found significantly higher than previously published results [12]. However, the present results are in good agreement with the capillary model presented for a tightly packed bank of rods. A power equation was proposed in order to approximate the results obtained for the constant volumetric heat source condition. Finally the difference between the two thermal boundary conditions was discussed. It was shown that heat transfer in the array of rods was insensitive to this condition for the highest values of $Re_d$, $Pr$, $\varepsilon$ and the lowest values of $Nu_{2e}$. To our best knowledge, it was not possible to compare these numerical computations with published

well-documented experimental data. Experiments on this topic are therefore keenly encouraged.

## Acknowledgements

This research was supported by the CNRS and Rhône-Alpes Region. The fellowship of G. Gamrat was supported by the French Ministry of Education and Research, which is gratefully acknowledged.

## References

[1] A. Zukauskas, Convective heat transfer in cross flow, in: S. Kakac, R.K. Shah, W. Aung (Eds.), Handbook of Single-phase Convective Heat Transfer, Wiley, New York, 1987 (Chapter 6).

[2] D.A. Edwards, M. Shapiro, P. Bar-Yoseph, M. Shapira, The influence of Reynolds number upon the apparent permeability of spatially periodic arrays of cylinders, Phys. Fluids A 2 (1) (1990) 45–55.

[3] D. Nagelhout, M.S. Bhat, J.C. Heinrich, D.R. Poirier, Permeability for flow normal to a sparse array of fibres, Mater. Sci. Eng. A 91 (1995) 203–208.

[4] D.L. Koch, A.J.C. Ladd, Moderate Reynolds number flows through periodic and random arrays of aligned cylinders, J. Fluid Mech. 349 (1997) 31–66.

[5] M. Firdaouss, J.-L. Guermond, P. Le Quere, Nonlinear corrections to Darcy’s law at low Reynolds numbers, J. Fluid Mech. 343 (1997) 331–350.

[6] T.D. Papathanasiou, B. Markicevic, E.D. Dendy, A computational evaluation of the Ergun and Forchheimer equations for fibrous porous media, Phys. Fluids 13 (2001) 2795–2804.

[7] B. Markicevic, T.D. Papathanasiou, On the apparent permeability of regular arrays of nonuniform fibers, Phys. Fluids 14 (2002) 3347–3349.

[8] D.R. Graham, J.J.L. Higdon, Oscillatory forcing of flow through porous media. Part 1. Steady flow, J. Fluid Mech. 465 (2002) 213–235.

[9] W.H. Zhong, I.G. Currie, D.F. James, Creeping flow through a model fibrous porous medium, Exp. Fluids 40 (2006) 119–126.

[10] M. Fourar, G. Radilla, R. Lenormand, C. Moyne, On the non-linear behavior of a laminar single-phase flow through two and three-dimensional porous media, Adv. Water Resour. 27 (2004) 669–677.

[11] A.R. Martin, C. Saltiel, W. Shyy, Frictional losses and convective heat transfer in sparse, periodic cylinder arrays in cross flow, Int. J. Heat Mass Transfer 41 (15) (1998) 2383–2397.

[12] F. Kuwahara, M. Shirota, A. Nakayama, A numerical study of interfacial convective heat transfer coefficient in two-energy equation model for convection in porous media, Int. J. Heat Mass Transfer 44 (2001) 1153–1159.

[13] D. Ghosh Roychowdhury, S. Kumar Das, T. Sundararajan, Numerical simulation of laminar flow and heat transfer over banks of staggered cylinders, Int. J. Numer. Meth. Fluids 39 (2002) 23–40.

[14] V.K. Mandhani, R.P. Chhabra, V. Eswaran, Forced convection heat transfer in tube banks in cross flow, Chem. Eng. Sci. 57 (2002) 379–391.

[15] A. Nakayama, F. Kuwahara, T. Hayashi, Numerical modelling for three-dimensional heat and fluid flow through a bank of cylinders in yaw, J. Fluid Mech. 498 (2004) 139–159.

[16] M.B. Saito, M.J.S. de Lemos, Interfacial heat transfer coefficient for non-equilibrium convective transport in porous media, Int. Commun. Heat Mass Transfer 32 (2005) 666–676.

[17] G. Stanescu, A.J. Fowler, A. Bejan, The optimal spacing of cylinders in free-stream cross-flow forced convection, Int. J. Heat Mass Transfer 39 (2) (1996) 311–317.

[18] R.P. Taylor, H.W. Coleman, B.K. Hodge, Prediction of turbulent rough-wall skin friction using a discrete element approach, J. Fluids Eng. 107 (1985) 251–257.
[19] R.P. Taylor, H.W. Coleman, B.K. Hodge, Prediction of heat transfer in turbulent flow over rough surfaces, J. Heat Transfer 111 (1989) 568-572.

[20] R. Bavière, G. Gamrat, M. Favre-Marinet, S. Le Person, Modelling of laminar flows in rough-wall microchannels, J. Fluids Eng. 128 (2006) 734-741.

[21] D.J. Tritton, Physical Fluid Dynamics, first ed., Van Nostrand Reinhold Company, New York, 1977.

[22] A. Dybbs, R.V. Edwards, A new look at porous media fluid mechanics – Darcy to turbulent, in: J. Bear, Y.M. Corapcioglu (Eds.), Fundamentals of Transport Phenomena in Porous Media, Martinus Nijhoff Publishers, 1984, pp. 199–254.

[23] M. Kaviany, Principles of Heat Transfer in Porous Media, second ed., Springer-Verlag, New York, 1995.

[24] R.K. Shah, A.L. London, Laminar Flow Forced Convection in Ducts, Advanced Heat Transfer, Academic Press, New York, 1978 (Chapter 6).

[25] S. Pahor, J. Strand, A note on heat transfer in laminar flow through a gap, Appl. Sci. Res., Sect. A 10 (1961) 81-84.

[26] C.C. Grosjean, S. Pahor, J. Strand, Heat transfer in laminar flow through a gap, Appl. Sci. Res., Sect. A 11 (1963) 292-294.

[27] A. Bejan, Convection Heat Transfer, third ed., John Wiley & Sons, New York, 2004.

[28] A. Bejan, E. Sciubba, The optimal spacing of parallel plates cooled by forced convection, Int. J. Heat Mass Transfer 35 (1992) 3259-3264.

[29] S.J. Kim, S.P. Jang, Effects of the Darcy number, the Prandtl number, and the Reynolds number on local thermal non-equilibrium, Int. J. Heat Mass Transfer 45 (2002) 3885-3896.