Local Thermal Nonequilibrium (LTNE) Modeling of a Partially Porous Channel With Spatial Variation in Biot Number

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 [10–13]. 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 [11–14,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.

Here, $Tf$ and $Ts$ are fluid and solid phase temperatures in the porous medium, respectively, whereas u_p, $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.

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

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

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), u_i is the interface pure fluid velocity, u_p,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].

where q_i 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.

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

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

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

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.

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

The final model considered here assumes an interface heat transfer coefficient h_int [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 h_int.

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.

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, η₁ = 0.98, Bi = 10, L_s = 0, L_f = 0, α* = 0.1, and Da = 10⁻⁵ for model A. In addition, Bi_i = 0.1 for model C. Note that η₁ 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 η = η₁ 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 = Bi₀(1−η²), Da = 10⁻⁵, Bi₀ = 30, α* = 0.1, L_s = 0, and L_f = 0, Figs. 3(a) and 3(b) plot fluid and solid temperature distributions for this work with different values of η₁. 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 η₁ 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, 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.

Figure 4 plots fluid and solid temperature distributions for a quadratically varying Biot number, given by Bi = Bi₀ (1 + η²). The effect of the magnitude of the Biot number function, Bi₀ is investigated first. Other problem parameters for this figure are k = 1, η₁ = 0.98, β = 0.9, Bi_i = 1, L_s = 0, L_f = 0, α* = 0.1, and Da = 10⁻⁵. 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 η = η₁. 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 Bi₀ = 0.5 and Bi₀ = 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 Bi₀ = 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.

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 Bi₀ 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 Bi₀, whereas the fluid temperature distribution goes up with increasing Bi₀. 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.

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, η₁ = 0.98, β = 0.9, Bi_i = 1, L_s = 5, L_f = 0, α* = 0.1, and Da = 10⁻⁵. 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.

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.

Temperature distributions predicted by the three models for different thicknesses of the porous medium, i.e., η₁, are plotted in Figs. 8 and 9, for quadratically varying (Bi = Bi₀ (1 + η²)) and periodically varying (Bi = 50(1 + cos(2πωη)) Biot number distributions, respectively. The problem parameters used for these figures are k = 10, Da = 10⁻⁵, Bi_i = 1, Bi₀ = 0.5, L_s = 0, L_f = 0, α* = 0.1, and ω = 1. As expected, for each model, the temperature difference between the solid and fluid phases reduces as η₁ 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 η₁ = 0.96. This is due to a reduction in the imposed heat flux at the interface as η₁ is reduced. As seen in Fig. 8, for the quadratically varying Biot number case, reducing η₁ 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 η₁, 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.

The interfacial Biot number, Bi_i is a key nondimensional parameter governing model C. Bi_i defines the rate of heat transfer from the pure fluid region to the solid phase at the interface, η = η₁, as shown in Eq. (52). Figure 10 plots the temperature distributions for model C to analyze the effect of Bi_i. 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 Bi_i, 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.

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 L_s, 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, η₁ = 0.98, β = 0.9, Bi_i = 1, ω = 2, L_f  =  0, α* = 0.1, and Da = 10⁻⁵. Three different values for L_s 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 L_s has a more significant impact on the solid temperature. As expected, the temperature difference is highest at points where the Biot number is minimum.

In Fig. 12, the effect of velocity slip coefficient on the temperature distribution is analyzed using model B for a quadratically varying (Bi = Bi₀ (1 + η²)) Biot number distribution. Problem parameters used to analyze the effect of the velocity slip coefficient are k = 1, η₁ = 0.98, β = 0.9, Bi₀ = 0.5, L_s = 0, L_f = 0, and Da = 10⁻⁵. 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.

Finally, temperature distributions for model B are plotted in Fig. 13 for different values of the Darcy number. Quadratically varying Biot number distribution (Bi = Bi₀ (1 + η²)) is used for this figure. Other parameters used to study the impact of Darcy number are k = 1, η₁ = 0.98, β = 0.9, Bi₀= 0.5, L_s = 0, L_f = 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⁻¹⁰.

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.

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

Bi =

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

Bi_i =

interfacial Biot number, $Bii=hiHks,eff$

C_p =

specific heat of the fluid (J kg⁻¹ K⁻¹)

Da =

Darcy number, $Da=KH2$

H =

half height of the channel (m)

h_i =

interfacial heat transfer coefficient at the interface between porous and nonporous regions (W m⁻² K⁻¹)

h_int =

interstitial solid–liquid heat transfer coefficient in the porous region (W m⁻² K⁻¹)

H₁ =

half height of the porous medium (m)

k =

thermal conductivity ratio, $k=kf,effks,eff$

K =

permeability (m⁻²)

k_eff =

effective thermal conductivity (W m⁻¹ K⁻¹)

k_f =

thermal conductivity of the fluid (W m⁻¹ K⁻¹)

k₁ =

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

L =

$nondimensional internal heat generation, L=SH1qi$

p =

pressure (N m⁻²)

q_i =

heat flux at the interface (W m⁻²)

q_w =

heat flux at the wall (W m⁻²)

S =

internal heat generation (W m⁻³)

T =

temperature (K)

u =

fluid velocity (m s⁻¹)

U =

dimensionless fluid velocity, $U=u−dPdx · H2μf$

u_i =

interface pure-fluid velocity (m s⁻¹)

U_i =

dimensionless interface pure fluid velocity

U_m =

dimensionless average fluid velocity over the entire channel

x =

horizontal coordinate (m)

y =

vertical coordinate (m)

α =

specific interfacial area of the porous medium (m⁻¹)

α* =

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(T−Ts,i)qwH$

μ =

dynamic viscosity (kg m⁻¹ s⁻¹)

ρ =

fluid density (kg m⁻³)

ω =

nondimensional frequency

Subscripts

f =

fluid

i =

interface

o =

open region

s =

solid

w =

wall