Max Jakob Award Paper

Conjugate Heat Transfer Within a Heterogeneous Hierarchical Structure PUBLIC ACCESS

[+] Author and Article Information
Ivan Catton

Department of Mechanical and Aerospace Engineering, School of Engineering and Applied Science,  University of California, Los Angeles, Los Angeles, CA, 90095

J. Heat Transfer 133(10), 103001 (Aug 11, 2011) (16 pages) doi:10.1115/1.4003576 History: Received October 12, 2010; Revised January 28, 2011; Published August 11, 2011; Online August 11, 2011

Optimization of heat exchangers (HE), compact heat exchangers (CHE) and microheat exchangers, by design of their basic structures is the focus of this work. Consistant models are developed to describe transport phenomena in a porous medium that take into account the scales and other characteristics of the medium morphology. Equation sets allowing for turbulence and two temperature or two concentration diffusion are obtained for nonisotropic porous media with interface exchange. The equations differ from known equations and were developed using a rigorous averaging technique, hierarchical modeling methodology, and fully turbulent models with Reynolds stresses and fluxes in the space of every pore. The transport equations are shown to have additional integral and differential terms. The description of the structural morphology determines the importance of these terms and the range of application of the closure schemes. A natural way to transfer from transport equations in a porous media with integral terms to differential equations with coefficients that could be experimentally or numerically evaluated and determined is described. The relationship between computational fluid dynamics, experiment and closure needed for the volume averaged equations is discussed. Mathematical models for modeling momentum and heat transport based on well established averaging theorems are developed. Use of a “porous media” length scale is shown to be very beneficial in collapsing complex data onto a single curve yielding simple heat transfer and friction factor correlations. The general transport equations developed for a single phase fluid in a heat exchange medium have many more integral and differential terms than the homogenized or classical continuum mechanics equations. Once these terms are dealt with by closure, the resulting equation set is relatively simple and their solution is obtained using simple numerical methods quickly enough for multiple parameter optimization using design of experiment or genetic algorithms. Current efforts to significantly improve the performance of an HE for electronic cooling, a two temperature problem, and of a finned tube heat exchanger, a three temperature problem, are described.

Many heat transfer devices can be called heterogeneous and hierarchical. We treat such problems everyday by assigning a heat transfer coefficient and a friction factor to a simple statement and are on our way to obtaining a solution. For many problems, this is sufficient, but did you ever wonder how this might relate back to the pointwise correct Navier–Stokes and energy equations for the fluid and solid side of the problem. An example of such a structure is the simple heat sink.

The simple heat sink requires many parameters to describe its geometry making optimization difficult. The simple equations are the only answer if we want to find the optimum configuration to this conjugate problem. Computational fluid dynamics (CFD) or experiment by itself is out of the question, and simple equations are needed but they need to be made more rigorous. This can be accomplished by incorporating more rigors into our choice of the coefficients that couple the fluid to the solid phase: the friction factor and the heat transfer coefficient. It is proposed that volume averaging theory (VAT) be used to develop the simple equations allowing clear rigorous statements to be made that define how the friction factor and heat transfer coefficient are to be determined.

For optimization to be accomplished, one needs to separate conduction and convection at its natural interface. This means friction and heat transfer at the interface need to be treated in a way that allows each side to be independently varied. If you increase the flow rate, the heat transfer coefficient will increase and you need to take advantage of this by changing the geometry to modify the conduction side of the problem. Of course you try to do this while keeping the friction factor as small as possible. All of this is done within a complex heterogeneous hierarchical structure.

Several issues arise when we set out to solve such a problem. One of the most interesting issues is how you compare CFD calculations to measured values of the two parameters that are needed. Another is how to directly use experimental data. The simple equations contain coefficients that are essentially averages of some kind and we are usually unsure of what they are. Measured friction factors are fairly straightforward to interpret although calculated values are not. Measured heat transfer coefficients are usually based on some kind of upper level assumptions and often not what we need, and calculated values can be difficult to interpret correctly. Using VAT to derive the simple equations leads us to a rigorous definition of the needed coefficients whether obtained from experiment or through the use of CFD.

Determination of flow-variables and scalar transport in a heterogeneous (or porous) media is difficult even when subject to simplifications allowing the specification of medium periodicity or regularity. Linear or linearized models fail to intrinsically account for transport phenomena, requiring dynamic coefficient models to correct for short-comings in the governing models. Allowing inhomogeneities to adopt random or stochastic character further confounds the already daunting task of properly identifying pertinent transport mechanisms and predicting transport phenomena.

Some aspects of the development of transport phenomena in porous media are now well understood and have seen substantial progress in thermal physics and fluid mechanics sciences, particularly in porous media transport phenomena. The basis for this progress is VAT first proposed in the sixties by Anderson and Jackson [1], Whitaker [2], and others. Many of the important details and examples of its application are found in books by Kheifets and Neimark [3] and most recently by Whittaker [4].

In most physically realistic cases, highly complex integral–differential equations result. The largest challenge is the insufficient development of closure theory, especially for integral–differential equations. The ability to accurately evaluate various kinds of medium morphology irregularities results from the modelling methodology once a porous medium morphology is assigned. Further, when attempting to describe transport processes in a heterogeneous media, the correct form of the governing equations remains an area where inattention to procedure by some researchers has led to significantly different equations for the same media (see Travkin and Catton [5-6]).

The VAT based approach has the following desirable features:
  • Effects of interfaces and boundaries can be included in the modeling.
  • The effect of morphology of the different phases (solid, liquid, and vapor) is incorporated. The morphology description is directly incorporated into the integrodifferential field equations.
  • Separate and combined fields and their interactions are described exactly through the integral terms appearing in the field equations. No assumptions about effective transport coefficients are required.
  • Effective coefficients are rigorously derived whereas many present results are often simply wrong.
  • Deliberate design and optimization of materials using hierarchical physical descriptions based on the VAT governing equations can be used to connect properties and morphological characteristics to heat exchanger features.

The method of volume averaging begins by associating an averaging volume ΔΩ with every point in space. The averaging volume, the representative elementary volume (REV), is illustrated in Fig. 1 and is represented in terms of the volumes of the individual phases by Display Formula

The volumes of the fluid and solid phases can depend on position and, in addition, the volume of the fluid phase can depend on time. Average values are defined in terms of these volumes and these averaged values are associated with the centroid of ΔΩ.

In the development of volume averaged transport equations, a superficial average and an intrinsic average are defined. The superficial average of f is given by Display Formula

and the intrinsic average of f is defined by Display Formula
The intrinsic and superficial averages are related by Display Formula
where m is the porosity in the averaging volume ΔΩ defined by Display Formula
The random porosity function in terms of the average value of m(x) in the REV and its fluctuations m̂(x) in various directions is Display Formula
The different types of averaging over the REV function f are defined by the following averaging operators: The average of the function f over the REV is given by Display Formula
the phase average of the function f in the fluid Display Formula
the phase average of the function f in the solid phase Display Formula
and the interphase average of each component is Display Formula
Display Formula
When the interface is fixed in space, averaged functions for the first and second phase within the REV and over the entire REV fulfill all four of the Reynolds conditions as well as the following four consequences Display Formula
The differential operator average Display Formula
results from using the Δ averaging theorem, see Whitaker [4]. The fourth condition implies an unchanging porous medium morphology. All the above are well known from earlier work by Slattery [7], Whittaker [8], and Gray et al.  [9]. At the same time, ff and fs fulfill neither the third Reynolds conditions nor all the consequences of the other Reynolds conditions.

Most of the following analysis is based on averaging techniques developed by Whitaker [4] who focused on solving linear diffusion problems and by Travkin et al.  [10] who focused on solving nonlinear turbulent diffusion problems. The divergent and nondivergent forms of the averaged convective term in the diffusion equation are Display Formula

where T̃=T-T̂.

Applying the VAT rules and definitions to the heat flux (or mass flux) gradient term with spatially dependent heat conductivity kT yields Display Formula

Consistent application of these averaging rules and definitions to the conservation equations will allow the VAT governing equations to be systematically and rigorously derived.

VAT Based Equations Describing Thermal Physics and Fluid Mechanics in Heterogeneous Media.

VAT based thermal physics and fluid mechanics governing equations in heterogeneous porous media developed from the Navier–Stokes equation and the thermal energy equation are the starting point and the basis for studying flow and heat transfer in porous media. The turbulent transport considered here is flow in a highly porous layer combined with a two energy equation model.

For flow in porous media, the instantaneous velocity for the turbulent regime is represented by Display Formula

where Display Formula
k signifies turbulence induced components that are independent of inhomogeneities in spatial dimensions and properties resulting from the multitude of porous medium channels (pores), and r denotes the fluctuation contributions due to the porous medium inhomogeneity. The fluctuation velocity ur(x,z) due to porous medium inhomogeneity is ur"(x,z)=u¯̂(x,z)

The temperature is represented by Display Formula


Development of the Continuity Equation.

For incompressible flow, the continuity equation is Display Formula

Using Eq. 12, the averaged continuity equation in the fluid phase for flow in a porous media is written as Display Formula

Development of the Momentum Equation.

Starting from the momentum transport equation for incompressible flow Display Formula

and using Eqs. 12,13,14, the averaged momentum equation is written as Display Formula
where Fi incorporates the impact of microroughness and augmentation of the previous level of the simulation hierarchy.

For many general flow situations, the Reynolds stress can be approximated by Display Formula

For one-dimensional parallel flow, under steady-state conditions and no flow penetration through Sw, Eq. 21 has the form Display Formula
An equation for one dimensional steady fully developed turbulent flow in a porous layer with regular packed bed structure characteristics and an impermeable interphase surface simplifies to Display Formula

Development of the Turbulent Kinetic Energy Equation.

Following Rodi [11], the turbulent kinetic is written as Display Formula

where CD,Cμ, and σb are empirical coefficients. CDCμ0.08 and σb=1 appear to be reasonable values for the empirical constants. Assuming that most of the mean motion kinetic energy lost due to interaction of the flow with the porous medium solid obstacles translates into increasing the turbulent fluctuation energy (see Ref. [11]) leads to the conclusion that Display Formula
It follows that the equation for the mean turbulent fluctuation energy b(z) can be written in the form Display Formula
where the mean eddy viscosity ν̃T can be determined using Prandtl’s mixing length theory. By choosing b as the velocity scale, ν̃T can be expressed as Display Formula
where L is the mixing length scale function defined by the assumed porous medium structure and Display Formula
is the turbulent kinetic energy. Equations 28,29 were used in the modeling of a regular porous medium filled channel transport by Travkin and Catton [5,12] and Gratton et al.  [13-14].

Development of the Energy Equations.

To derive the volume averaged energy equation, the energy transport equation with turbulent transport Display Formula

is the starting point. Using Eqs. 13,14,15, the averaged energy equation in fluid phase becomes Display Formula
The term ST represents a heat source.

By analogy with the concept of eddy viscosity, the turbulent fluctuation terms are proportional to the mean temperature gradients. For two dimensional heat transfer in a channel and no flow penetration through Sw, the above equation simplifies to Display Formula

where k̃T is the turbulent heat conductivity.

The solid phase energy equation is Display Formula

where ks is the effective heat conductivity in the solid phase. For two dimensional heat transfer in the solid phase, the volume averaged solid phase thermal energy equation is written as Display Formula
where ds1=-ds.

Closure theories for transport equations in heterogeneous media have been the primary measure of advancement and for measuring success in research on transport in porous media. It is believed that the only way to achieve substantial gains is to maintain the connection between porous medium morphology and the rigorous formulation of mathematical equations for transport. There are only two well known types of porous media morphologies for which researchers have had major successes, namely, straight parallel pores and equal size spherical inclusions. But even for these morphologies, not enough evidence is available to state that the closure problems for them “are fully closed.”

Before moving to development of the closure relationships, it is interesting to note that using a particular length scale leads to a parameter that is very beneficial when scaling heat transfer and friction factor results. It was shown by Travkin and Catton [15] that media of globular morphologies can be described in terms of Sw , m, and dp and can generally be considered to be spherical particles with Display Formula

This expression has the same dependency on equivalent pore diameter as found for a one diameter capillary morphology leading naturally to Display Formula
This observation leads to defining a simple “universal” porous media scale Display Formula
that meets the needs of both morphologies: capillary and globular. This was also recognized by Whitaker [16] when he used a very similar (differing by a constant) length scale when correlating heat transfer for a wide variety of morphologies. Here, we will show results for correlation of friction factor and heat transfer coefficient for a fin tube heat exchanger.

The closure terms in the VAT momentum equations are (see Eq. 26) Display Formula

The first and second terms are the laminar and turbulent flow contributions to skin friction, the third term represents the spatial flow oscillations induced by the porous media morphology, and the fourth term is the form drag.

A natural way to close the integral terms in transport equations is to integrate over the interphase surface or of some other outlined areas of this surface. From a physical point of view, the integration terms in Eq. 39 represent momentum loss due to the friction resistance and pressure (form) drag over all interphase surfaces. It should be noted that this is relatively easy with modern commercial CFD packages.

The pressure drag resistance integral terms appearing in Eq. 39 are closed in a manner similar to that for a pressure resistance coefficient over an obstacle or hydrofoil, where according to its definition, the form drag is defined as Display Formula

where Swp is the cross flow projected area. Now, using this definition, the pressure drag resistance integral term can be closed with Display Formula
Similarly, the skin friction resistance integral terms appearing in Eq. 25 are closed using the definition of skin friction coefficient defined by Display Formula
Travkin and Catton [12] showed that using the proper length scale, like the one represented by Eq. 38, enables one to write the friction factor of porous media in the following form Display Formula
where ff is the Fanning friction factor and the Reynolds number is defined as Display Formula
with A = 100/3 and B = 7/12 for the Ergun equation for packed bed porous media. These constant values were shown to be very close to the correlation developed by Klevon and Matros [17] with A = 36.3 and B = 0.45 for spherical particles and A = 37.6 and B = 0.585 for cylindrical particles. Travkin and Catton [12] were able to show that by choosing the above form, it allows one to collapse true capillary flow and flow in a bed of spheres. This is a significant accomplishment, since it spans the physical description from globular to capillary geometry with a single length scale.

The friction factor presented by Eq. 43 has two distinctive terms, each corresponding to different pressure drop mechanisms. One is attributed to the viscous flow, which is linearly proportional to velocity, and the other term is due to convection of the fluid momentum that is proportional to the velocity squared. Furthermore, Travkin and Catton [12] stated that friction factor is related to a closure for VAT momentum equations and they showed that Display Formula

It was shown that the closure equation (Eq. 39) is an exact definition of friction factor. Note that for fully developed flow, Eq. 43 is more strictly defined as Display Formula
where A and B are constants corresponding to the different types of morphology of porous media.

The fluid energy equation can be closed in various ways and in general will depend on how many of the integrals appearing in the VAT equation one uses and lumps into a single transport coefficient, see Travkin and Catton [15]. The nature of the equation shows that the energy transferred from the surface is integrated over an area and then divided by a REV volume; therefore, the heat transfer coefficient is defined in terms of porous media morphology, usually described by specific surface and porosity. The simplest most intuitive and most often used variation includes only the intraphase heat transfer and is defined as Display Formula

By the same token, one could include more terms, such as spatial velocity and temperature fluctuations, yielding a new form defined by Display Formula
Lastly, a third form could be defined that uses all of the fluid energy equation integral terms to yield the most complete form Display Formula
The number and the nature of the closure terms in the macroscale equations were rigorously derived from microscale governing equations and they are clearly defined. The way one defines a heat transfer coefficient can vary as shown by the three definitions of the heat transfer coefficient and therefore, one has to be careful on which variation is applicable to the problem that one tries to solve. When in doubt, one should use the complete form given by Eq. 47. In most engineered devices, the geometry is regular and a well chosen REV will lead to only the first type being needed.

The first and third terms in Eq. 49 can be estimated using scaling arguments. It is assumed that the length scale within the REV is l and that the whole REV length scale is L. With these assumptions, the first term can be estimated as Display Formula

while the third term is scaled as Display Formula
Since the macrolength scale is much larger than the microscale, L>l, the third term is much smaller than the first and can therefore be ignored yielding Display Formula
For conditions of fully developed flow, the spatial flow oscillations within the REV are constant or nearly constant, generally depending on the boundary conditions: either constant flux or isothermal. Thus, the product of flow oscillations is also constant and the second term in Eq. 48 or Eq. 49 is essentially zero, yielding Display Formula
This is an important result as it greatly simplifies the scope of the problem. The use of scaling arguments showed that for fully developed flow, small scale heat transfer is the dominating term in the closure formulation and it is much greater than the heat transfer on the larger macroscale. Furthermore, for spatially periodic media morphology, the velocity and temperatures fluctuation term is equal to zero, which further states that the main mechanism of intraphase transport is the small scale intraphase heat transfer defined by the energy equation closure (Eq. 51).

How the REV is selected is very important. The selection for a finned tube exchanger, Fig. 2, is seen to repeat in both the cross-stream and flow directions. The fluctuation terms in Eq. 50 were evaluated and found to be zero. This, not unexpected, result simplifies the CFD task.

An example of how closure is obtained from experimental data is given by Zhou et al.  [18]. The porosity and specific surface are geometrically defined terms. Figure 2 depicts the fin geometry studied in their work. For a plate finned tube heat exchanger, the porosity for the air side is Display Formula

m 1=1-δfFp-πDc2(Fp-δf)4PlPtFp
and for the water side is Display Formula
The specific surface area for the air side is given by Display Formula
and for the water side Display Formula
Wang et al.  [19] developed precise correlations for the friction factor for the air side of plain finned tube heat exchangers. They were scaled with the fin collar outside diameter and the maximum velocity. The friction factor results are shown on the left of Fig. 3 and the rescaled results are shown on the right clearly demonstrating the value of the VAT based length scale. The experimental data of Wang et al.  [19] were rescaled using the length scale given by Eq. 38. The results were surprisingly good, collapsing all the data onto a single curve. This, in turn, with the help of JMP 8, an available statistical analysis tool, enabled us to develop a simple correlation for the friction factor, using the form suggested by Eq. 43Display Formula
A comparison of some values of A and B with other morphologies is shown in Table 1. The type of correlation given by Eq. 43 is very powerful, allowing one to estimate the result when no data are available.

The VAT closure for the fluid energy equation is found in terms of the heat transfer coefficient h, (see Eq. 49). Experimental data from Xie et al.  [20] for N = 7 and by Tang et al.  [21] for N = 9 and N = 12 were plotted on a single graph for each of the length scales, Dc and dh , and are shown in Fig. 4. It can be seen that all the data collapse very close to a single curve when scaled using the VAT based length scale given by Eq. 38. Using JMP 8, a correlation of the rescaled data was found and is Display Formula

It is worth noting that the correlation is similar to the Dittus–Boelter correlation (Nu=0.023Redh0.8Pr1/3) for forced convection and, in some sense, a testament to the power of the hydraulic diameter approach to scaling.

Closure Using CFD Results.

CFD can be used to obtain detailed solutions of flow through a representative element of some morphology, and these results can be used to evaluate the closure terms needed for a fast running VAT based model by integrating the VAT based closure formula for drag coefficient and heat transfer coefficient over the surface of the solid parts of the geometry.

Numerical calculations were performed by Vadnjal [22] at different Reynolds numbers for a simple body center cubic packing arrangement of spheres, see Fig. 5. This configuration was chosen because it has been well studied, it is not too complex, and, although not always in agreement, a great deal of data exists.

Once the CFD results were obtained, the integrations needed to evaluate the terms in Eq. 39 were used to obtain drag resistance coefficients (friction factor) and compared to the well established results of Ergun [23], Morcom [24], Burke and Plumer [25], and a correlation developed by Ergun [23]. The calculated value included all three terms given by Eq. 39 including the fluctuating terms. The experimental data are representative of porous media morphology composed of different size spheres, sand, and pulverized coke.

The comparison is shown in Fig. 6 and the results show the validity of the VAT calculations for laminar flow regimes where the Reynolds number, Reh , ranged from 0.4 to 124. Agreement is good for Reynolds numbers up to approximately 50, whereas for higher values, the calculation yielded slightly higher values of the friction factor. It is interesting to note that the curve starts to change at the particular Reynolds number where the flow is transitioning from a laminar to a turbulent regime. The purpose of this work was to model laminar flow and to avoid CFD turbulence modeling and the approximations that are associated with it.

Closure for the momentum equation shows that the drag resistance coefficient has two distinctive contributions: one is due to pressure and the other is due to viscous forces. The relative contributions to the total drag resistance coefficient are shown in Fig. 7, where it can be seen that higher Reynolds numbers lead to a pressure contribution that is much greater than the viscous contribution. This demonstrates the need to describe the morphology in terms of both the total specific area for calculation of viscous forces and the specific area normal to the flow direction to determine the form drag. Further, the form of the drag coefficient lends itself to a single porous media type code for a wide variety of internal structures.

To obtain energy equation closure, the heat transfer coefficient is obtained by resolving the fluid flow and the heat transfer on the lower microscale. The use of scaling arguments showed that for fully developed flow, small scale heat transfer is the dominating term in the closure formulation and it is much greater than the heat transfer on the larger macroscale. Furthermore, for spatially periodic media morphology, the velocity and temperatures fluctuation term is equal to zero, which further states that the main mechanism of intraphase transport is the small scale intraphase heat transfer defined by the energy equation closure (Eq. 53).

The numerical results for heat transfer closure are shown in Fig. 8 and are compared to the various experimental results presented by the correlation given by Whitaker [16]. The result obtained by VAT closure included only the first term and show fair agreement with the data; therefore, one can confirm that the heat transfer coefficient when scaled correctly can be computed numerically with satisfactory results.

The local and average Nusselt number shown in Fig. 8 differed from the average as expected and therefore one has to pay special attention when using heat transfer correlations and especially when one uses the heat transfer coefficient with the VAT governing equations to solve large scale problems. When compared to the available experimental data, one can conclude that calculated values, either local or average heat transfer coefficients, are within the experimental bounds. However, it is not so clear how the experimental data were processed, whether the values are local or averaged and therefore the differences are obvious, and even a couple of orders apart when Reynolds number is small.

Having obtained solutions at the microscale, the entrance effects in a porous media were addressed. The interest here was to determine how averaged temperatures and heat transfer coefficients vary as the flow is developing from its initial uniform entrance condition to its thermally fully developed state. The microscale solution of a five REV long morphology was averaged over a single REV. Average temperatures for fluid and solid phase were computed for a flow with Reh  = 31.2 at different locations and the results are shown on Fig. 9. As expected, the average temperatures of fluid and solid are increasing from the inlet to the exit. The largest temperature difference is at the inlet while the smallest is at the exit. Intuitively, one would to expect the heat transfer coefficient to be largest at the entrance and that it would gradually asymptote to some value until thermally fully developed flow is reached.

By inspecting the heat transfer coefficient definition, Eq. 53, one can see that it has two parts: one in the numerator and the one in the denominator, each representing the physics of interfacial heat transfer. The integral term represents the heat transfer across the interface while the temperature difference represents its potential. Figure 1 shows their values and their development from the inlet to the exit and it is interesting to see that their slope is almost identical. With both numerator and denominator computed, the local heat transfer coefficient can be computed and the results are plotted on Fig. 1. It is surprising to see that the local coefficient is almost constant throughout the domain and it is not largest at the inlet where the temperature difference between solid and fluid is the largest. Those results were achieved at fairly low Reynolds number, and to make this result more convincing, the same steps were repeated at a higher Reynolds number, Reh  = 440.8, and are shown on Figs.  1213. Again, it is shown that the local heat transfer coefficient is essentially constant from the inlet to the exit.

In general, at the entrance, the point local and average heat transfer coefficients will both be some large number and they both asymptote to corresponding numbers as the Graetz number goes to a large value. For the microscale solution obtained for the whole domain, a local point value heat transfer coefficient on a solid–fluid interface was found to be largest at the entrance and gradually decreases as the flow progresses further downstream. However, the local VAT average heat transfer is essentially independent of entrance region. This remarkable result is very helpful when developing a simple model based on VAT.

To explore how various heat sources and boundary conditions affect the heat transfer coefficient, different heating methods were chosen and a comparison of local heat transfer coefficients was made. The choices are constant temperature, constant heat flux, and constant volumetric heat flux. The magnitude for each boundary condition is irrelevant since the heat transfer coefficient does not depend on how much heat is applied to the fluid, or in other words, by varying energy sources or sinks, the temperature difference will adjust accordingly, and therefore a boundary condition with a constant heat flux of 20 W/m2 (BC:A), a constant temperature at T = 70 °C (BC:B), and a constant volumetric heat of Q·=1.5kW/m3 (BC:C) was chosen. The values of the local heat transfer coefficients with different types of heating are shown on Fig. 1. The local Nusselt number, Nu = hdh /k, was essentially a function only of the magnitude of the Reynolds and Prandtl numbers and is not dependent on the method of heating or boundary conditions.

A finned tube heat exchanger geometry was investigated by Zhou et al.  [18] to make sure that these simplifying aspects of our search for closure could be extended to other geometries. This was a two step process where the CFD results were first compared to experimental data (see Fig. 1), and then the VAT closure terms were calculated.

Local values of the heat transfer coefficient were calculated using Eq. 53. The calculated values were again found to be independent of location and only a function of Reynolds number (see Fig. 1). This is because the REV repeats itself in all directions leading to zero value for the fluctuation terms.

Local values are important because they are the only values that have a physical meaning when describing transport phenomena with VAT macroscale equations. The average values represent values that are valid over some length at a location where the flow is fully developed and they represent the total intraphase heat transfer. The VAT local values, however, represent a local heat transfer coefficient that is based on a local average fluid temperature and a local average solid temperature. To be consistent with the VAT theory, only the local values of heat transfer closure can be used in calculations to accurately predict total heat transfer in a porous media, although average values are more important for comparison with experimental data.

An important aspect of what we have done is to show how CFD can be used to close the VAT based equations. There are a number of terms that contribute to the closure, see Eqs. 40,41,42,43. They are calculated using results from the CFD solution to the problem. The local Nusselt number was found to depend only on the Reynolds number and Prandtl number and to be independent of the type of heating. The friction factor was shown to be almost independent of the geometry with a universal closure for pin fins in both inline and staggered configuration and for various packed bed porous media. A table showing closure values of A and B for various physical configurations with their range of validity is shown in Table 1.

Most importantly, one can use heat transfer and friction factor data from any experiment that reports the values for fully developed flow providing the local average velocity and local values of the fluid and wall temperatures are used in their development.

The governing equations used in the analysis are the volume averaged momentum and energy equations describing transport phenomena within and between the fluid and solid phase (Travkin and Catton [15]). The averaged turbulent momentum equation is Display Formula

The terms on the right hand side represent pressure drop, internal drag (two terms), viscous dissipation, and Reynolds stresses. Internal drag is a sum of friction and form drag, due to the solid fluid interfaces inside of the channel. The turbulent kinetic energy of a porous media flow from Gratton et al.  [14] Display Formula
is used to close the momentum equation. The relationship between the eddy viscosity and the turbulent kinetic energy is Display Formula
where l(z) is a turbulent scale function, defined by an assigned porous medium structure from Travkin and Catton [5,26], is the turbulent exchange coefficient, and b(z) is the mean turbulent fluctuation energy function.

The energy equation for the fluid phase is given by Display Formula

and for the solid phase by Display Formula
Both the momentum and the energy equations need closure: the friction, ff, for the momentum equation, and the internal heat transfer coefficient, h, for the energy equations. The needed closure is obtained experimentally from available data for a specific morphology. The remaining parameters appearing in the VAT conservation equations, 〈m〉 and Sw , are the porosity and the specific surface area. For engineered devices, these are easily calculated. The set of equations, Eqs. 60,61,62,63,64, is valid for all geometries needing only two parameters to describe the morphology (〈m〉 and Sw ) and two parameters to describe the transport phenomena (h and f).

The energy and momentum equations are rewritten in finite difference form using numerical schemes found in many numerical methods texts and rearranged into a form suitable for solving with an Alternating Direction Implicit (ADI) scheme. Figure 1 shows the numerical simulation flow diagram for flow and heat transfer in porous media using VAT model.

Method of Solution.

For ADI schemes, an implicit equation in each direction is solved at each psuedotime step (or iteration). Hence, for the present two-dimensional case, each time pseudostep is divided into two steps. For example, the solid phase energy equation, the finite difference equation, for the sweep in the vertical direction is Display Formula

and for the horizontal direction Display Formula
Similarly, for the fluid phase, the finite difference equation for the sweep in the vertical direction is Display Formula
and for the horizontal direction Display Formula
In these four equations, the effects of the nonuniform grid, relaxation, morphology, and thermophysical properties of the working solid and liquid materials are grouped together and accounted for in the coefficients A, B, C, D, and F.

Code Validation.

Several different morphologies were tested and the comparisons for both heat transfer and pressure drop were quite good. For optimization, extreme accuracy is not vital because variation in the parameter being optimized can be as much as an order of magnitude. The comparison for heat transfer is shown by Figs.  181920. Note that the different results almost fall on top of one another and that a single correlation could cover the entire range of geometries examined here. This again points to the strength of the VAT or porous media formulation based length scale used in the correlations.

Agreement was also very good for pressure drop although there were differences in the friction factor not absorbed by the use of the VAT based length scale, see Table 1. For a more detailed description of the results of the verification, see Vadnjal [22].

Optimization when there are multiple parameters to deal with can be difficult even with a fast running code. Among the many approaches, two are of interest here. The first is design of experiment (DOE) where a statistically based selection of a number of experiments is established from given ranges of parameters. The second is using genetic algorithms. A future effort will be made to fully incorporate the genetic algorithm method. For this work, however, only DOE will be used. The number of experiments, or in our case the number of computations, depends on the number of parameters that are important.

Once experiments are run, the results are input to FUSION PRO ® DOE software (S-Matrix® Corporation) and the software assigns a weight to each variable and finds a global maximum or minimum for each result. The resulting importance of the various parameters is shown in Fig. 2 where it can be seen that most of the parameters are important making optimization an arduous task. Results are shown in Fig. 2 for staggered pin fins where the approach to an optimum configuration can be seen to be very steep and could easily be missed.

To increase the performance of a plane fin heat sink, the flat surfaces of the fins could be augmented to have a greater heat transfer coefficient and a correspondingly higher friction factor. Of the several possibilities, scales were chosen for their high performance. For scale-roughened walls (see Fig. 2), a correlation from Chang et al.  [26] was used, where x is the distance along the channel and d is the channel diameter. For forward and backward flow, respectively, the correlations are Display Formula

and Display Formula
The Defense Advanced Research Projects Agency (DARPA) Microtechnologies for Air-Cooled Exchangers (MACE) goals were to obtain a thermal resistance of 0.05 C/W, with a pumping power of 33 W, at an effectiveness of 30 and a maximum midpoint temperature of 65 °C. The inlet air temperature is 25 °C and the power applied to the 4 × 4 surface is 1000 W. Plane fin heat sinks were optimized with Re = 70,000 and a fin tip width of 1/3 of the fin base width. The flow was less than 200 cubic feet per minute.

Plane fins, even when tapered, could not meet all the criteria. It was necessary to augment the fin surface with scales. Results from optimization of an augmented plane fin heat sink are shown in Fig. 2. The white area on the figure is where the DARPA MACE criteria are predicted to be met. It was found that the scales had enough increase in heat transfer coefficient without an excessive increase of the friction factor. This result lacked experimental confirmation because the configuration would have been very difficult to achieve.

Note the extent of the white region. There is a great deal of room for variation in the geometric parameters while remaining in the white zone. Further the thermal resistance is very low. It should be noted, however, that these results are overly optimistic because the analysis was done without a base plate. This analysis is presently being redone with the base thickness as another parameter.

Plane fins, even when tapered, could not meet all the criteria. It was necessary to augment the fin surface with scales. Results from optimization of an augmented plane fin heat sink are shown in Fig. 2. The white area in the figure is where the DARPA MACE criteria are predicted to be met. It was found that the scales had enough increase in heat transfer coefficient without an excessive increase of the friction factor. This result lacked experimental confirmation because the configuration would have been very difficult to achieve. Note the extent of the white region. There is a great deal of room for variation in the geometric parameters while remaining in the white zone. Further the thermal resistance is very low. It should be noted, however, that these results are overly optimistic because the analysis was done without a base plate. This analysis is presently being redone with the base thickness as another parameter.

A technique for 3D printing recently developed by Lyons et al.  [27] is being used to construct augmented heat sinks for experimental confirmation of the results shown in Table 2.

The VAT based equations result from averaging over the REV shown in Fig. 2. The REV average velocities are constant values making the equations particularly simple in appearance. It is not necessary to solve the momentum equations once the closure terms have been evaluated. All that must be solved are the energy equations for the two fluids and for the solid material.

The energy equation for fluid 1 is Display Formula

for fluid 2 Display Formula
and for the solid material Display Formula
The simplicity of the equations describing a heat exchanger are actually much simpler than those for a heat sink and do not require the solution of the momentum equation. The reason is that the choice of the REV leads to only needing an average velocity and a properly scaled heat transfer coefficient.

A finned tube heat exchanger was chosen from a Heathcraft Heat Transfer Division brochure that describes their products. Their quotation request form gives the purchaser a great deal of freedom when selecting a design. A flow rate of hot fluid, water, of 1 kg/s with an inlet temperature of 60 °C was chosen. The cold side, the finned side, is air with an inlet temperature of 30 °C. The dimensional choices are wide, and it is not clear how to arrive at a sensible choice when there are ten parameters to work with. The tenth parameter was eliminated by choosing standard tube sizes and their accompanying wall thicknesses. The parameters and their limiting values were taken from a Heatcraft Heat Transfer Division order sheet for a vertical feed fin tube heat exchanger and are given in Table 3.

Each row is a single tube whose serpentine shape becomes the columns of the heat exchanger. FUSION _PRO was used to determine what and how many experimental determinations would be needed for optimization. The number was 110 experiments. These computations were carried out using the simple VAT based computer model. The resulting limiting values of some response variables are shown in Table 4 for aluminum.

It can be seen from Table 4 that you can get whatever you want if you are willing to pay for it. This is further demonstrated in Fig. 2, where it can be seen that the lowest weight results when the pumping power is the highest. Choosing a set of parameters to describe the geometry is a function of what you want. In many cases, power consumption is not an issue and the light weight high power unit would be chosen. The results for copper were very similar.

The lack of a clear optimal design can only be treated with a more meaningful parameter such as cost. This was demonstrated by Peng and Ling [28] using a relatively simple heat exchanger analysis combined with a robust treatment of cost as a measure of an optimal design based on a combination of neural networks augmented with genetic algorithms.

VAT is little more than a judicious application of the theorems by Green and Stokes to carry out the integration needed to average the pointwise conservation equations in a rigorous way. Many everyday engineered devices are hierarchical and heterogeneous and can be effectively treated by application of VAT. It is an approach that can be applied to many different types of transport phenomena, see Travkin and Catton [15]. In this paper, the focus is on conjugate heat transfer using the simple heat sink and finned tube heat exchanger as examples. A printed circuit heat exchanger optimization will follow.

Conservation equations based on rigorous VAT methods were used to determine surface transport processes in support of a search for meaningful heat transfer augmentation. How a two scale heterogeneous heat transfer problem can be solved using exact procedures for closure of additional differential and integral VAT terms was shown. The key to successfully applying VAT is good treatment of the resulting closure. As was seen, this led to some interesting results.

Optimizing a heat sink or a heat exchanger using solutions to the full pointwise energy and momentum equations presents an impossible task given the number of parameters that must be dealt with. Experiments are time consuming and expensive and can only support theory and verify the findings. Without a fast running computer efficient method of calculation, optimization will just not be done.

The importance of the various geometric parameters was a surprise. The importance of the Reynolds number came as no surprise. The importance of all the parameters was unexpected and is a clear demonstration of both the need for the type of modeling developed here and the need to exercise it.

Being able to optimize a heat sink or heat exchanger with complex geometry has typically been avoided because of cost of producing the resulting design. This is no longer the case with 3D printers becoming quite reasonable in cost. Further, 3D printing has been shown by Lyons et al.  [27] to be an effective means to produce any geometry that can be described. We are getting close to being able to develop the ultimate heat sink or heat exchanger.

In a heat transfer problem, the heat transfer coefficient defined locally was shown to be very effective for two reasons. First, it is independent of the type of heating or boundary condition. Second, it was shown that the local heat transfer coefficient is constant throughout the domain and that it is independent of the entrance length. This enables one to broaden the availability of experimental data or facilitate the use of CFD because results from fully developed flow are directly usable.

With closure of the friction factor and the heat transfer coefficient, the problem is closed and the porous media governing equation derived from VAT are

where M̃ stands for averaged momentum equation variables. From the statements above, the macroscale equations are functions only of porous media morphology, represented by porosity and specific surface area, and its closure. Further, it was shown that by proper scaling, closure is a function of the porous media as well, which further generalizes macroscale porous media equations. The simple applications described in this work show that multiple parameters can be dealt with and that heat sinks and heat exchangers have a great deal of room for improvement.

The support of a Department of Energy NERI grant, Award No. DE-FC07-07ID14827, and a DARPA grant in support of the MACE program are gratefully acknowledged. The optimization computations were performed using a commercial version of Fusion Pro® with SC/Tetra® used for needed CFD calculations. The views, opinions, and/or findings contained in this article/presentation are those of the author/presenter and should not be interpreted as representing the official views or policies, either expressed or implied, of the Defense Advanced Research Projects Agency or the Department of Defense.


mean turbulent fluctuation energy (m2 /s2 )


mean drag resistance coefficient in REV


mean form resistance coefficient in REV


mean skin friction coefficient over laminar region in REV


specific heat (J/(kg K))


empirical coefficient


empirical coefficient


tube diameter (m)


internal surface in the REV (m2 )


effectiveness (1/K)


friction factor


averaged over ΔΩf value f


gravitational constant (1/m2 )


height of roughness element (m)


thermal conductivity (W/(m K))


sand grain roughness (cm)


effective thermal conductivity of solid phase (W/(m K))


turbulence mixing length (m)


average porosity


mass flow rate (kg/s)


Nusselt number based on diameter


pitch (nondimensionalized)


pressure (Pa) and pitch in regular porous 2D and 3D medium (m)


Prandtl number


wall heat flux (W/m2 )


heat transfer rate (W)


rib height (nondimensionalized)


Reynolds number based on diameter


specific surface of a porous medium Sw/ΔΩ (1/m)


cross flow projected area of obstacles (m2 )


temperature (K)


velocity in x direction (m/s)


velocity in y direction (m/s)


velocity in z direction (m/s)


averaged heat transfer coefficient over Sw (W/(m2 K))


representative elementary volume (m3 )


pore volume in REV (m3 )


solid phase volume in REV (m3 )


turbulent coefficient exchange ratio


kinematic viscosity (m2 /s)


density (kg/m3 )


fluid phase


component of turbulent vector variable


cross section




solid phase





value in fluid phase averaged over REV


value in solid phase averaged over REV


mean turbulent quantity

turbulent fluctuation quantity

Copyright © 2011 by American Society of Mechanical Engineers
View article in PDF format.



Grahic Jump Location
Figure 12

Development of averaged temperatures at Reh  = 440.8

Grahic Jump Location
Figure 13

Local heat transfer coefficient at Reh  = 440.8

Grahic Jump Location
Figure 15

Comparison of CFD results with experimental data

Grahic Jump Location
Figure 16

Local heat transfer coefficient at different flow rates from Zhou [30]

Grahic Jump Location
Figure 17

VAT solution flow diagram (Hu [29])

Grahic Jump Location
Figure 18

Nusselt number comparison for inline pin fins (Vadnjal [22])

Grahic Jump Location
Figure 19

Nusselt number comparison for staggered pin fins (Vadnjal [22])

Grahic Jump Location
Figure 20

Nusselt number comparison for a bed of spheres (Vadnjal [22])

Grahic Jump Location
Figure 21

Pareto chart for a pin fin heat sink

Grahic Jump Location
Figure 22

Optimum geometric parameters for a staggered pin fin geometry (Hu [29])

Grahic Jump Location
Figure 23

Surface augmentation using scales

Grahic Jump Location
Figure 24

Results for a plane fin with scales

Grahic Jump Location
Figure 25

Response curve for 25.575 mm tubes

Grahic Jump Location
Figure 1

Representative elementary volume (REV)

Grahic Jump Location
Figure 2

REV for a finned tube heat exchanger

Grahic Jump Location
Figure 3

Overlay plot of all the data from [19] using length scales Dc and Dh

Grahic Jump Location
Figure 4

Comparison of heat transfer coefficients for a finned tube heat exchanger using Dc and Dh as the length scales using data from [20] and [21]

Grahic Jump Location
Figure 5

Computational domain morphology

Grahic Jump Location
Figure 6

Comparison of VAT closure term cd by CFD to data

Grahic Jump Location
Figure 7

Contribution of viscous and pressure forces to the drag resistance coefficient

Grahic Jump Location
Figure 8

Comparison of VAT closure Nu obtained by CFD with experiment [31-33]

Grahic Jump Location
Figure 9

Development of averaged temperatures of fluid and solid phase at Reh  = 31.2

Grahic Jump Location
Figure 10

Temperature differences at Reh  = 31.2

Grahic Jump Location
Figure 11

Local heat transfer coefficient at Reh  = 31.2

Grahic Jump Location
Figure 14

Local heat transfer coefficient with different types of heating


Table Grahic Jump Location
Table 3
Variables studied and their ranges
Table Grahic Jump Location
Table 4
Variable responses over range of variables tested using aluminum
Table Grahic Jump Location
Table 1
Closure coefficients for different porous media morphology
Table Grahic Jump Location
Table 2
MACE goals and optimal fin geometry performance


Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In