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

Abstract 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.


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 6 Þ.
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][3][4][5][6][7][8][9][10]. These models are very attractive for numerical simulations, 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][12][13][14][15][16]. In the context of porous media, one of the issues is that of local thermal equilibrium 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 Prandtl number and on the thermal conditions on the solid surfaces. The range of parameters and boundary conditions found in [11][12][13][14][15][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 heating and for low Reynolds number flows. The motivations were twofold. Firstly, the thermal condition in heat exchangers is often neither uniform flux nor uniform temperature 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, 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 microheat 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 discreteelement 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-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 parallelepipeds 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 microchannel. 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 sources uniformly distributed inside the rods. Computations were also performed with the condition of uniform temperature heating for comparison with published data.

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 e ¼ 1 À d 2 L 2 . The flow is considered as two-dimensional with the mean direction along the xaxis. 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 k x in the main flow direction along the x-axis and of only one half wavelength k y in the y-direction (k 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: Energy equation for the fluid phase qc p ðŨ Á rT Þ ¼ k f r 2 T : ð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 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 arẽ Hðx; yÞ ¼Hðx þ k x ; yÞ ð 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 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 v d 2 =2 _ Mc p L , 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 whereñ 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 Additionally, the Darcian Reynolds number is defined by where K is the permeability in the limit of creeping flow. Note that K is a function of d and e. The Nusselt number is defined by where u; 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.

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.

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 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 % 10-40, depending on the authors and increases up 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].

Hydrodynamics
For 2D cross-flow through an array of rods, the momentum equation is reduced to when inertia is neglected. In the limit of creeping flow, the permeability does not depend on the flow velocity and may be related to e by the Carman-Kozeny equation 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 app ¼ lu Dcomp: À dp dxcomp: : 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 e ¼ 0:44. However, 1=C significantly decreases when e is increased beyond e ¼ 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 e ¼ 0:98 up to 24% for e ¼ 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 (%1-10) roughly independent of e. 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 where C F is the Forchheimer coefficient. Fig. 4 presents the variations of the Forchheimer coefficient as a function of 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.  [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.

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 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 e 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  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 where x * is the dimensionless distance along the channel defined by x Ã ¼ x 2e 1 Pe and Nu 0Àx ¼ 7:54. The dimensionless temperature r 0 at the exit of the computation domain may be estimated by Eq. (19) with x ¼ 3d. The parameter r 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 a, we can compute from Eq. (19) the value of Pe giving r 0 ¼ a. We can then estimate the order of magnitude of the minimal Pe for which the computations are possible The result is not very sensitive to a. The present numerical experiments suggested to take a % 10 À5 . For this value, Eq. (20) gives Pe min % 15 for e ¼ 0:36. This analysis shows that the computations are limited to a range of Pe larger than about 10, especially for low values of e. 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 r 0 ¼ 2:6 Â 10 À79 for x ¼ 3d, Re 2e ¼ 1, e ¼ 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 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.

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.

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).

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 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 Computations were carried out for Pr in the range 1-100. Fig. 9 shows the variations of Nu d as a function of Re m d 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 for the aligned and staggered arrangements, respectively.   In order to assess the heat performance of banks it is interesting to relate the Nu d number to the hydraulic resistance. The Colburn analogy between the wall friction and the heat transfer coefficient may be written: 2St Re d Pr is the Stanton number and C d is the dimensionless resistance force exerting on the cylinders inside the control domain. Writing the permeability K app ¼ 2L 2 C d Re d , one obtains the expression When regarding the flow on a flat plate, the ratio 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 D , 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 D 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.

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  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 krT k. 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. 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 fullydeveloped 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 which shows that local thermal equilibrium is favoured by low values of Re d , Pr, e 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.

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, e 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.