Flow and Heat Transfer in Graded Porous Media and Its Application in Aeroengine Cooling

Abstract

In recent decades, many studies on flow and heat transfer in porous media have been conducted by researchers to take advantage of its high specific surface area and good heat transfer performance. What is more, the graded porous media have also drawn great attention since the graded arrangement enhances local heat transfer coefficient, so as to improve the overall heat transfer performance and temperature uniformity. The new forms of structure provide a new design of cooling system for some high-power heat sources, such as electronic components, compact heat exchangers, and hot components in aeroengine. In this chapter, the problems of channels filled with vertically graded porous media and axially graded porous media have been introduced, respectively. The flow and heat transfer characteristics of graded porous media have been studied based on porous media theory, considering the parameters of porosity, permeability, pore diameter, etc. The profiles of velocity and temperature vary with different porosity arrangements present. The maximum heat transfer coefficient was obtained for the case of high porosity at the center and low porosity near the wall. Furthermore, the possibility of application in the aeroengine cooling has been discussed.

1. Introduction

Porous media (e.g. metal foam) has drawn great attention due to its potential for application in various scenarios, such as compact heat sinks, regenerative heat exchangers, and combustor-incinerators [1]. Due to its high specific surface area, the metal foam matrix is regarded as a promising material and structure to improve energy absorption and heat transfer in hot component system. In past decades, many relative investigations were conducted theoretically and experimentally and reported for metal foam matrices [2, 3, 4, 5].

Functionally graded porous materials are those porous media with a gradual change of porosity (solidity) and pore size. Thanks to the recent advances in manufacturing technologies, the design and fabrication of such non-uniform structures have become possible. Since the concept of functionally graded materials was proposed in the literature [6] in the 1980s, the studies on making continuous changes in the composition, microstructure, porosity, and other properties have been continuously conducted, since these changes improved its mechanical performance or thermal performance in some reasonable ranges. In recent years, some ceramic-made functionally graded materials have been produced by using a position-dependent chemical composition [7].

The graded porous media may achieve a better heat transfer performance since the graded arrangement enhances local heat transfer coefficient, so as to improve the overall heat transfer performance and temperature uniformity. For example, we may consider some functionally graded metal foam structures that are under equal mean porosity. Firstly, the metal foam matrix may be arranged with higher porosity near the wall subject to constant heat flux and lower porosity in the core regions. This arrangement may create a higher velocity close to the wall and lower velocity in the core, which differs from a regular velocity profile in a channel. On the contrary, lower porosity near the wall may result in thicker thermal boundary layer, which prompts the thermal conductivity near the heated wall but weakens convective heat transfer. Therefore, each case of graded porosity arrangement brings favorable and unfavorable factors on heat transfer performance, but an optimal spatial distribution of the porosity may exist.

On the other hand, axially graded porous materials in a channel with lower porosity (i.e. higher metal volume fraction) in the entrance region and higher porosity in the exit region may further enhance the convective heat transfer coefficient in developing region, however, the convective heat transfer in the downstream maybe weaken because the temperature of the coolant air is increased by mixing. On the contrary, the higher porosity (i.e. lower metal volume fraction) in the entrance region and lower porosity in the exit region may enhance the heat transfer performance of the downstream, since the mixing of the coolant air improves its cooling quality in the downstream. Incidentally, the increase of the heat transfer rate downstream accelerated the heat transfer on the wall and improved the temperature uniformity over the wall.

Wang et al. [8] conducted numerical investigation on the heat transfer of a pipe embedded with functionally graded materials with larger diameter particles near the walls. Zheng et al. [9] have numerically evaluated the heat transfer performance of a pipe filled with a porous foam with radially variable sizes of particles under equal porosity. Bai et al. [10] analytically investigated the heat transfer performance of a channel filled with functionally graded foam matrix in a vertical direction from the central line to the heated walls. They concluded that the graded arrangement of foam matrices may achieve a substantial increase in heat transfer coefficient at the expense of an acceptable consume of pressure loss. Subsequently, Iasiello et al. [11] numerically investigated the effects of variable porosity and cell size on the thermal performance of functionally graded foams, whose porosity and cell size vary according to different power-law functions in the direction of normal to the wall. As a result, a 42% increase in Performance Evaluation Criterion (PEC) is found in foams that account for both variable porosity and cell size. Mauro et al. [12] carried out a comprehensive optimization framework to design a graded foam structure by mono-and multi-objective optimization method. They proposed optimum design parameters that achieve the best thermal performance as high as 1.5 performance evaluation criterion (PEC). Recently, Bai et al. [13] considered an occasion of heat transfer in a channel filled with axially graded porous materials to investigate the temperature distribution on the heated walls. They found that the wall temperature variation is quite sensitive to the axial gradient of the local porosity. Furthermore, the functionally graded foams with variable PPI enhance pool boiling heat transfer and that with linearly-variable porosity enhance heat transfer in phase change materials were presented by Xu and Zhao [14] and Yang et al. [15], respectively. In short, the enhancement of heat transfer by exploiting graded porous materials has been comprehensively substantiated.

With the increase of the demand on the large thrust of the engine, the temperature value of gas in front of turbine increases. For the military aircraft, the afterburner is adopted to increase the transient trust of the aeroengine. However, the gas temperature at the afterburning state far exceeds the limitation of the material. In this case, the graded porous structure may be useful to adjust the coolant distribution and to improve the temperature uniformity.

In this chapter, the flow and heat transfer in channels filled with graded porous materials (radially graded and axially graded) have been investigated, and the flow structures and heat transfer characteristics have been revealed for both cases with different parameters. The possibility of the application of graded porous structure in afterburner heat shield has been discussed.

2. Flow and heat transfer in a channel filled with vertically graded porous media

The analytical investigation of heat transfer characteristics of a channel filled with radially graded porous media has been conducted in this section. A series of graded metal foam matrices of the same solidity were arranged with their porosity either increasing or decreasing toward the heated walls following a parabolic function. As a result, the maximum heat transfer coefficient of the channel filled with a functionally graded metal foam matrix was achieved for the case that low porosity near the heated wall and high porosity in the central line, which is 20–50% higher than that of the uniform distribution of metal foam matrix. Therefore, the effectiveness of functionally graded metal foam matrices to achieve higher heat transfer coefficient with an acceptable pressure loss was substantiated.

2.1 Physical model and analysis

As illustrated in Figure 1, there are two cases of channels filled with a functionally graded metal matrix, one case has high porosity near the wall and low porosity in the core region of channel, the other case has reversed arrangement. Both upper and lower walls were subject to a constant heat flux $q_w$. We mainly focus on the region of hydro-dynamically and thermally fully developed flow in these cases, because the developing region, which may take several hydraulic diameters, is short enough due to the addition of the metal foam matrix.

Figure 1.

Physical models of channels with graded porous media. (a) The case with high porosity near the wall and low porosity in the core and (b) the case with low porosity near the wall and high porosity in the core.

The volume-averaged version of the momentum equation was derived by exploiting the Brinkman-extended-Darcy momentum equation Eq. (1), likewise, the volume-averaged energy equations for the fluid and solid phases were obtained as Eqs. (2) and (3).

$$-\frac{\mathit{dp}}{\mathit{dx}}+\frac{d}{\mathit{dy}}\left(\frac{\mu }{\epsilon }\frac{\mathit{du}}{\mathit{dy}}\right)-\frac{\mu }{K}u=0$$E1

$${\mathrm{\rho }}_f{\mathrm{c}}_{\mathit{pf}}u\frac{\partial T_f}{\partial x}=\frac{\partial }{\partial y}\left({\epsilon }^{\ast }{k}_f+{\mathit{\epsilon k}}_{\mathit{dis}}\right)\frac{\partial T_f}{\partial y}-h_v\left(T_f-T_s\right)$$E2

$$\frac{\partial }{\partial y}\left(1-{\epsilon }^{\ast }\right)k_s\frac{\partial T_s}{\partial y}+h_v\left(T_f-T_s\right)=0$$E3

where $u$ is the Dacian bulk velocity, which is related to the local porosity and intrinsic velocity. $\mathrm{\epsilon }\left(y\right)$ and $K\left(\epsilon \right)$ are porosity and permeability, respectively, which vary spatially across the channel in vertical direction(y), whose origin is located on the lower heated wall. $T_f$ and $T_s$ are the intrinsic volume average fluid and solid phase temperatures, respectively. $h_v$ is the volumetric interstitial heat transfer coefficient. In general, the thermal dispersion term is often expressed by the gradient diffusion hypothesis with the thermal dispersion conductivity $k_{\mathit{dis}}$ [16]. The concept of the effective porosity ${\mathrm{\epsilon }}^{\ast }$ was introduced to account for the tortuosity of the structure, which is theoretically given as

$${\mathrm{\epsilon }}^{\ast }\equiv \frac{k_s-k_e}{k_s-k_f}\cong \frac{2+\epsilon }{3}$$E4

where$k_e$ is the effective stagnant thermal conductivity, which is given as

$$k_e\left(\epsilon \right)={\epsilon }^{\ast }{k}_f+\left(1-{\epsilon }^{\ast }\right)k_s=\frac{2+\epsilon }{3}{k}_f+\frac{1-\epsilon }{3}{k}_s$$E5

The two energy equations, namely Eqs. (3) and (4), may be added to yield

$${\mathrm{\rho }}_f{\mathrm{c}}_{\mathit{pf}}u\frac{\partial T_f}{\partial x}=\frac{\partial }{\partial y}\left(\left({\epsilon }^{\ast }{k}_f+{\mathit{\epsilon k}}_{\mathit{dis}}\right)\frac{\partial T_f}{\partial y}+\left(1-{\epsilon }^{\ast }\right)k_s\frac{\partial T_s}{\partial y}\right)$$E6

Then, Eq. (6) may be reduced under the local thermal equilibrium assumption and neglecting the comparatively small term of the thermal dispersion conductivity$k_{\mathit{dis}}$ for the case of aluminum-air combination [17] to

$${\mathrm{\rho }}_f{\mathrm{c}}_{\mathit{pf}}u\frac{\partial T}{\partial x}=\frac{\partial }{\partial y}\left(k_e\left(\epsilon \right)\frac{\partial T}{\partial y}\right)$$E7

Eq. (7) should be solved with the momentum Eq. (1) under the following boundary conditions:

For $y=0$:

$$u=0\text{and}{q}_w=-k_e\left({\epsilon }_w\right)\frac{\partial T}{\partial y}$$E8

For $y=H$:

$$\frac{\mathrm{\partial u}}{\partial y}=0\text{and}\frac{\partial T}{\partial y}=0$$E9

The porosity of the functionally graded metal foam matrix in a channel varies following a parabolic function across the channel as follows:

$$\epsilon \left(y\right)={\epsilon }_c+\left({\epsilon }_w-{\epsilon }_c\right){\left(1-\frac{y}{H}\right)}^2$$E10

where the subscripts of w and c indicate the wall and the center of the channel, respectively.

As for the permeability, the correlation provided by Calmidi [18] for metal foam was adopted:

$$K\left(\epsilon \right)=0.000365{\left(1-\epsilon \right)}^{-0.224}{\left(\frac{1.18}{1-e^{-\left(1-\epsilon \right)/0.04}}\sqrt{\frac{1-\epsilon }{3\pi }}\right)}^{-1.11}{H}^2$$E11

where $K\left(\epsilon \right)/H^2$ is function of the porosity, which typically ranges from 3 × 10^–3 to 8 × 10^–3.

Subsequently, the non-dimensional momentum equation and energy equation were solved using a standard Runge-Kutta-Gill scheme, with the following dimensionless boundary conditions:

$${\left.u^{\ast }\right|}_{\eta =0}=0$$E12

and

$${\left.\frac{{\mathit{du}}^{\ast }}{\mathit{d\eta }}\right|}_{\eta =1}=0$$E13

Once the velocity profile was determined and fed into the energy Eq. (6), the following expression can be yielded by integrating from y = 0 to H

$$H{\rho }_f{c_p}_f{u}_B\frac{d{T}_B}{\mathit{dx}}=q_w$$E14

where

$$T_B=\frac{{\rho }_f{c}_{\mathit{pf}}{\int }_0^H\mathit{\text{uTdy}}}{{\rho }_f{c}_{\mathit{pf}}{u}_B}$$E15

is the bulk mean temperature. Since $d{T}_B/\mathit{dx}=\partial T/\partial x$ for the thermally fully developed flow, the energy Eq. (6) simplifies to

$$\frac{u^{\ast }}{u{\ast }_B}=\frac{d}{\mathit{d\eta }}\left(\frac{k_e\left(\epsilon \right)}{k_e\left({\epsilon }_B\right)}\frac{{\mathit{dT}}^{\ast }}{\mathit{d\eta }}\right)$$E16

where

$$u{\ast }_B={\int }_0^1{u}^{\ast }d\eta$$E17

and

$$T^{\ast }=\frac{k_e\left({\epsilon }_B\right)\left(T-T_w\right)}{q_{w}H}$$E18

are the dimensionless bulk mean velocity and local dimensionless temperature, respectively. Integrate the Eq. (16) from $\eta =0$ can attain

$$\frac{{\mathit{dT}}^{\ast }}{\mathit{d\eta }}=\frac{{\int }_0^{\eta }\frac{u^{\ast }}{u{\ast }_B}\mathit{d\eta }-1}{\frac{k_e\left(\epsilon \right)}{k_e\left({\epsilon }_B\right)}}$$E19

The dimensionless temperature can be found by further integrating the foregoing Eq. (19) by using the Runge-Kutta-Gill scheme. The following boundary conditions on the temperature and its gradient on the wall are implemented:

$${\left.T^{\ast }\right|}_{\eta =0}=0\text{and}{\left.\frac{{\mathit{dT}}^{\ast }}{d\eta }\right|}_{\eta =0}=-\frac{k_e\left({\epsilon }_B\right)}{k_e\left({\epsilon }_w\right)}$$E20

Thus, the Nusselt number of our primary interest can be obtained from

$${\mathit{Nu}}_H\equiv \frac{q_{w}H}{k_e\left({\epsilon }_B\right)\left(T_w-T_B\right)}=-\frac{1}{{\int }_0^1\frac{u^{\ast }}{u{\ast }_B}{T}^{\ast }\mathit{d\eta }}=\frac{1}{{\int }_0^1\frac{\left({\int }_0^{\eta }\frac{u^{\ast }}{u{\ast }_B}\mathit{d\eta }\right)\left({\int }_0^{\eta }\frac{u^{\ast }}{u{\ast }_B}\mathit{d\eta }-1\right)}{\frac{k_e\left(\epsilon \right)}{k_e\left({\epsilon }_B\right)}}\mathit{d\eta }-{\left.T^{\ast }\right|}_{\eta =1}}$$E21

The definition of ${\mathit{Nu}}_H$ contains $k_e\left({\epsilon }_B\right)$ accounted for the effect of increase in the effective thermal conductivity on the Nusselt number. It is meaningful to compare ${\mathit{Nu}}_H$ for all cases of different weights, even for the case of the channel without a metal foam, based on the definition given by Eq. (21).

2.2 Velocity and temperature profiles

The dimensionless velocity profiles obtained for the metal foams of the same weight of ${\epsilon }_B$ = 0.933 were presented in Figure 2. The wall porosity ${\epsilon }_w$ varies from 0.840–0.980, and the case of ${\epsilon }_w$ = 0.933 corresponds to the case of uniform porosity metal foam. As we can see, the increase in the ${\epsilon }_w$ creates a higher velocity value close to the wall and lower velocity with a concave profile in the core region.

Figure 2.

Dimensionless velocity profiles.

To the contrary, the decrease of the ${\epsilon }_w$ results in a higher velocity value in the core region and lower velocity value near the wall comparing with the velocity distribution for the case of the uniform porosity. Naturally, more metal near the wall has favorable influence on conductive wall heat transfer, but more metal may block the near-wall flow, decreasing the flow velocity. Therefore, these favorable and unfavorable factors form a trade-off for the heat transfer enhancement.

The corresponding dimensionless temperature profile for the case of $k_s/k_f$ = 10⁴ was shown in Figure 3, which assumes a combination of air and aluminum.

Figure 3.

Dimensionless temperature profiles.

As can be seen, the dimensionless temperature profile of the uniform porosity case of ${\epsilon }_w$ = ${\epsilon }_w$ = 0.933 has been presented and indicated. As increasing the value of ${\epsilon }_w$, the dimensionless temperature gradient at the wall become steeper, which can also be expected from Eq. (20). While decreasing the ${\epsilon }_w$, the dimensionless temperature profile near the wall changes its curvature to negative one, while the protrusion of the dimensionless temperature in the core becomes more extensive.

2.3 Heat transfer characteristics

Figure 4 presents the curves of ${\mathit{Nu}}_H$ against the porosity distribution for the given bulk mean porosity ${\epsilon }_B$ = 0.900, 0.933, and 0.950. The level of ${\mathit{Nu}}_H$ for the channel without a metal foam, namely, ${\mathit{Nu}}_H$ = 2.06, was also indicated. The figure indicates that the increase in the effective thermal conductivity enables the convective heat transfer coefficient in a channel filled with a metal foam to increase about 20–50% more than 2.06$k_e\left({\epsilon }_B\right)/H$. However, for the cases of positive ${\epsilon }_w$-${\epsilon }_B$, the values of ${\mathit{Nu}}_H$ go below the value 2.06.

Figure 4.

Nusselt number ${\mathit{Nu}}_H$ against $\left({\epsilon }_{\mathrm{w}}-{\epsilon }_B\right).$

As a result, the maximum values of ${\mathit{Nu}}_H$ are found in the negative range of ${\epsilon }_w–{\epsilon }_B$, where the porosity of the functionally graded metal foam decreases toward the wall (i.e. more metal near the wall).

In comparison with the case of the channel without a metal foam, the addition of the metal foam achieves high heat transfer coefficient but at expense of a larger pressure drop, or pressure gradient $\left(-\mathit{dp}/\mathit{dx}\right)$, which can be related to its pumping power per unit axial length $H{u}_B\left(-\mathit{dp}/\mathit{dx}\right)$ as

$$-\frac{\mathit{dp}}{\mathit{dx}}=\sqrt{\frac{1}{u{\ast }_B}\frac{\mu }{H^3}\left({\mathit{Hu}}_B\left(-\frac{\mathit{dp}}{\mathit{dx}}\right)\right)}$$E22

Thus, the pressure drop ratio in the channel with a metal foam to that of an empty channel (where ${u^{\ast }}_f=1/3$) under equal pumping power is given by

$$\mathrm{\Delta }p/\mathrm{\Delta }{p}_f=1/\sqrt{3u{\ast }_B}$$E23

As shown in Figure 5, the pressure drop ratios for the corresponding three cases are presented. In comparison to the case of the channel without a metal foam, the addition of the metal foam matrix results in about 7–8 times increase in pressure drop under equal pumping power. Moreover, the bulk mean porosity of ${\epsilon }_B$ brings more sensitive effects than $\left({\epsilon }_{\mathrm{w}}-{\epsilon }_B\right)$.

Figure 5.

Pressure drop ratio $∆p/∆p_f$ against $\left({\epsilon }_{\mathrm{w}}-{\epsilon }_B\right).$

2.4 Conclusions

In this section, the heat transfer characteristics of forced convection in channels filled with functionally graded metal foam in vertical direction normal to the heated walls have been analytically investigated. The fully developed set of Brinkman-extended Darcy momentum equation and the energy equation associated with local thermal equilibrium assumption were taken into account. The porosity of parabolic distribution was assumed across the channel height. The integrated energy equation associated with momentum equation was solved by using the Rung-Kutta-Gill integration scheme. According to the analytical results, for the metal foams of the same weight, the maximum heat transfer coefficient can be obtained for the case of the high porosity in the core and low porosity near the wall, at the expense of an acceptable increase in pressure drop.

3. Flow and heat transfer in a channel filled with axiallygraded porous media

In this section, the problem of flow and heat transfer in a channel filled with axially graded porous materials based on local thermal non-equilibrium has been investigated to find a way to control its wall temperature distribution [13]. A series of axially graded foams under an equal solid fraction were considered. The porosity decreasing along the downstream of the channel following a power function of axial coordinate with an exponent of n was assumed. The Brinkman-Forchheimer extended Darcy model and the local thermal non-equilibrium assumption were exploited to obtain the velocity profiles and temperature distributions on the heated walls, respectively.

3.1 Physical model and analysis

A channel of height 2H filled with a fluid-saturated axially graded porous media subject to constant heat flux is presented in Figure 6.

Figure 6.

Channel model filled with an axially graded porous medium.

The indexes of ${\delta }_f$ and ${\delta }_s$ represent thermal boundary layer thicknesses for fluid phase and solid phase, respectively. Both the fluid phase thermal boundary layer and solid phase thermal boundary layer grow from the entrance to the downstream. However, the thermal boundary layer of solid grows thicker than that of the fluid due to the higher thermal conductivity of the solid than that of the fluid. In the present study, the thermally developing region of the entrance is focused so that ${\delta }_f$ and ${\delta }_s$ are less than half the channel height, H.

In the present case, the Brinkman-Forchheimer-extended Darcy model is still valid since the velocity change in such an axially variable porosity-arranged porous media is so moderate that the convective inertia terms are almost negligible. Hence, the velocity field can be determined from the momentum and continuity equations as follows:

$$-\frac{\mathit{dp}}{\mathit{dx}}+\frac{{\mu }_f}{\epsilon \left(x\right)}\frac{{\partial }^{2}u}{\partial y^2}-\frac{{\mu }_f}{K\left(x\right)}u-{\rho }_{f}b\left(x\right)u^2=0$$E24

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$E25

where $u$ and $v$ are the local volume average velocity components. $p$ is the intrinsic pressure, ${\mu }_f$ is the fluid viscosity, $K\left(x\right)$ is the permeability, and $b$ is the inertial coefficient. Moreover, x is the axial coordinate while y is the coordinate measured vertically from the lower wall surface of the channel.

Following Amiri and Vafai [19], the energy equations for the fluid phase and the solid phase can be given as

$${\mathrm{\rho }}_f{\mathrm{c}}_{p_f}\frac{\partial {\mathit{uT}}_f}{\partial x}+{\mathrm{\rho }}_f{\mathrm{c}}_{p_f}\frac{\partial {\mathit{vT}}_f}{\partial y}=\frac{\partial }{\partial y}{k}_{\mathit{fe}}\left(x\right)\frac{\partial T_f}{\partial y}-h_v\left(x\right)\left(T_f-T_s\right)$$E26

$$\frac{\partial }{\partial y}{k}_{\mathit{se}}\left(x\right)\frac{\partial T_s}{\partial y}+h_v\left(x\right)\left(T_f-T_s\right)=0$$E27

where $k_{\mathit{fe}}$ and $k_{\mathit{se}}$ are the effective thermal conductivities of the fluid and the solid, respectively. $h_v$ denotes the interstitial volumetric heat transfer coefficient.

Considering the symmetry of the geometry, the lower half of the channel is analyzed with corresponding boundary conditions given as follows:

$$x=0:T_f=T_0=\text{const}.$$E28

$$y=0:u=v=0,T_f=T_s=T_w\left(x\right)\text{and}{q}_w=-k_{\mathit{fe}}{\left.\frac{\partial T_f}{\partial y}\right|}_{y=0}-k_{\mathit{se}}{\left.\frac{\partial T_s}{\partial y}\right|}_{y=0}=\text{const}.$$E29

$$y=H:v=0,{\left.\frac{\partial T_f}{\partial y}\right|}_{y=H}={\left.\frac{\partial T_s}{\partial y}\right|}_{y=H}=0$$E30

where $q_w$is the wall heat flux. $T_w\left(x\right)$ is the wall temperature.

According to Yi et al. [20], the Eq. (24) can be normalized by using the constant bulk velocity $u_B$ with axially variable porosity and permeability as

$$\frac{{\partial }^2}{\partial {\eta }^2}\left(\frac{u}{u_B}\right)-\left(\frac{\epsilon \left(x\right)}{\mathit{Da}\left(x\right)}\right)\left(\frac{u}{u_B}\right)-{\mathit{Re}}_b\left(x\right){\left(\frac{u}{u_B}\right)}^2=-\left(-\frac{\epsilon H^2}{{\mu }_f{u}_B}\frac{\mathit{dp}}{\mathit{dx}}\right)$$E31

where

$$\eta =\frac{y}{H}$$E32

$$\mathit{Da}\left(x\right)=\frac{K\left(x\right)}{H^2}\left(\text{Darcy number}\right)$$E33

$${\mathit{Re}}_b\left(x\right)=\frac{\epsilon \left(x\right){\rho }_{f}b\left(x\right)H^2{u}_B}{{\mu }_f}$$E34

where ${\mathit{Re}}_b\left(x\right)$ is the Forchheimer Reynolds number based on the scale of length $b{H}^2.$ The following velocity profile $u\left(x,\eta \right)$ may be assumed with an unknown parameter$\zeta \left(x\right)$:

$$\frac{u}{u_B}=\frac{\zeta +1}{\zeta }\left(1-{\left(1-\eta \right)}^{\zeta }\right)$$E35

which automatically satisfies the boundary conditions (29) and ${\int }_0^1\left(\frac{u}{u_B}\right)\mathit{d\eta }=1$. It is obvious that the case of $\zeta =2$ in Eq. (35) denotes the fully developed profile in an empty channel flow, while for the case of $\zeta \to \infty$, the velocity profile becomes the Darcy(plug) flow.

Evaluating Eq. (31) on the wall (at $\eta =0$) using the foregoing velocity profile (35), one obtains

$$-\frac{\epsilon H^2}{{\mu }_f{u}_B}\frac{\mathit{dp}}{\mathit{dx}}={\zeta }^2-1$$E36

Thus, Eq. (31) may be integrated across the channel from $\eta =0$ to 1 as

$${\left.-\frac{\partial }{\partial \eta }\left(\frac{u}{u_B}\right)\right|}_{\eta =0}-\left(\frac{\epsilon }{\mathit{Da}}\right){\int }_0^1\left(\frac{u}{u_B}\right)\mathit{d\eta }-{\mathit{Re}}_b{\int }_0^1{\left(\frac{u}{u_B}\right)}^2\mathit{d\eta }=-\left({\zeta }^2-1\right)$$E37

Based on the Eq. (37) associated with the Eq. (31), a cubic equation for the unknown parameter $\zeta$ can be yielded as follows:

$$2{\zeta }^3-{\zeta }^2-\left(2\frac{\epsilon }{\mathit{Da}}+2{\mathit{Re}}_b+5\right)\zeta -\left(\frac{\epsilon }{\mathit{Da}}+2{\mathit{Re}}_b+2\right)=0$$E38

The local value of the $\zeta \left(\epsilon /\mathit{Da},{\mathit{Re}}_b\right)$ in this cubic equation can be easily evaluated by substituting the local values of $\epsilon \left(x\right)/\mathit{Da}\left(x\right)$ and${\mathit{Re}}_b\left(x\right)$. As for the root of the foregoing cubic equation, a highly accurate expression in a closed form based on the two asymptotes was obtained by Yi et al. [20]:

$$\zeta \left(\left(\frac{\epsilon }{\mathit{Da}}\right),{\mathit{Re}}_b\right)=\frac{{\mathit{Re}}_b\frac{3+\sqrt{25+16\left({\mathit{Re}}_b+\frac{\epsilon }{\mathit{Da}}\right)}}{4}+\left(\frac{\epsilon }{\mathit{Da}}\right)\frac{1+\sqrt{9+4\left(\frac{\epsilon }{\mathit{Da}}+{\mathit{Re}}_b\right)}}{2}}{\left(\frac{\epsilon }{\mathit{Da}}\right)+{\mathit{Re}}_b}$$E39

Eq. (39) reduces to the two distinct exact asymptotes as:

$$\zeta \left(\frac{\epsilon }{\mathit{Da}},0\right)=\frac{1+\sqrt{9+4\frac{\epsilon }{\mathit{Da}}}}{2}$$E40

and

$$\zeta \left(0,{\mathit{Re}}_b\right)=\frac{3+\sqrt{25+16{\mathit{Re}}_b}}{4}$$E41

Naturally, the velocity profile function with Eqs. (40) and (41) gives $\frac{u}{u_B}=\frac{3}{2}\left(1-{\left(1-\eta \right)}^2\right)$ for $\frac{\mathit{Da}}{\epsilon }\to \infty$ (or${\mathit{Re}}_b\to 0$) (i.e. empty channel) and $\frac{u}{u_B}=1$for $\frac{\mathit{Da}}{\epsilon }\to 0$ (or${\mathit{Re}}_b\to \infty$) (i.e. Darcy flow).

Once $\zeta \left(x\right)$ is determined, the local friction coefficient may be estimated from

$$\frac{H}{{\rho }_f{u_B}^2}\left(-\frac{\mathit{dp}}{\mathit{dx}}\right)=\frac{{\zeta }^2-1}{\epsilon {\mathit{Re}}_H}$$E42

where

$${\mathit{Re}}_H=\frac{{\rho }_f{u}_{B}H}{{\mu }_f}$$E43

is the Reynolds number based on $u_B$ and $H$.

The vertical velocity component $v\left(x,\eta \right)$ may be determined by substituting Eq. (31) into the continuity Eq. (25) as

$$\begin{array}{rcl}\frac{v}{u_B}& =& -H\frac{\partial }{\partial x}{\int }_0^{\eta }\frac{u}{u_B}\mathit{d\eta }=-H\frac{\partial }{\partial x}\frac{\zeta +1}{\zeta }{\int }_0^{\eta }\left(1-{\left(1-\eta \right)}^{\zeta }\right)\mathit{d\eta }\\ & =& -\frac{H}{{\zeta }^2}\frac{\mathit{d\zeta }}{\mathit{dx}}\left(1-\eta \right)\left(1-{\left(1-\eta \right)}^{\zeta }\left(1-\zeta \mathit{ln}\left(1-\eta \right)\right)\right)\end{array}$$E44

which satisfies the boundary conditions (29). It is interesting to note that the velocity component $v$ in the lower half of the channel is negative for the channel of $\frac{\mathit{d\zeta }}{\mathit{dx}}$ > 0, which is met when either $\frac{\epsilon }{\mathit{Da}}=\frac{\epsilon \left(x\right)H^2}{K\left(x\right)}$ or ${\mathit{Re}}_b$=$\frac{\epsilon \left(x\right){\rho }_{f}b\left(x\right)H^2{u}_B}{{\mu }_f}$ increases downstream, while the vertical velocity component $v$ vanishes in the case of channels filled with uniform porous materials ($\mathrm{i}.\mathrm{e}.\frac{\mathit{d\zeta }}{\mathit{dx}}=0$).

For the thermally fully developed flows in a channel filled with a uniform porous medium, Nakayama et al. [4] obtained a temperature field under the L.T.E. assumption as

$$\begin{array}{rcl}\frac{T_f-T_w}{T_{\mathit{fB}}-T_w}& =& \frac{T_s-T_w}{T_{\mathit{fB}}-T_w}\\ & =& \frac{3\left(\zeta +3\right)\left(2\zeta +3\right)}{2\zeta \left(\zeta +2\right)\left(2{\zeta }^2+13\zeta +17\right)}(\left(\zeta +1\right)\left(\zeta +2\right)\left(1-{\left(1-\eta \right)}^2\right)\\ & & -2\left(1-{\left(1-\eta \right)}^{\zeta +2}\right))\end{array}$$E45

For the case of the channel filled with axially graded porous materials, analytical expressions for its wall temperature under the L.T.N.E. assumption were obtained by exploiting the above temperature profile (Eq. (45)) and velocity profile (Eq. (35)), which prevail across both fluid and solid thermal boundary layers (${\delta }_f<{\delta }_s<H$):

$$\frac{T_f-T_0}{T_w-T_0}=\left\{\begin{array}{c}1-\frac{\left(\zeta +1\right)\left(\zeta +2\right)\left(1-{\left(1-\frac{\eta }{{{\delta }_f}^{\ast }}\right)}^2\right)-2\left(1-{\left(1-\frac{\eta }{{{\delta }_f}^{\ast }}\right)}^{\zeta +2}\right)}{\zeta \left(\zeta +3\right)}:0\le \eta \le {{\delta }_f}^{\ast }\\ 0:{{\delta }_f}^{\ast }\le \eta \le 1\end{array}\right.$$E46

and

$$\frac{T_s-T_0}{T_w-T_0}=\left\{\begin{array}{c}1-\frac{\left(\zeta +1\right)\left(\zeta +2\right)\left(1-{\left(1-\frac{\eta }{{{\delta }_s}^{\ast }}\right)}^2\right)-2\left(1-{\left(1-\frac{\eta }{{{\delta }_s}^{\ast }}\right)}^{\zeta +2}\right)}{\zeta \left(\zeta +3\right)}:0\le \eta \le {{\delta }_s}^{\ast }\\ 0:{{\delta }_s}^{\ast }\le \eta \le 1\end{array}\right.$$E47

where

$${{\delta }_f}^{\ast }={\delta }_f/H$$E48

and

$${{\delta }_s}^{\ast }={\delta }_s/H$$E49

These profiles satisfied the boundary conditions (29) and (30) and naturally recovered the fully-developed profile given by Eq. (45) as both ${{\delta }_f}^{\ast }$ and ${{\delta }_s}^{\ast }$ grow downstream to reach unity. The fluid bulk mean temperature $T_{\mathit{fB}}$ in the entrance region may be estimated using the profiles (35) and (46) as

$$\frac{T_{\mathit{fB}}-T_0}{T_w-T_0}\equiv {\int }_0^1\frac{u}{u_B}\left(\frac{T_f-T_0}{T_w-T_0}\right)d{y}^{\ast }$$

$$={\int }_0^{{{\delta }_f}^{\ast }}\frac{\zeta +1}{\zeta }\left(1-{\left(1-y^{\ast }\right)}^{\zeta }\right)\left(1-\frac{\left(\zeta +1\right)\left(\zeta +2\right)\left(1-{\left(1-\frac{y^{\ast }}{{{\delta }_f}^{\ast }}\right)}^2\right)-2\left(1-{\left(1-\frac{y^{\ast }}{{{\delta }_f}^{\ast }}\right)}^{\zeta +2}\right)}{\zeta \left(\zeta +3\right)}\right)d{y}^{\ast }$$

$$\cong {\int }_0^{\mathit{min}\left(\frac{1}{\zeta },{{\delta }_f}^{\ast }\right)}\frac{\zeta +1}{\zeta }\left(\zeta y^{\ast }\right)\left(1-\frac{\left(\zeta +1\right)\left(\zeta +2\right)\left(1-{\left(1-\frac{y^{\ast }}{{{\delta }_f}^{\ast }}\right)}^2\right)-2\left(1-{\left(1-\frac{y^{\ast }}{{{\delta }_f}^{\ast }}\right)}^{\zeta +2}\right)}{\zeta \left(\zeta +3\right)}\right)d{y}^{\ast }$$

$$+{\int }_{\mathit{min}\left(\frac{1}{\zeta },{{\delta }_f}^{\ast }\right)}^{{{\delta }_f}^{\ast }}\frac{\zeta +1}{\zeta }\left(1-\frac{\left(\zeta +1\right)\left(\zeta +2\right)\left(1-{\left(1-\frac{y^{\ast }}{{{\delta }_f}^{\ast }}\right)}^2\right)-2\left(1-{\left(1-\frac{y^{\ast }}{{{\delta }_f}^{\ast }}\right)}^{\zeta +2}\right)}{\zeta \left(\zeta +3\right)}\right)d{y}^{\ast }$$

$$\begin{array}{l}=\frac{\left(\zeta +1\right)\left(\zeta +5\right)\left({\zeta }^2+5\zeta +10\right)}{12{\left(\zeta +3\right)}^2\left(\zeta +4\right)}{{{\delta }_f}^{\ast }}^2\\ -\frac{{\left(\zeta +1\right)}^2\left(\zeta +2\right)}{12\mathrm{\zeta }\left(\zeta +3\right)}{\left(1-\frac{\mathit{min}\left(\frac{1}{\zeta },{{\delta }_f}^{\ast }\right)}{{{\delta }_f}^{\ast }}\right)}^4\left(1-\frac{24{\left(1-\frac{\mathit{min}\left(\frac{1}{\zeta },{{\delta }_f}^{\ast }\right)}{{{\delta }_f}^{\ast }}\right)}^{\zeta }}{\left(\zeta +1\right)\left(\zeta +2\right)\left(\zeta +3\right)\left(\zeta +4\right)}\right){{{\delta }_f}^{\ast }}^2\end{array}$$E50

Since the fluid phase thermal boundary layer is thinner in this entrance developing region, the velocity field across the thermal boundary layer may be described by a linearized velocity profile.

The interstitial heat transfer term $h_v\left(x\right)\left(T_f-T_s\right)$ can be canceled by adding the energy Eqs. (26) and (27). Moreover, integrating the resulting equation across the channel and over the axial distance from 0 to x with the boundary conditions (28) and (29), the energy balance relationship can be obtained as follows:

$${\mathrm{\rho }}_f{\mathrm{c}}_{p_f}{u}_{B}H\left(T_{\mathit{fB}}-T_0\right)=q_{w}x$$E51

where

$$q_w=-k_{\mathit{fe}}{\left.\frac{\partial T_f}{\partial y}\right|}_{y=0}-k_{\mathit{se}}{\left.\frac{\partial T_s}{\partial y}\right|}_{y=0}=\frac{\left(T_w-T_0\right)}{H}\frac{2\left(\zeta +2\right)}{\left(\zeta +3\right)}\left(\frac{k_{\mathit{fe}}}{{{\delta }_f}^{\ast }}+\frac{k_{\mathit{se}}}{{{\delta }_s}^{\ast }}\right)$$E52

Furthermore, integrating the solid phase heat conduction Eq. (4) from y = 0 to H with the boundary condition (30) may yield

$$-k_{\mathit{se}}\left(x\right){\frac{\partial T_s}{\partial y}⌋}_{y=0}+h_v\left(x\right){\int }_0^H\left(T_f-T_s\right)\mathit{dy}=0$$E53

The temperature profile functions (46) and (47) are fed into the foregoing equation to find

$${{\delta }_s}^{\ast }\left({{\delta }_{\mathrm{s}}}^{\ast }-{{\delta }_f}^{\ast }\right)=\frac{6}{\mathit{Bi}\left(x\right)}\frac{\left(\zeta +2\right)\left(\zeta +3\right)}{\left({\zeta }^2+6\zeta +11\right)}$$E54

which gives the solid phase thermal boundary layer thickness ${{\delta }_s}^{\ast }$:

$$\frac{{{\delta }_s}^{\ast }}{{{\delta }_f}^{\ast }}=\frac{1}{2}\left(1+\sqrt{1+\frac{24}{\mathit{Bi}\left(x\right)}\frac{\left(\zeta +2\right)\left(\zeta +3\right)}{\left({\zeta }^2+6\zeta +11\right){{{\delta }_f}^{\ast }}^2}}\right)$$E55

where

$$\mathit{Bi}\left(x\right)=\frac{h_v\left(x\right)H^2}{k_{\mathit{se}}\left(x\right)}$$E56

is the local Biot number accounts for the interstitial convective conductance and the solid phase conductance.

Eqs. (50)-(52), and (55) are combined to eliminate$\left(T_{\mathit{fB}}-T_0\right)$, $q_w$ and ${{\delta }_s}^{\ast }$ resulting in the following relationship between ${{\delta }_f}^{\ast }$ and $\xi$:

$$\xi =\frac{\frac{\left(\zeta +1\right)\left(\zeta +5\right)\left({\zeta }^2+5\zeta +10\right)}{12{\left(\zeta +3\right)}^2\left(\zeta +4\right)}{{{\delta }_f}^{\ast }}^2-\frac{{\left(\zeta +1\right)}^2\left(\zeta +2\right)}{12\mathrm{\zeta }\left(\zeta +3\right)}{\left(1-\frac{\mathit{min}\left(\frac{1}{\zeta },{{\delta }_f}^{\ast }\right)}{{{\delta }_f}^{\ast }}\right)}^4\left(1-\frac{24{\left(1-\frac{\mathit{min}\left(\frac{1}{\zeta },{{\delta }_f}^{\ast }\right)}{{{\delta }_f}^{\ast }}\right)}^{\zeta }}{\left(\zeta +1\right)\left(\zeta +2\right)\left(\zeta +3\right)\left(\zeta +4\right)}\right){{{\delta }_f}^{\ast }}^2}{\frac{2\left(\zeta +2\right)}{\left(\zeta +3\right){{\delta }_f}^{\ast }}\left(\frac{k_{\mathit{fe}}\left(x\right)}{k_{\mathit{ref}}}+\frac{k_{\mathit{se}}\left(x\right)}{k_{\mathit{ref}}}\frac{2}{1+\sqrt{1+\frac{24}{\mathit{Bi}}\frac{\left(\zeta +2\right)\left(\zeta +3\right)}{\left({\zeta }^2+6\zeta +11\right){{{\delta }_f}^{\ast }}^2}}}\right)}$$E57

where

$$\xi =\frac{\left(x/H\right)}{\left(\frac{{\mathrm{\rho }}_f{\mathrm{c}}_{p_f}{u}_{B}H}{k_{\mathit{ref}}}\right)}$$E58

is the Graetz number. Moreover, $k_{\mathit{ref}}$ is a reference thermal conductivity.

For given properties of the graded porous materials, namely $\frac{\epsilon \left(x\right)}{\mathit{Da}\left(x\right)}$ and ${\mathit{Re}}_b\left(x\right)=\frac{\epsilon \left(x\right){\rho }_{f}b\left(x\right)H^2{u}_B}{{\mu }_f}$, the velocity profile parameter $\zeta \left(x\right)$ can be determined from Eq. (37). Once $\zeta \left(x\right)$ is known, the fluid phase boundary layer thickness ${{\delta }_f}^{\ast }$ can be found by using Eq. (57) as a function of $\xi =\frac{\left(x/H\right)}{\left(\frac{{\mathrm{\rho }}_f{\mathrm{c}}_{p_f}{u}_{B}H}{k_{\mathit{ref}}}\right)}$ for given local values of $\frac{k_{\mathit{fe}}\left(x\right)}{k_{\mathit{ref}}},\frac{k_{\mathit{se}}\left(x\right)}{k_{\mathit{ref}}},\mathit{Bi}\left(x\right)=\frac{h_v\left(x\right)H^2}{k_{\mathit{se}}\left(x\right)}$. With $\zeta \left(x\right)$ and ${{\delta }_f}^{\ast }\left(x\right)$ thus determined, the development of the dimensionless wall temperature $\frac{k_{\mathit{ref}}\left(T_w\left(x\right)-T_0\right)}{H{q}_w}$ of great interest may readily be found from

$$\frac{k_{\mathit{ref}}\left(T_w\left(x\right)-T_0\right)}{H{q}_w}=\frac{\frac{\left(\zeta +3\right)}{2\left(\zeta +2\right)}{{\delta }_f}^{\ast }}{\frac{k_{\mathit{fe}}\left(x\right)}{k_{\mathit{ref}}}+\frac{k_{\mathit{se}}\left(x\right)}{k_{\mathit{ref}}}\frac{2}{1+\sqrt{1+\frac{24}{\mathit{Bi}\left(x\right)}\frac{\left(\zeta +2\right)\left(\zeta +3\right)}{\left({\zeta }^2+6\zeta +11\right){{{\delta }_f}^{\ast }}^2}}}}$$E59

which is based on Eqs. (52) and (55). Correspondingly, the local Nusselt number may be calculated from

$$\begin{array}{rcl}{\mathit{Nu}}_H& \equiv & \frac{q_{w}H}{\left(T_w-T_{\mathit{fB}}\right)k_{\mathit{ref}}}=1/\left(\frac{k_{\mathit{ref}}\left(T_w\left(x\right)-T_0\right)}{H{q}_w}-\frac{k_{\mathit{ref}}\left(T_{\mathit{fB}}\left(x\right)-T_0\right)}{H{q}_w}\right)\\ & =& \frac{1}{\frac{\frac{\left(\zeta +3\right)}{2\left(\zeta +2\right)}{{\delta }_f}^{\ast }}{\frac{k_{\mathit{fe}}\left(x\right)}{k_{\mathit{ref}}}+\frac{k_{\mathit{se}}\left(x\right)}{k_{\mathit{ref}}}\frac{2}{1+\sqrt{1+\frac{24}{\mathit{Bi}\left(x\right)}\frac{\left(\zeta +2\right)\left(\zeta +3\right)}{\left({\zeta }^2+6\zeta +11\right){{{\delta }_f}^{\ast }}^2}}}}-\xi }\end{array}$$E60

3.2 Porosity variation and velocity profiles

As described in Eq. (61), the porosity of the graded metal foam in the channel of length L was assumed to vary axially according to the function as follows:

$$\epsilon \left(x\right)=1-\left(1-{\epsilon }_{\mathit{ref}}\right)\left(1+n\right){\left(\frac{x}{L}\right)}^n$$E61

such that the solidity remains constant irrespective of the value of the exponent n:

$$\frac{2H{\int }_0^L\left(1-\epsilon \left(x\right)\right)\mathit{dx}}{\left(2\mathit{HL}\right)}=\left(1-{\epsilon }_{\mathit{ref}}\right):\text{Solidity}$$E62

where ${\epsilon }_{\mathit{ref}}$ is the mean porosity of the channel. The porosity decreases (i.e. the solidity increases) steeply downstream as increasing the value of n from 0.

The correlations for the permeability $K$ and interstitial heat transfer coefficients $h_v$ can be found in Calmid [18], in which the hydrodynamic and thermal characteristics of flows in metal foams were investigated. A typical correlation for the permeability of metal foam was proposed in their publications:

$$\mathit{Da}\left(\epsilon ;\frac{d_m}{H}\right)=\frac{K}{H^2}=0.00073{\left(1-\epsilon \right)}^{-0.224}{\left(\frac{1.18}{1-{\mathrm{e}}^{-\left(1-\epsilon \right)/0.04}}\sqrt{\frac{1-\epsilon }{3\pi }}\right)}^{-1.11}{\left(\frac{d_m}{H}\right)}^2$$E63

In order to overcome the deficiency around the inlet in the Calmid correlation, a simple power law relationship is proposed, which is as follows:

$$\mathit{Da}={\mathit{Da}}_{\mathit{ref}}{\left(\frac{\epsilon }{{\epsilon }_{\mathit{ref}}}\right)}^m$$E64

Wong et al. [21] theoretically recommended the range of the exponent m in three-dimensional porous media. Furthermore, the power law (64) with m = 4 and ${\mathit{Da}}_{\mathit{ref}}=\mathit{Da}\left({\epsilon }_{\mathit{ref}};\frac{d_m}{H}\right)$ was expressed as

$$\mathit{Da}=\left(0.00073{\left(1-{\epsilon }_{\mathit{ref}}\right)}^{-0.224}{\left(\frac{1.18}{1-{\mathrm{e}}^{-\left(1-{\epsilon }_{\mathit{ref}}\right)/0.04}}\sqrt{\frac{1-\epsilon }{3\pi }}\right)}^{-1.11}{\left(\frac{d_m}{H}\right)}^2\right){\left(\frac{\epsilon }{{\epsilon }_{\mathit{ref}}}\right)}^4$$E65

gives the maximum value at the inlet and then roughly follows the urves based on the Calmid correlation, which continuously decreases downstream.

Figure 7a and b illustrate the axial development of the velocity field whose porosity varies following Eq. (61) with Da (x) evaluated from the power law (65) for the case of ${\epsilon }_{\mathit{ref}}=0.8$, $d_m/H=0.1,L/H=10,{\mathit{Re}}_H=1500,$ and n = 1.5.

Figure 7.

Axial development of velocity field. (a) Streamwise velocity component and (b) normal velocity components.

The streamwise velocity profiles $u/u_B$ clearly indicate that the viscous boundary layer becomes “thinner” downstream, as both porosity and permeability decrease downstream. The corresponding normal velocity component $-v/u_B$ at the same locations are presented in Figure 7b. Figure 7a and b together suggest that some part of the fluid mass is shifted from the core toward the heated wall. This weak secondary flow results in a comparatively higher streamwise velocity layer near the wall, which is favorable in view of cooling the heated wall.

3.3 Boundary layer thickness and wall temperature variation

It is found that the wall temperature and Nusselt number variations are quite insensitive to the difference in the local Darcy number observed in the above two correlations. This may be attributed to the forgetfulness of the inlet conditions, which is inherent to the boundary layer theory. Therefore, only the results based on the Calmid correlation (59) are presented.

A series of calculations were conducted analytically and numerically for forced convective flows in a channel filled with a graded metal foam under equal solidity, for the cases of ${\epsilon }_{\mathit{ref}}=0.8,d_m/H=0.1,L/H=10$, ${\mathit{Re}}_H=1500,\mathit{Pr}=0.7,$ and $\frac{k_f}{k_s}=0.01$. The exponent n for controlling the gradient of the porosity was set to 0, 0.5, 1.0, and 1.5, respectively. For given exponent n and other parameters, the velocity shape parameter $\zeta \left(x\right)$ is evaluated using Eq. (39) to determine the velocity field. The local value of ${{\delta }_f}^{\ast }\left(x\right)={\delta }_f/H$ also can be found by using a standard shooting method once the $\zeta \left(x\right)$ is locally determined, and then some other quantities of interest, such as the wall temperature, can readily be provided.

As illustrated in Figure 8, the thermal boundary layer thickness of the solid is always greater than that of the fluid phase due to the thermal conductivity of the solid phase is much greater. In fact, the thermal boundary layer of the fluid phase just started to grow at the inlet, namely ${\left.{\delta }_f\right|}_{x=0}=0$, while the boundary layer thickness of solid phase is non-zero as given by ${\delta }_{\mathrm{s}}=\sqrt{\frac{6}{\mathit{Bi}}\frac{\left(\zeta +2\right)\left(\zeta +3\right)}{\left({\zeta }^2+6\zeta +11\right)}}H$. Therefore, the boundary layer thickness ratio ${\delta }_f/{\delta }_s$ increases from zero at the inlet along downstream. Since ${\delta }_f/{\delta }_s=1$ is required for both fluid and solid temperature profiles to be identical, the L.T.E. assumption fails over a substantial distance from the inlet as can be seen from Figure 8.

Figure 8.

Axial development of thermal boundary layer thickness ratio.

The effects of the exponent n on the wall temperature development are illustrated in Figure 9.

Figure 9.

Axial development of wall temperature.

A typical thermal boundary layer development can be observed in Figure 9 for the case of uniform porosity n = 0, that is, under the condition of constant heat flux on the wall, the wall temperature increases from the inlet in proportion to $x^{\frac{1}{3}\sim \frac{1}{2}}$, then to $x$ far downstream, according to the classical boundary layer theory. The axial development of the wall temperature changes rather drastically depending on the value of n. When n is 0.5, the wall temperature near the inlet rises faster than the uniform porosity case, while the downstream wall temperature is suppressed. According to the correlation (61), the porosity for the case of n = 0.5 varies from 1 to 0.7 from the inlet to the exit based on the average porosity of 0.8, which changes in a moderate gradient but results in a remarkable temperature distribution variation over the axial distance. For the case of n = 1, the wall temperature was kept nearly at constant temperature of reasonably low level over most of the channel in axial direction of the channel. Thus, the results indicate that a uniform temperature can be achieved if the porosity is set to decrease linearly from the inlet to the outlet.

Furthermore, for the case of n = 1.5, the wall temperature steeply increases from the inlet and then decreases toward the exit after attaining its maximum value. The rather drastic increase in the effective thermal conductivity toward the exit explains the high sensitivity of the porosity gradient to the wall temperature distribution. Even a small increase in the metal volume fraction will result in a substantial increase in the effective thermal conductivity because the conductivity ratio of the metal to the fluid was set too large. Thus, the heat conducts away from the wall to the core region in the channel, keeping the wall temperature comparatively lower there.

3.4 Conclusions

In this section, the velocity profile and wall temperature development for the channel filled with axially graded metal foam were investigated analytically and numerically. The variation of the porosity follows power function of axial coordination with an exponent n based on an equal solid volume fraction. Velocity profiles for different cases were obtained based on the Brinkman-Forchheimer extended Darcy model. The solutions indicate that there exists a weak secondary flow toward the heated wall of the channel, enhancing cooling of the wall. It is found that the wall temperature is very sensitive to the axial gradient of the local porosity, even a moderate degree of the axial gradient with n = 0.5 changes the wall temperature distribution rather drastically. As a result, the wall temperature becomes almost uniform throughout the channel when n = 1.0, while, for n = 1.5, the wall temperature rises from the inlet, attains its maximum, and then decreases toward the exit. Thus, axially graded porous materials can be quite useful for designing effective cooling systems and controlling the desired working temperature and its uniformity.

4. Application of graded porous media in aeroengine cooling

The afterburner is an important component, especially in military aeroengine. It substantially improves the transient thrust of the aeroengine. However, the gas temperature at the afterburning state rises up to 2200 K, which is far higher than the temperature that the outer wall material can stand. Therefore, an efficient cooling system is definitely needed. The afterburner heat shield is often used to protect the outer wall from thermal load. However, the cooling area of the heat shield is large, film cooling holes are too many, but the coolant available for the cooling of the heat shield is relatively less. Therefore, the effective cooling of the heat shield is difficult, the high efficient cooling structure for the heat shield is definitely needed. However, due to the non-uniform pressure difference of the coolant and the mainstream, over the surface of the heat shield, the distribution of the coolant over the heat shield is also not uniform, which may lead to the thermal stress concentration on the heat shield itself. Therefore, the graded porous media may be adopted in the afterburner heat shield structure to improve its cooling effectiveness and temperature uniformity.

For convenience, in this section, four cases of impingement/effusion cooling structures embedded with uniform and non-uniform distribution pin fin array structures were considered as models to study the effects of graded structure on its heat transfer performance and cooling effectiveness. With several rib turbulator shapes, Rhee et al. [22] studied the motion and local heat/mass transfer parameters. They discovered that the location of the rib turbulators and the impingement/effusion holes affected the level of local heat/mass transfer, which was found to be boosted by the rib turbulators.

The cooling processes can be considerably impacted by lattice frame structures, which lead to higher levels of airflow disruption and an increase in the intensity of heat transfer within the flow channel thanks to the enlargement of the heat transfer surface area. Pins were shown by Funazaki [23] to be essential in improving the overall heat transfer of double-layer cooling. Li and Mongia [24] looked into how a double-walled cooling system affected heat transfer. Researchers suggested that for double wall cooling with a thin wall or thick wall, the addition of pins results in an increase of about 0.03 in overall cooling performance.

4.1 Calculation model

In this section, four cases of impingement/effusion cooling double-layer structures embedded with uniform or non-uniform (graded) distributed pin fin structures were treated as calculation models, as can be seen in Figure 10. Moreover, the relative location of impingement holes, effusion holes, and pin fins were also indicated.

Figure 10.

Calculation models: (a) case 1 (uniform distribution; 121 pins); (b) case 2 (almost uniform distribution; 96 pins); (c) case 3 (non-uniform distribution from highest density to lowest; 166 pins); (d) case 4 (non-uniform distribution from lowest density to highest; 166 pins).

It is significant to note that the shape of the hollow elements, as well as their diameter, remained constant. The location of the holes for both the upper and lower plates, as well as their diameters and angle of inclination, are also unchanged. Periodical boundary conditions are imposed on the side parts of the computational model. The CFD software was employed to create an unstructured mesh of the computational area, and a depiction of intricacy of the mesh is presented in Table 1.

Turbulence model	SST k-$\omega$
Mesh type	Three-dimensional
Mesh form	Orthogonal

Table 1.

The main parameters of the simulation.

Table 1 shows the main parameters of the simulation, and Table 2 illustrates the number of grid cells for each model.

Name of a case	Number of cells
Case 1 (uniform distribution)	1,809,286
Case 2 (almost uniform distribution)	1,543,251
Case 3 (non-uniform distribution from highest density to lowest)	2,304,339
Case 4 (non-uniform distribution from lowest density to highest)	2,294,669

Table 2.

The number of cells for each case.

4.2 Flow structure and temperature distribution

Figure 11a–d demonstrates the characteristic features of the flow structure for each model.

Figure 11.

Flow structures regarding the 4 cases with lattice frame elements of different distribution densities: (а) case 1, (b) case 2, (c) case 3 and (d) case 4.

In comparison to Figure 11a and b, it can be noticed that the most striking manifestations of vortices are revealed in the model whose structural elements are slightly deviated from a uniform distribution (Figure 11b). The difference between cases 3 and 4 (Figure 11c and d) is also noticeable. It is tangible that in the region of the lowest density of solid elements, the flow structure swirls throughout the volume of the occupied part of the space. However, otherwise, when the distribution density is high, there is a restriction of the free zone of flow propagation. In this case, good mixing of the flow is provided, which is accompanied by small vortex structures throughout the volume of the area that is occupied.

The evaluation of the flow behavior analysis is very necessary, as it affects further heat exchange processes. So, Figure 12 shows the temperature distributions for these 4 models.

Figure 12.

Temperature distribution regarding the different cases with different distribution densities of lattice elements: (a) case 1, (b) case 2, (c) case 3, and (d) case 4.

Going deeper into the analysis of the flow temperature inside the model in the interlayer space, it is noted that the part of the computational model where there is a high density of the distribution of elements of the lattice structure leads to a high intensity of heat exchange. This is due to the fact that with an increase in the distribution density, the interlayer space for mixing air flows decreases.

4.3 Overall cooling effectiveness

It is important to note that the overall cooling effectiveness is influenced by the interplay between different cooling mechanisms and their respective effectiveness values, and it is crucial to optimize the design and configuration of each cooling mechanism to achieve high overall cooling effectiveness. Figure 13 demonstrates the distribution of the overall cooling effectiveness on 3-D models. It is vital to note that the range of values ranges from 0.3 to 0.8 according to the indicator scale above with intervals of 0.325.

Figure 13.

3-D distribution of the overall cooling effectiveness for four models with different distribution densities of lattice elements: (a) case 1; (b) case 2; (c) case 3; (d) case 4.

From the scientific point of view, it can be noticed that with a decrease in the number of pins, the intensification of heat exchange decreases, while an increase in the same elements is accompanied by an increase in heat exchange processes. That is, the density of the distribution of lattice elements is directly proportional to the overall cooling effectiveness. This is the conclusion that can be traced from Figure 14.

Figure 14.

The overall cooling effectiveness for four models with different distribution densities of lattice elements: red line—case 1; green line—case 2; blue line—case 3; black line—case 4.

However, to assess the overall cooling effectiveness, a graph was also constructed, presented in Figure 14. The blue line characterizes model 3, that is, the one whose intensification was maximal already in the first segment of the x/d ratio and whose overall cooling effectiveness is 0.515. It is revealed that the decrease in the pin distribution density is directly proportional to the decrease in the overall cooling effectiveness.

The worst initial data was shown by model 4 (black line), where there was minimum distribution density at the first stage (0.47). The overall cooling effectiveness across the entire model ∆Ф in this case was 12%. Nevertheless, the black line at its minimum value began to grow sharply and reached a maximum of 0.71, while the difference in overall cooling effectiveness was as much as 25%. Such an effective indicator was achieved thanks to a rise in pins in the model under study, which only provoked enhancement in heat exchange intensification.

Summing up, the best indicators were identified in the case 3. Also, it boils down to the fact that it is necessary to introduce such lattice frame structures in which the tendency of the distribution of pins from a constructive point of view will be accompanied by an increase in the density of their distribution.

4.4 Conclusions

The simulation work of the graded lattice structures with different arrangements of pins was numerically investigated using the appropriate software (CFD). Four different cases with various distributions of the ligament elements were assumed, as well as the characteristics of the mesh construction were presented. The evaluation information showed that additional indicators that improve heat transfer processes are directly affected by the distribution density. In addition, a decrease in the distribution density leads to a decline in the overall cooling effectiveness. It is proved that when designing a model with a high density of the distribution of elements at the initial stage, there is a sharp increase in the efficiency of heat exchange. Thus, functionally graded lattice structures can be used to achieve significantly high overall cooling effectiveness while at the same time having a high distribution density of lattice structure elements.

5. Conclusions

In this section, the flow and heat transfer in channels filled with vertical graded porous media and longitudinal graded porous media have been analytically investigated, respectively. The heat transfer characteristics and wall temperature variation principles for both cases have been presented. According to the analytical results, the vertical graded porous media in a channel enables to change the local velocity in the regions of near the wall or in the core, prompts the flow field in the low porosity region more distorted while accelerates the flow filed in the high porosity region. The highest heat transfer coefficient is obtained in the case of where the porosity decreases toward the heated wall, namely, more metal near the wall and less metal in the core region. For the case of a channel filled with axially graded porous foam structure, which follows a power function of axial distance with an exponent n, however, the velocity field is found that exist a weak secondary flow toward the heated wall, enhancing cooling of the wall. It is found that the exponent n controlling the axial gradient of the local porosity has a strong influence on the axial development of the wall temperature. A uniform wall temperature can be achieved in the case of n = 1.0, while a maximum temperature can be attained at the middle of the wall in the case of n = 1.5, instead of the fact that temperature keeps rising from the inlet to the end under the constant heat flux, which all the other researchers hold. The finding can be useful for designing effective cooling system, controlling the desired working temperature, and its uniformity.

On the other hand, the possibility of application of graded porous media in aeroengine cooling has been discussed. The numerical simulation work on the cooling of the afterburner heat shield with graded pin fin structures has been conducted. Four different cases of double-layer heat shields with various distributions of the ligament elements were constructed and examined their flow structure and cooling effectiveness. As a result, the decrease in the distribution density leads to decline in the overall cooling effectiveness, whereas a high density of pin fin elements arranged at the initial stage prompts the heat exchange and improves the overall cooling effectiveness upstream of the heat shield, although sacrifices part of that at the downstream, the gradient arrangement of pin fins improves the temperature uniformity over the surface of the heat shield. Thus, functionally graded lattice structures can be used to achieve significantly high overall cooling effectiveness while at the same time having a high distribution density of lattice structure elements.

Acknowledgments

The authors acknowledge gratefully the support from the National Science and Technology Major Project (J2019-III-0019-0063), the Innovation Capacity Support Plan in Shaanxi Province of China (Grant No.2023-CX-TD-19), the Grant of National Science Foundation of China (No.52006179), and the Fundamental Research Funds for the Central Universities (No.501XTCX2023146001).

Nomenclature

h	heat transfer coefficient
H	half channel height
k	thermal conductivity
ke	effective stagnant thermal conductivity
K	local permeability
NuH	Nusselt number NuH=hHkeεB
p	pressure
qw	wall heat flux
T	temperature
T∗	dimensionless temperatureT∗=T−TwkeεB/Hqw
u	Darcian velocity (apparent velocity) in x direction
u∗	dimensionless velocityu∗=u/−H2μdpdx
uB	bulk mean velocity
x	coordinate in the flow direction
y	coordinate normal to the channel wall
ε	porosity
ε∗	effective porosity
εB	bulk mean porosity
η	dimensionless vertical coordinate η=y/H
μ	viscosity
B	bulk mean
c	channel center-line
f	fluid phase, channel without a metal foam matrix
s	solid phase
w	wall