Abstract

The local thermal nonequilibrium (LTNE) model has been used widely for analyzing heat transfer during internal flow through porous media, including when a channel is only partially filled with a porous medium. In such problems, the Biot number describes the rate of convective heat transfer between solid and fluid phases. While uniform Biot number models are commonly available, recent advances in functionally graded materials necessitate the analysis of spatially varying Biot number in such geometries. This paper presents LTNE-based heat transfer analysis for fully developed flow in a channel partially filled with porous medium and with spatially varying Biot number to describe solid-fluid interactions in the porous medium. Fully uncoupled ordinary differential equations for solid and fluid temperature distributions are derived under three different boundary condition models. Solid and fluid temperature fields are presented for a variety of Biot number distributions, including quadratically and periodically varying functions. An explanation of the nature of temperature distribution predictions for such problems is provided. For special cases, the results presented here are shown to reduce to past work on constant Biot number. This work improves the theoretical understanding of porous media heat transfer and facilitates the use of such theoretical models for functionally graded materials.

1 Introduction

Fluid flow and heat transfer in porous media are of much practical and theoretical interest, with applications including micro-electronics cooling [1,2], geothermal energy extraction [3], thermofluidic/thermochemical energy systems [4], oil/gas extraction [5], bio-engineering [6] and thermal energy storage [7]. Heat transfer and fluid flow in porous media are relevant for both naturally occurring materials, such as soil [8], as well as engineered materials such as metal foams [2]. In both cases, the porous nature of the material results in increased specific interfacial area for heat transfer and facilitates performance unavailable with nonporous materials.

Several complications exist in the theoretical analysis of porous media heat transfer. Specifically, modeling of local thermal interactions between solid and fluid phases that co-exist in the porous medium must be carefully developed. The simplest approach is to ignore local heat exchange between solid and fluid phases by assuming local thermal equilibrium between the two phases [9]. This may not be valid beyond a limited number of specific, simplified problems. In a more general approach, local thermal nonequilibrium (LTNE) between the two phases may be assumed [10], with distinct temperature fields for the solid and fluid phases linked to each other through a local convective heat transfer coefficient. Due to its improved accuracy, the LTNE technique has been used widely to study heat transfer characteristics in a channel filled with porous medium [1013]. For instance, the LTNE technique was used to investigate constant wall heat flux boundary conditions in porous media [14] and to model two-phase flows with phase-change in porous media [15]. Partially porous flow in a pipe [16] and between parallel plates [17] has been analyzed. Solutions for conjugate heat transfer problems under the LTNE assumptions have been derived [18]. In addition, a numerical and experimental study of liquid-cooled aluminum foam to cool power electronics was reported recently [19]. The study validated the enhancement of heat transfer as a result of using a porous medium and the accuracy of LTNE over local thermal equilibrium [19]. In addition, an experimental investigation to highlight the heat transfer performance of a channel filled with metal foam has also been carried out [2].

In general, a porous medium within a channel aids in heat transfer, but also results in increased pressure drop. A possible mitigation strategy is to fill the channel only partially with a porous medium, resulting in fully developed fluid flow in the remainder of the channel. This approach has been studied in several theoretical and experimental papers. An experimental and numerical study on heat transfer enhancement for gas heat exchangers fitted with porous medium confirms the consequence of increased pressure drop and increased pumping power [20]. Numerical analysis of a partially filled porous channel has been presented [21]. In the modeling of a partially filled channel, the boundary condition at the porous-fluid interface must be modeled carefully. Several models have been proposed for this purpose. In the first model, the temperatures of solid and fluid phases are assumed to be equal at the interface because of a very large rate of heat transfer (model A) [12]. The second model in turn assumes different temperatures for the solid and fluid phases at the interface and distribution of heat flux between both the phases, through an interface thermal parameter (model B) [12]. Finally, a flux jump condition has been proposed based on an interface heat transfer coefficient (model C) [12,22].

Most of the past work on porous media heat transfer, including in partially filled channel, assumes a uniform Biot number (Bi) throughout the channel [1114,23], which is defined based on of the local convective heat transfer coefficient, the specific interfacial area in the porous medium, the effective thermal conductivity of the solid phase and the half-width of the channel. While this may be a reasonable assumption for several applications, a nonuniform Bi may be encountered in specific cases. For example, change in the specific interfacial area across the channel may result in a nonuniform Bi. This could occur in a functionally graded material [24] or due to variations caused by inadvertent manufacturing defects. A numerical investigation has shown that nonlinear gradient in porous material may help improve heat transfer without pressure drop penalty [25]. Given the relative lack of literature in this direction, however, more detailed theoretical modeling, including analysis of the fluid-porous interface is needed for a partially porous channel, in which Bi varies across the channel. While an LTNE-based theoretical model has been presented to account for spatial variation in Bi in a fully porous channel [10], a similar analysis for a partially filled channel is desirable.

This paper presents an LTNE-based analysis of heat transfer in a channel partially filled with porous medium with spatially varying Biot number. In this work, spatially varying Biot number is modeled as a consequence of variation in the specific interfacial area. Uncoupled fourth-order ordinary differential equations for solid and fluid temperature fields are derived for three different interface heat transfer models. Solutions of the uncoupled ODEs are used to understand the effect of Biot number variation and other nondimensional parameters on the temperature distribution of the solid and fluid phases in the channel. Results presented here improve the fundamental understanding of porous media heat transfer, and may also help in the design and optimization of functionally graded materials in heat transfer applications.

2 Mathematical Modeling

Figure 1 presents a schematic of the problem under consideration. A symmetric channel with half-width H is partially filled with a porous medium of half-width H1. The channel experiences a uniform heat flux qw at the wall. Steady, incompressible, and laminar fluid flow in the porous region governed by the Darcy model is assumed across the channel, with thermally and hydrodynamically developed conditions. Assumption of local thermal nonequilibrium between the solid and fluid phases, along with heat generation in the porous region of the channel results in the following governing energy equations [12]:
kf,eff2Tfy2+hintα(TsTf)+Sf=ρCpupTfx
(1)
ks,eff2Tsy2hintα(TsTf)+Ss=0
(2)

Here, Tf and Ts are fluid and solid phase temperatures in the porous medium, respectively, whereas up, Sf, Ss, kf,eff, ks,eff, ρ, Cp, and hint are the porous region fluid velocity, fluid and solid internal heat generation, effective fluid and solid thermal conductivities, fluid density, fluid specific heat, and interstitial heat transfer coefficient, respectively. α is the specific interfacial area in the porous medium. Hydrodynamic analysis of this problem has been presented before [12] and is not affected by the spatially varying Bi considered in this work. These hydrodynamic equations are summarized below.

Using the Darcy model, the constant velocity of the fluid flowing in the porous region is [26]
up=dPdxKμf
(3)

where K is the permeability of the porous medium, and μf and P are the dynamic viscosity of the fluid and pressure, respectively.

The corresponding governing equations for fluid flow in the open region of the channel are as follows [12]:
μfd2uody2=dPdx
(4)
kf2Tfy2=ρCpuoTfx
(5)

In Eqs. (4) and (5), u0 is the velocity of the fluid in the open region.

The associated boundary conditions are as follows:
(uoy)y=0=0
(6)
(uo)y=H=0
(7)
(uy)y=H1+=α*K(uiup,i)
(8)
(Tfy)y=0=0
(9)
(Tsy)y=0=0
(10)
kf(Tfy)y=H=qw
(11)

Equations (6) and (9) arise from the assumed symmetry in the problem. Equation (7) represents the no-slip condition at the channel wall. Subsequently, the interface condition represented by Eq. (8) corresponds to the slip velocity condition [27], which implies relative motion between the fluid in the open region and porous medium interface. In Eq. (8), ui is the interface pure fluid velocity, up,i is the fluid velocity in the porous region at the interface and α* is the velocity slip coefficient. The velocity slip coefficient is a dimensionless quantity that mainly depends on the structure of the porous material near the interface [27]. Velocity slip coefficient impacts the interface velocity, as well as the rate of heat transfer in the transition region. Further, the slip velocity coefficient is known to have a strong dependence on the porosity and permeability of the porous medium [28].

While Eq. (8) fully defines the hydrodynamic interfacial condition at the interface between porous and nonporous regions, a similar thermal interfacial condition remains to be defined. In the past work, three distinct interface models have been proposed [12]. The assumption of a very high rate of heat transfer between the solid and fluid results in model A, in which, the two phases are at the same temperature at the interface. The assumption of a moderate rate of heat transfer between the solid and fluid gives rise to model B, in which, heat flux is distributed between the solid and fluid phases by an interface thermal parameter. Finally, by assuming an interface heat transfer coefficient to determine the rate of heat transfer between the solid and fluid, one may obtain model C. Nondimensionalization is carried out as shown below:
θ=ks,eff(TTs,i)qwH;η=yH;η1=H1H;k=kf,effks,eff;k1=kfks,eff;Bi(η)=hintαH2ks,effLs=SsH1qi;Lf=SfH1qi;Bii=hiHks,eff;Da=KH2;U=udPdxH2μf;γ=qiqw
(12)

where qi is the heat flux at the interface. Note that Bi in this work is assumed to be a function of η.

The resulting equations and solutions are discussed below for three specific interfacial models next.

2.1 Model A.

Model A assumes that temperatures of the two phases are equal at the interface due to a sufficiently large rate of heat transfer between the solid and fluid at the interface. The corresponding boundary and interface conditions for this model are as follows:
(Tfy)y=0=(Tsy)y=0=0
(13)
(Tf)y=H1=(Ts)y=H1=(Tf)y=H1+
(14)
kf,eff(Tfy)y=H1+ks,eff(Tsy)y=H1=kf(Tfy)y=H1+=qi
(15)

Here, H1 and H1+refer to the porous and pure-fluid sides of the interface, respectively.

Expressions for the advection terms in Eqs. (1) and (5) are given as shown below [12]:
ρCpupTfx=qiH1+Ss+Sf
(16)
ρCpum,oTfx=qwqiHH1
(17)
Using the nondimensionalization scheme described in Eq. (12), the expressions for the fluid velocity in the porous medium and open region are as shown below [12]:
Up=Da0ηη1
(18)
Uo=(ηη1)22+α*Da(UiDa)(ηη1)+Uiη1<η1
(19)
where Ui is the dimensionless interface pure fluid velocity and is given by
Ui=(1η1)22+α*Da(1η1)1+α*Da(1η1)
(20)
The dimensionless average velocity over the entire channel is as follows:
Um=η1Up+(1η1)Um,o
(21)
In Eq. (21), Um,o is the dimensionless average velocity of the fluid in the open region, given by
Um,o=(1η1)26+α*2Da(UiDa)(1η1)+Ui
(22)
Now, dividing Eq. (17) with Eq. (16) and carrying out nondimensionalization, along with the use of Eq. (21) results in an expression for the dimensionless heat flux at the interface
γ=qiqw=1(1+Lf+Ls)(UmUpη11)+1
(23)
Note that setting Lf = Ls = 0 in Eq. (23) results in the same expression for the dimensionless interfacial heat flux as Yang and Vafai's work [12], in which, heat generation was neglected. The following set of nondimensional equations and boundary/interface conditions are obtained from Eqs. (1)(2), (5), (11), and (13)(15):
kθf(η)+Bi(η)(θS(η)θf(η))=γη1(1+Ls)0ηη1
(24)
θs(η)Bi(η)(θS(η)θf(η))+Ls=00ηη1
(25)
k1θf(η)=Uo(UmUpη1)+(Upη11+Lf+Ls)η1<η1
(26)
θf(0)=θs(0)=0
(27)
θf(η1)=θs(η1)=θf(η1+)=0
(28)
kθf(η1)+θs(η1)=k1θf(η1+)=γ
(29)
θf(1)=1k1
(30)

It can be noted that the governing equation for the open region differs slightly from previous work [12], where heat generation was not assumed. Setting Ls and Lf equal to zero in Eq. (26) reduces it to the form appropriate for zero heat generation [12].

θs and θf are coupled to each other in Eqs. (24) and (25). In order to uncouple, θs is written in terms of θf from Eq. (24) and substituted in Eq. (25), resulting in
θf(kBi)+θf(1+k+2k(Bi)2Bi3kBiBi2)+Ls(1γη1+2γ(Bi)2η1Bi3γBiη1Bi2)+(2γ(Bi)2η1Bi3γBiη1Bi2γη1)=0
(31)
Similarly, an uncoupled ordinary differential equation for θs is obtained as follows:
θs(kBi)+θs(2kBiBi2)+θs(1+k+kBiBi22k(Bi)2Bi3)+Ls(1γη12k(Bi)2Bi3+kBiBi2)γη1=0
(32)

Note that unlike past work, Bi is treated to be a function of η here, which is why Eqs. (31) and (32) are more complicated than similar results for constant Bi in past work [12].

Additional boundary conditions needed to solve Eqs. (31) and (32) are obtained as follows: First, differentiating Eqs. (24) and (25) and evaluating at the center of the channel results in
θf(0)=Bi(0)(θs(0)θf(0))k
(33)
θs(0)=Bi(0)(θs(0)θf(0))
(34)
Further, Bi(0)=0 due to channel symmetry, and, therefore, Eqs. (33) and (34) result in
θf(0)=0
(35)
θs(0)=0
(36)
The remaining boundary conditions are obtained by simply evaluating Eqs. (24) and (25) at the interface and using Eq. (28)
θf(η1)=γη1k(1+Ls)
(37)
θs(η1)=Ls
(38)
Temperature distribution in the nonporous region remains unaffected by the Bi generalization carried out in this work. The temperature distribution in the nonporous region is [12]
θf=D0(ηη1)4+D1(ηη1)3+D2(ηη1)2+D3(ηη1)
(39)
where
K=(UmUpη1)+Upη11+Lf+Ls;D0=124Kk1;D1=α*6Kk1Da(UBDa);D2=Ui2Kk1;D3=1k14D0(1η1)33D1(1η1)22D2(1η1)
(40)

In summary, Eqs. (31) and (32) are the uncoupled ODEs for θf and θs, respectively, that must be solved in order to determine solid and fluid temperature distributions in the porous region in the presence of spatially varying Biot number. Associated boundary conditions for θf are given by Eqs. (27), (28), (35), and (37), and for θs are given by Eqs. (27), (28), (36), and (38). While a general analytical solution is not likely to be possible for the general case of Bi(η), these equations can be easily solved numerically.

2.2 Model B.

Model B assumes that heat flux distribution in the solid and fluid phases of the porous medium at the interface is proportional to the porosity of the porous medium [12]. The governing energy equations for model B are the same as model A, given by Eqs. (24)(26). The boundary condition given by Eq. (27) remains unchanged. Other associated boundary and interface conditions are as follows:
θf(η1)=θf(η1+)andθs(η1)=0
(41)
kθf(η1)=βγ
(42)
θs(η1)=(1β)γ
(43)
k1θf(η1+)=γ
(44)
The uncoupled fourth-order ordinary differential equations for model B are the same as for model A, given by Eqs. (31) and (32), since these models differ only in the treatment of the interfacial boundary condition. Further, boundary conditions at the center of the channel are also the same, given by Eqs. (35) and (36). Other boundary conditions for model B are derived. Evaluating Eq. (24) at the interface results in
kθf(η1)Bi(η1)θf(η1)=γη1(1+Ls)
(45)
Now, differentiating Eq. (25) and evaluating at the interface gives the final boundary condition as follows:
θs(η1)Bi(η1)(θs(η1)βγk)+Bi(η1)θf(η1)=0
(46)
The temperature distribution in the open region for this model can be obtained from the expression shown below [12]:
θf=D0(ηη1)4+D1(ηη1)3+D2(ηη1)2+D3(ηη1)+θf(η1)
(47)

where D0, D1, D2, and D3 are given by Eq. (40). The uncoupled ODEs for θf and θs, given by Eqs. (31) and (32) are solved numerically using corresponding boundary conditions. The boundary conditions for θf are given by Eqs. (27), (42), (35), and (45), whereas the associated boundary conditions for θs are given by Eqs. (27), (41), (36), and (46).

2.3 Model C.

The final model considered here assumes an interface heat transfer coefficient hint [12] to distribute the heat flux among the solid and fluid phases (model C). Similar to model B, temperatures of the solid and fluid phases at the interface are not equal, and are, in fact, governed by the value of hint.

Similar to model B, Eqs. (24)(26) remain the governing differential equations for the fluid and solid temperature distributions, and as usual, the nondimensionalization scheme helps transform the boundary/interface conditions into the following.
kθf(η1)=γBii(θf(η1)θs(η1))
(48)
θs(η1)=Bii(θf(η1)θs(η1))
(49)
where
Bii=hiHks,eff
(50)
The boundary and interface are conditions given by Eqs. (27), (41), and (44) remain the same.
The uncoupled ordinary differential equations and the set of additional boundary conditions do not change for this model as well. Similar to model B, the remaining boundary/interface conditions are obtained as follows:
kθf(η1)Bi(η1)θf(η1)=γη1(1+Ls)
(51)
θs(η1)Bi(η1)(θs(η1)γk+Biiθf(η1))+Bi(η1)θf(η1)=0
(52)

In this case, the relevant boundary conditions for θf are given by Eqs. (27), (48), (35), and (51) and for θs are given by Eqs. (27), (41), (36), and (52). Finally, the temperature distribution for the open region can be determined using the expression available in the previous section.

3 Results and Discussion

3.1 Comparison With Past Work.

A fourth-order boundary value problem solver is used to solve the governing equations associated with models A, B, and C using appropriate boundary conditions. The three-stage Lobatto IIIa formula is used by the solver. Results from the present model are first compared with past work for the special case of uniform Bi. Figure 2 compares fluid and solid temperature distributions predicted by this work with results from Yang and Vafai [12] for a special case of constant Biot number, Bi =10. Results are plotted for Models A and C in Figs. 2(a) and 2(b), respectively. Table 1 lists the values of various parameters for this and all subsequent figures. In this case, k = 0.01, η1 = 0.98, Bi = 10, Ls = 0, Lf = 0, α* = 0.1, and Da = 10−5 for model A. In addition, Bii = 0.1 for model C. Note that η1 is very close to 1, indicating that nearly the entire channel is porous. Figures 2(a) and 2(b) both show excellent agreement between the present analysis and past work. The solid phase temperature at the interface η = η1 is shown to be zero for both models A and C, which is consistent with the boundary condition given by Eqs. (28) and (41), respectively. In addition, for model A, the solid and fluid temperatures match at the interface, which is consistent with model assumptions discussed in Sec. 2.1. On the other hand, the curves do not necessarily match with each other at the interface for model C, as shown in Fig. 2(b). Finally, for both models, the curves exhibit zero slope at the center of the channel, which is justifiable by the symmetry conditions assumed in the problem.

Further, a special case of this work is compared with a recent paper that reported LTNE modeling for spatially varying Biot number in a fully filled channel [10]. For k =1, Bi = Bi0(1−η2), Da =10−5, Bi0 = 30, α* = 0.1, Ls = 0, and Lf = 0, Figs. 3(a) and 3(b) plot fluid and solid temperature distributions for this work with different values of η1. For comparison, results corresponding to the fully porous channel case for the same set of parameters reported in a recent paper [10] are also plotted. As η1 becomes closer and closer to 1, the fluid and solid temperature distributions predicted by this work get closer and closer to results for the previously reported fully porous channel case. At η1 = 1, the two are identical, thereby providing additional verification of this work for a specific special case. Similar to Fig. 2, the results plotted in Fig. 3 are consistent with the various boundary conditions considered in the problem.

While the model is, in general, capable of accounting for any Biot number distribution, two specific Biot number distributions are of particular interest and, therefore, are analyzed next. Polynomial and sinusoidal functions are commonly used to approximate and represent more general functions encountered in practical applications. Therefore, the effect of Bi represented by such functions is analyzed next.

3.2 Analysis of Quadratically Varying Biot Number Distribution.

Figure 4 plots fluid and solid temperature distributions for a quadratically varying Biot number, given by Bi = Bi0 (1 + η2). The effect of the magnitude of the Biot number function, Bi0 is investigated first. Other problem parameters for this figure are k =1, η1 = 0.98, β = 0.9, Bii = 1, Ls = 0, Lf = 0, α* = 0.1, and Da = 10−5. Results for models A, B, and C are plotted in Figs. 4(a)4(c), respectively. In each case, the fluid temperature distributions approach those of the solid close to the interface. Due to the nature of Biot number distribution considered here, the rate of local solid-to-fluid heat transfer increases going from the centerline to the interface. In addition, for model A, the boundary condition also causes the fluid and solid temperatures to converge with each other at the interface, as seen in Fig. 4(a). In each model, the temperature of the solid phase increases as η increases, which is attributable to the heat flux imposed at the interface η = η1. A similar increase in fluid temperature as η increases is also attributed to heat flux at the interface as well as localized solid-to-liquid heat transfer that increases with η. In the case of model C (Fig. 4(c)), the quadratic Biot number near the interface is dominant over the interface Biot number that defines the amount of heat transferred from the fluid to the solid phase in the porous region. Therefore, the rate of heat transfer to the solid phase closer to the interface is much lesser than the rate of heat transfer to the fluid phase, which is the reason behind the greater fluid temperature near the interface. For model B, most of the heat flux at the interface flows through the fluid phase because of the large value of β. Therefore, the temperature of the fluid phase near the interface is higher than that of the solid phase for Bi0 = 0.5 and Bi0 = 0.7 as a result of increased heat transfer to the fluid phase. In addition, the quadratically increasing Biot number also contributes toward the increased heat transfer to the fluid phase. In contrast, for the Bi0 = 0.3 cases the Biot number is not as high as the other cases. Hence, the temperature of the fluid phase is lower than the solid phase temperature even near the interface.

3.3 Analysis of Thermal Conductivity Ratio for Quadratically Varying Biot Number Distribution.

Figure 5 investigates the effect of the ratio of thermal conductivities, k on the temperature distribution. In contrast with k =1 in Fig. 4, k is taken to be 10 in Fig. 5. All other parameters, including the quadratic Biot number distribution, are the same as Fig. 4. The impact of the greater thermal conductivity ratio is clearly seen in the resulting solid and fluid temperature plots for models A, B, and C plotted in Figs. 5(a)5(c), respectively. Similar to Fig. 4(a), the fluid and solid temperatures come close to each other as η increases in Fig. 5(a) due to the boundary condition. In model A, the effective thermal conductivity determines the amount of heat flowing through each of the phases. As a result of the greater value of k in Fig. 5(a), the effective thermal conductivity of the fluid phase is much higher than that of the solid phase. This causes higher Bi0 values to have a diminished impact on the fluid temperature throughout the channel. Interestingly, this reduces the solid phase temperature throughout the channel. In the case of model B, the solid phase temperature distribution is largely invariant of the value of Bi0, whereas the fluid temperature distribution goes up with increasing Bi0. Both temperature distributions exhibit weak dependence on η because of the high value of fluid thermal conductivity and β. In comparison with Fig. 4(c), the behavior of the plot for model C largely remains the same with only a change in the scale. The behavior does not change because the localized heat transfer still has the same effect on the temperature distributions.

3.4 Analysis of Periodical Biot Number Distribution.

Figure 6 presents the temperature distribution plot for a periodically varying Biot number, given by Bi = 50(1 + cos(2πωη)). Other problem parameters are k =1, η1 = 0.98, β = 0.9, Bii = 1, Ls = 5, Lf = 0, α* = 0.1, and Da = 10−5. Temperature distributions are plotted for models A, B, and C in Figs. 6(a)6(c), respectively. In each case, three different frequencies, ω = 1, 2, and 3 are considered. In each model, and for different values of ω, the Biot number is maximum at the center of the channel and very close to the maximum value at the porous-fluid interface. As expected, the temperature distributions appear close to each other at η = 0. In the case of model A alone, the temperatures of the solid and fluid phases converge at the interface because of the boundary condition. For models B and C, the fluid phase temperature curves near the interface are located above the solid phase temperature curves. For model B, this is due to increased heat transfer to the fluid phase at the interface. In the case of model C, the dominance of the periodic Biot number over the interfacial Biot number near the interface results in greater fluid temperature near the interface. In the rest of the porous medium away from the interface, the solid phase temperature is significantly higher because of heat generation in the solid phase. For ω = 2, Biot number is minimum at η = 0.25 and 0.75. This can be noticed through an increase in the difference between the temperature distributions at these points, with the fluid temperature staying below the solid temperature. For ω = 3, the Biot number has two additional maxima, in addition to the ones at the center of the channel and the interface. At these points, a reduction in temperature difference can be observed, with the fluid phase distribution lying below the solid phase curve.

3.5 Analysis of Thermal Conductivity Ratio for Periodic Biot Number Distribution.

Figure 7 presents results for a similar periodic Biot number distribution, but with a greater thermal conductivity ratio, k =10. Similar to Fig. 6(a), the temperature distributions converge at the interface for model A. For model B, the increased heat transfer to the fluid phase as seen in the previous case (Fig. 6) causes the temperature of the fluid phase to be higher than the solid phase temperature near the interface. In the case of model C, the Biot number still dominates over the interfacial Biot number, which causes the temperature of the fluid phase to be greater than that of the solid phase near the interface. In each model, the effect of the Biot number on fluid temperature variation is negligible compared to the effect on solid temperature distribution. This is attributable to the relatively higher fluid thermal conductivity in this case due to the larger value of k.

3.6 Analysis of Porous Region Thickness.

Temperature distributions predicted by the three models for different thicknesses of the porous medium, i.e., η1, are plotted in Figs. 8 and 9, for quadratically varying (Bi = Bi0 (1 + η2)) and periodically varying (Bi = 50(1 + cos(2πωη)) Biot number distributions, respectively. The problem parameters used for these figures are k =10, Da =10−5, Bii = 1, Bi0 = 0.5, Ls = 0, Lf = 0, α* = 0.1, and ω = 1. As expected, for each model, the temperature difference between the solid and fluid phases reduces as η1 decreases, as shown in Fig. 8. In Fig. 9, at the location of the minimum Biot number, the temperature difference between the solid and fluid phase is the lowest for η1 = 0.96. This is due to a reduction in the imposed heat flux at the interface as η1 is reduced. As seen in Fig. 8, for the quadratically varying Biot number case, reducing η1 has a significant effect on the fluid phase temperature distribution. On the contrary, for the periodic Biot number case, the effect is uniform across both phases. Beyond a certain value of η1, the temperature difference will reduce to zero and the temperature distributions in the porous region will be a constant zero because the imposed heat flux at the interface will be negligible.

3.7 Analysis of Interfacial Biot Number.

The interfacial Biot number, Bii is a key nondimensional parameter governing model C. Bii defines the rate of heat transfer from the pure fluid region to the solid phase at the interface, η = η1, as shown in Eq. (52). Figure 10 plots the temperature distributions for model C to analyze the effect of Bii. The interfacial Biot number essentially competes with the local Biot number to determine the overall rate and direction of heat transfer. For lower values of Bii, the quadratically increasing Biot number dominates near the interface. Therefore, the rate of heat transfer to the solid phase closer to the interface is much lesser than the rate of heat transfer to the fluid phase, thereby raising the fluid phase temperature beyond the solid phase temperature. As the interfacial Biot number is increased, this effect is countered, and the fluid phase curve shifts to the right. This can be noticed in Fig. 10. For a particular interfacial Biot number that is sufficiently large, the temperature of the fluid phase is equal to the temperature of the solid phase at the interface.

3.8 Analysis of Internal Heat Generation in the Solid Phase.

The impact of heat generation is also investigated. Heat generation is more likely to occur in the solid phase, for example, due to Joule heating. Figure 11 presents the effect of Ls, the solid phase internal heat generation on the temperature distribution for the case of a periodic Biot number. The problem parameters used for this figure are k =1, η1 = 0.98, β = 0.9, Bii = 1, ω = 2, Lf  =  0, α* = 0.1, and Da = 10−5. Three different values for Ls are considered – 0, 3, and 5. As expected, in each model, the solid and fluid phase temperatures in the porous medium go up as the internal heat generation increases. The difference between solid and fluid temperatures also increases, as increasing the value of Ls has a more significant impact on the solid temperature. As expected, the temperature difference is highest at points where the Biot number is minimum.

3.9 Analysis of Velocity Slip Coefficient.

In Fig. 12, the effect of velocity slip coefficient on the temperature distribution is analyzed using model B for a quadratically varying (Bi = Bi0 (1 + η2)) Biot number distribution. Problem parameters used to analyze the effect of the velocity slip coefficient are k =1, η1 = 0.98, β = 0.9, Bi0 = 0.5, Ls = 0, Lf = 0, and Da = 10−5. As α* increases, the velocity of the fluid at the interface reduces. This causes decreased advection in the transition region. As a result, the temperature of the fluid phase is higher for α* = 0.3 near the interface. At the same time, the temperature of the solid phase and fluid phase in the porous region is lower for α* = 0.3 because of lesser heat flow into the porous region.

3.10 Analysis of Darcy Number.

Finally, temperature distributions for model B are plotted in Fig. 13 for different values of the Darcy number. Quadratically varying Biot number distribution (Bi = Bi0 (1 + η2)) is used for this figure. Other parameters used to study the impact of Darcy number are k =1, η1 = 0.98, β = 0.9, Bi0= 0.5, Ls = 0, Lf = 0, and α* = 0.1. As the Darcy number reduces, the temperature of the solid phase starts to rise due to the diminishing effect of advection in the porous region. If the Darcy number is reduced beyond a certain point, the effect of advection in the porous region is negligible. This results in the solid and fluid phase distributions aligning with the solid phase temperature at the interface. In the figure, this can be noticed when Da =10−10.

4 Conclusions

Heat transfer during internal flow through porous media has been widely analyzed using the LTNE assumptions. While much of the past work assumed spatially invariant Biot number, recent progress in functionally graded porous materials has necessitated theoretical analysis of problems with spatially varying Biot number, which may occur, for example, due to variation in the specific interfacial area of the porous medium. The analysis presented in this work in the context of a partially porous channel provides the theoretical groundwork for analyzing heat transfer problems in functionally graded porous materials.

While results are derived here for general Biot number functions, two specific functions—quadratic and sinusoidal—are discussed in more detail due to their practical importance. The resulting solid and fluid temperature curves for three different interfacial models are consistent with constraints imposed on the problem by various boundary conditions. It is expected that the results presented in this work improve the theoretical understanding of heat transfer in a porous medium, and may also help in the design and optimization of functionally graded porous materials in a wide variety of applications.

Funding Data

  • National Science Foundation (CAREER Award No. CBET-1554183; Funder ID: 10.13039/100000001).

Nomenclature

Bi =

Biot number in the porous region, Bi(η)=hintαH2ks,eff

Bii =

interfacial Biot number, Bii=hiHks,eff

Cp =

specific heat of the fluid (J kg1 K1)

Da =

Darcy number, Da=KH2

H =

half height of the channel (m)

hi =

interfacial heat transfer coefficient at the interface between porous and nonporous regions (W m2 K1)

hint =

interstitial solid–liquid heat transfer coefficient in the porous region (W m2 K1)

H1 =

half height of the porous medium (m)

k =

thermal conductivity ratio, k=kf,effks,eff

K =

permeability (m2)

keff =

effective thermal conductivity (W m1 K1)

kf =

thermal conductivity of the fluid (W m−1 K−1)

k1 =

ratio of fluid thermal conductivity to solid effective thermal conductivity, k1=kfks,eff

L =

nondimensionalinternalheatgeneration,L=SH1qi

p =

pressure (N m−2)

qi =

heat flux at the interface (W m−2)

qw =

heat flux at the wall (W m2)

S =

internal heat generation (W m−3)

T =

temperature (K)

u =

fluid velocity (m s1)

U =

dimensionless fluid velocity, U=udPdx·H2μf

ui =

interface pure-fluid velocity (m s−1)

Ui =

dimensionless interface pure fluid velocity

Um =

dimensionless average fluid velocity over the entire channel

x =

horizontal coordinate (m)

y =

vertical coordinate (m)

α =

specific interfacial area of the porous medium (m−1)

α* =

velocity slip coefficient

β =

ratio of fluid phase heat flux to the total heat flux at the interface

γ =

dimensionless interfacial heat flux, γ=qiqw

η =

nondimensional coordinate, η=yH

η1 =

nondimensional thickness of porous region, η1=H1H

θ =

nondimensional temperature, θ=ks,eff(TTs,i)qwH

μ =

dynamic viscosity (kg m−1 s−1)

ρ =

fluid density (kg m−3)

ω =

nondimensional frequency

Subscripts
f =

fluid

i =

interface

o =

open region

s =

solid

w =

wall

References

1.
Hetsroni
,
G.
,
Gurevich
,
M.
, and
Rozenblit
,
R.
,
2006
, “
Sintered Porous Medium Heat Sink for Cooling of High-Power Mini-Devices
,”
Int. J. Heat Fluid Flow
,
27
(
2
), pp.
259
266
.10.1016/j.ijheatfluidflow.2005.08.005
2.
Panse
,
S. S.
,
Singh
,
P.
, and
Ekkad
,
S.
,
2019
, “
Air-Based Cooling in High Porosity, Aluminum Foams for Compact Electronics Cooling
,” 18th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (
ITherm
), Las Vegas, NV, May 28–31, pp.
376
383
.10.1109/ITHERM.2019.8757458
3.
Aliyu
,
M. D.
, and
Chen
,
H. P.
,
2018
, “
Enhanced Geothermal System Modelling With Multiple Pore Media: Thermo-Hydraulic Coupled Processes
,”
Energy
,
165
, pp.
931
948
.10.1016/j.energy.2018.09.129
4.
Torabi
,
M.
,
Karimi
,
N.
,
Torabi
,
M.
,
Peterson
,
G. P.
, and
Simonson
,
C. J.
,
2020
, “
Generation of Entropy in Micro Thermofluidic and Thermochemical Energy Systems-A Critical Review
,”
Int. J. Heat Mass Transfer
,
163
, p.
120471
.10.1016/j.ijheatmasstransfer.2020.120471
5.
Tiab
,
D.
, and
Donaldson
,
E. C.
,
2016
, “
Porosity and Permeability
,”
Petrophysics
,
4
, pp.
67
186
.
6.
Hristov, J.
,
2019
, “
Bio-Heat Models Revisited: Concepts, Derivations, Nondimensalization and Fractionalization Approaches
,”
Front. Phys.
,
7
, p.
189
.10.3389/fphy.2019.00189
7.
Pal
,
S.
,
Hajj
,
M.
,
Wong
,
W.
, and
Puri
,
I.
,
2014
, “
Thermal Energy Storage in Porous Materials With Adsorption and Desorption of Moisture
,”
Int. J. Heat Mass Transfer
,
69
, pp.
285
292
.10.1016/j.ijheatmasstransfer.2013.10.023
8.
Gonçalvès
,
J.
,
Adler
,
P. M.
,
Cosenza
,
P.
,
Pazdniakou
,
A.
, and
de Marsily
,
G.
,
2015
, “
Semipermeable Membrane Properties and Chemomechanical Coupling in Clay Barriers
,”
Dev. Clay Sci.
,
6
, pp.
269
327
.10.1016/B978-0-08-100027-4.00008-5
9.
Hunt
,
M.
, and
Tien
,
C.-L.
,
1988
, “
Non-Darcian Convection in Cylindrical Packed Beds
,”
ASME J. Heat Transfer-Trans. ASME
,
110
(
2
), pp.
378
384
.10.1115/1.3250495
10.
Parhizi
,
M.
,
Torabi
,
M.
, and
Jain
,
A.
,
2021
, “
Local Thermal Non-Equilibrium (LTNE) Model for Developed Flow in Porous Media With Spatially-Varying Biot Number
,”
Int. J. Heat Mass Transfer
,
164
, p.
120538
.10.1016/j.ijheatmasstransfer.2020.120538
11.
Yang
,
K.
, and
Vafai
,
K.
,
2010
, “
Analysis of Temperature Gradient Bifurcation in Porous Media – An Exact Solution
,”
Int. J. Heat Mass Transfer
,
53
(
19–20
), pp.
4316
4325
.10.1016/j.ijheatmasstransfer.2010.05.060
12.
Yang
,
K.
, and
Vafai
,
K.
,
2011
, “
Restrictions on the Validity of the Thermal Conditions at the Porous-Fluid Interface—An Exact Solution
,”
ASME J. Heat Transfer-Trans. ASME
,
133
(
11
), p.
112601
.10.1115/1.4004350
13.
Amiri
,
A.
, and
Vafai
,
K.
,
1994
, “
Analysis of Dispersion Effects and Non-Thermal Equilibrium non-Darcian, Variable Porosity Incompressible Flow Through Porous Medium
,”
Int. J. Heat Mass Transfer
,
37
(
6
), pp.
939
954
.10.1016/0017-9310(94)90219-4
14.
Alazmi
,
B.
, and
Vafai
,
K.
,
2002
, “
Constant Wall Heat Flux Boundary Conditions in Porous Media Under Local Thermal Non-Equilibrium Conditions
,”
Int. J. Heat Mass Transfer
,
45
(
15
), pp.
3071
3087
.10.1016/S0017-9310(02)00044-3
15.
Alomar
,
O. R.
,
2019
, “
Analysis of Variable Porosity, Thermal Dispersion, and Local Thermal Non-Equilibrium on Two-Phase Flow Inside Porous Media
,”
Appl. Therm. Eng.
,
154
, pp.
263
283
.10.1016/j.applthermaleng.2019.03.069
16.
Ahmed
,
H. E.
,
Fadhil
,
O. T.
, and
Salih
,
W. A.
,
2019
, “
Heat Transfer and Fluid Flow Characteristics of Tubular Channel Partially Filled With Grooved Metal Foams
,”
Int. Comm. Heat Mass Transfer
,
108
, p.
104336
.10.1016/j.icheatmasstransfer.2019.104336
17.
Xu
,
H. J.
,
2020
, “
Thermal Transport in Microchannels Partially Filled With Micro-Porous Media Involving Flow Inertia, Flow/Thermal Slips, Thermal Non-Equilibrium and Thermal Asymmetry
,”
Int. Comm. Heat Mass Transfer
,
110
, p.
104404
.10.1016/j.icheatmasstransfer.2019.104404
18.
Dehghan
,
M.
,
Valipour
,
M. S.
, and
Saedodin
,
S.
,
2016
, “
Conjugate Heat Transfer Inside Microchannels Filled With Porous Media: An Exact Solution
,”
J. Thermophys. Heat Transfer
,
30
(
4
), pp.
814
824
.10.2514/1.T4767
19.
Li
,
Y.
,
Gong
,
L.
,
Ding
,
B.
,
Xu
,
M.
, and
Joshi
,
Y.
,
2021
, “
Thermal Management of Power Electronics With Liquid Cooled Metal Foam Heat Sink
,”
Int. J. Therm. Sci.
,
163
, p.
106796
.10.1016/j.ijthermalsci.2020.106796
20.
Pavel
,
B. I.
, and
Mohamad
,
A. A.
,
2004
, “
An Experimental and Numerical Study on Heat Transfer Enhancement for Gas Heat Exchangers Fitted With Porous Media
,”
Int. J. Heat Mass Transfer
,
47
(
23
), pp.
4939
4952
.10.1016/j.ijheatmasstransfer.2004.06.014
21.
Al-Aabidy
,
Q.
,
Alhusseny
,
A.
, and
Al-Zurfi
,
N.
,
2021
, “
Numerical Investigation of Turbulent Flow in a Wavy Channel Partially Filled With a Porous Layer
,”
Int. J. Heat Mass Transfer
,
174
, p.
121327
.10.1016/j.ijheatmasstransfer.2021.121327
22.
Ochoa-Tapia
,
J. A.
, and
Whitaker
,
S.
,
1997
, “
Heat Transfer at the Boundary Between a Porous Medium and a Homogeneous Fluid
,”
Int. J. Heat Mass Transfer
,
40
(
11
), pp.
2691
2707
.10.1016/S0017-9310(96)00250-5
23.
Amiri
,
A.
,
Vafai
,
K.
, and
Kuzay
,
T. M.
,
1995
, “
Effect of Boundary Conditions on non-Darcian Heat Transfer Through Porous Media and Experimental Comparisons
,”
Numer. Heat Transfer J. A
,
27
(
6
), pp.
651
664
.10.1080/10407789508913724
24.
Matuła
,
I.
,
Dercz
,
G.
, and
Barczyk
,
J.
,
2020
, “
Titanium/Zirconium Functionally Graded Materials With Porosity Gradients for Potential Biomedical Applications
,”
Mater. Sci. Technol.
,
36
(
9
), pp.
972
977
.10.1080/02670836.2019.1593603
25.
Ma
,
P.
,
Wang
,
B.
,
Chen
,
S.
,
Zhang
,
X.
,
Tao
,
C.
, and
Xing
,
X.
,
2018
, “
Numerical Investigation of Heat Transfer Enhancement Inside the Pipes Filled With Radial Pore-Size Gradient Porous Materials
,”
ASME J. Therm. Sci. Eng. Appl.
,
10
(
5
), p.
054502
.10.1115/1.4040276
26.
Kaviany
,
M.
,
1995
,
Principles of Heat Transfer in Porous Media
, Springer, Berlin.https://link.springer.com/book/10.1007/978-1-4612-4254-3
27.
Beavers
,
G. S.
, and
Joseph
,
D. D.
,
1967
, “
Boundary Conditions at a Naturally Permeable Wall
,”
J. Fluid Mech.
,
30
(
1
), pp.
197
207
.10.1017/S0022112067001375
28.
Sahraoui
,
M.
, and
Kaviany
,
M.
,
1992
, “
Slip and No-Slip Velocity Boundary Conditions at Interface of Porous, Plain -Media
,”
Int. J. Heat Mass Transfer
,
35
(
4
), pp.
927
943
.10.1016/0017-9310(92)90258-T