A Pore-Scale Investigation of the Transient Response of Forced Convection in Porous Media to Inlet Ramp Inputs

Comprehending the transient response of porous media has grown in importance with its widespread use in emerging technologies [1–5]. Cases can be found in the use of heat exchangers [6], micro chemical reactors [7–9], fuel cells [10,11], and electrochemical systems [12,13]. Understanding the response of these systems is vital as the inlet fluid flowrate can be highly time-dependent [14]. It is fundamental to have the ability to predict the thermal response of the pertinent model subject to time-dependent disturbances at the inlet [15]. Steady flow has been the predominant focus of majority of existing studies of forced convection in porous media. This pattern is clearly visible in numerical investigations that were conducted using a pore-scale approach (e.g., Refs. [16–18]) and especially in those studies that took a macroscopic modeling approach (e.g., Refs. [19–22]).

A smaller amount of macroscopic investigations have analyzed unsteady and pulsating flows in porous media as will be discussed in the core part of this literature review. Fewer studies reported results using the pore-scale approach, with most studies being in the two-dimensional domain [8,13,23,24] and even less in the three-dimensional domain [9,25]. In order to fully comprehend steady and unsteady macroscopic modeling, a brief examination of literature on macroscopic studies in porous media is carried out using the Darcy model and its extensions. It is generally understood that a macroscopic model might be lacking for the analysis of periodic flow systems [12]; therefore, a pore-scale approach is employed in this paper.

Saito and de Lemos [26] proposed two-macroscopic-energy equation models for conduction and convection for packed beds of porous media assuming local thermal non-equilibrium. The authors modeled a two-dimensional unit cell model with square cylinders in a staggered arrangement, replicating Kuwahara et al.’s [27] physical model. Saito and de Lemos [26] numerically solved for the interfacial heat transfer coefficient, discretizing the macroscopic governing equations using the control volume method. The authors of Ref. [26] found their work to be in agreement with Kuwahara et al. [27]. They concluded that a higher porosity medium could reduce the energy transfer between the solid and fluid phase [26]. Lu et al. [28] simulated temperature disturbances during the mixing process of a T-junction filled with and without porous media. They implemented local thermal equilibrium model and concluded that the addition of porous media within the T-junction model lowers the temperature fluctuations during the mixing process [28].

Kim et al. [29] numerically simulated a two-dimensional channel filled with porous media to investigate the effects of forced oscillatory flow on the Nusselt number. Isothermal boundary conditions were set to the channel walls where the porosity of the porous media was fixed to 0.6. The authors of Ref. [29] compared the temperature and flow fields with ranges of pertinent parameters and discovered oscillatory flow to reduce heat transfer in the entrance region of the model when compared to steady flow. They further reported enhanced heat transfer in the central region of the model. A similar study was conducted by Guo et al. [30] who considered the effects of pulsating flow on heat transfer in a pipe partially filled with porous media. The system was modeled as two-dimensional with porous media attached to the walls of the pipe. Guo et al. [30] found that the Nusselt number monotonically increases with the growth in the porous layer thickness.

Bhargava et al. [31] present finite element solutions for a non-Newtonian fluid experiencing pulsation. The analysis was conducted by adopting a one-dimensional transient model depicting a porous medium conduit. The pulsatile effects of the pressure gradient on the model flow were examined by altering parameters such as Reynolds number, pulsating amplitude, and frequency. Bhargava et al. [31] discovered the velocity to augment with an increase in permeability and a reduction of non-Newtonian behavior. Huang and Yang [32] studied forced oscillatory convection in a parallel plate channel with porous blocks mounted to the lower wall. They described the effects of varying the amplitude, frequency, and porous blockage ratio on the flow field and heat transfer of the model. An oscillatory flow was introduced to improve flow mixing and greater thermal transport, revealing that heat transfer augments from the heat sources with the rise of the oscillation amplitude.

A numerical study was conducted by Dhahri et al. [33] to find the effects of pulsating flow on entropy generation within a pipe partially filled with porous media. The authors monitored the response of the amplitude, frequency, porous layer thickness, and thermal conductivity ratio with respect to entropy generation. They discovered the temperature difference in the fluid region to be greater than that of the porous region. Farooghi et al. [34] studied steady and oscillating flow in a parallel plate channel partially composed with porous media in the upper and lower walls of the channel, subjecting to a constant heat flux. The authors employed the Brinkman-Forchheimer-extended Darcy model and implemented the local thermal non-equilibrium (LTNE) method and reported that local thermal equilibrium (LTE) method overestimates the heat transfer. Farooghi et al. [34] suggested that augmenting the amplitude of flow oscillation could result in an increase in the time-averaged Nusselt number. They further noted the thermal conductivity ratio to increase heat transfer between the solid and the fluid, if the porous media thickness is kept constant.

Ghafarian et al. [35] computed pulsating flow through a channel filled with a porous medium to analyze the heat transfer. Various studies were conducted by altering the frequency, amplitude, and Reynolds number of the working fluid as well as a change of porosity. The local Nusselt number was calculated for a laminar model to monitor heat transfer, showing that the Nusselt number to augment with the use of a high amplitude and high frequency. Targui and Kahalerras [6] numerically simulated oscillatory flow through a double pipe heat exchanger coupled with porous baffles. It was found that the oscillation has a profound effect on heat transfer and maximizes it when only the hot fluid flow oscillates. Zaman et al. [36] numerically examined the flow of blood induced with pulsatile flow within a porous-saturated overlapping stenosed artery. The study depicted a pathological solution where the blood travels through an artery comprising of blood clots and fatty cholesterol. The numerical simulation showed the effects of flowrate, velocity, and wall shear stress on the porous medium. It was found that with a larger permeability the resistance of the flow is smaller whilst the flowrate and velocity augment.

It is observed that the existing studies on unsteady forced convection are mostly concerned with the macroscopic analysis by using the extended Darcy models. Further, most studies include two-dimensional simulations and gave little attention to the transient response of the system. Therefore, there is currently a gap in understanding the thermal transient response at pore scale in a three-dimensional domain. The current investigation aims to address this issue.

A representative elementary volume (REV) similar to the structure used in Ref. [17] is taken into consideration. The REV contains a single pore around an obstacle [38]. Such a unit is periodically repeated in each direction to create the large-scale porous structure. Therefore, it is safe to assume nearly identical flows in each REV unit. As a result, the boundaries of this unit element are treated as periodic and symmetric.

A series of steady-state investigations were performed on the terrain model to study the effects of main variables affecting the fluid flow inside the system. Three porosities of 0.87, 0.8, and 0.72 were considered. Uniform velocity and 300 K constant temperature were set at the inlet and atmospheric pressure was set at the outlet. Three wall temperatures of 500, 600, and 700 K were also tested with 300 K inlet temperature. The effects of circular obstacle geometries were investigated, and the Nusselt number was calculated separately on each obstacle to enable a more detailed evaluation of the results.

Along with the aforementioned Re and porosity values, three ramp durations of 10, 20, and 30 s were investigated along with three ramp-up levels of 10, 20, and 30% of the base value. All the investigations were tested for ramp-ups in inlet temperature and velocity. Furthermore, three fluids as air, CO₂, and H₂ were tested. A total of around 500 different settings were investigated in the course of this investigation.

A numerical model was developed using the commercial software star-ccm+. A satisfactory refined mesh of polyhedral cells was formed for the bulk region (Fig. 3(a)). In total, 15 high-quality, tetrahedral, or “prism layer” cells were generated per solid surface to ensure higher accuracy of gradient calculations. All equations were solved using a second-order discretization technique where a coupled solver was implemented for greater stability. The boundary conditions were set according to each case as described in Sec. 2.1. Stagnant flow along with ambient pressure and temperature were set as the initial conditions for all cases except for the unsteady simulations. For these cases, the converged solution from similar steady-state simulation was used instead. The central time-step chosen was that in which the magnitude was lesser by a factor of 100 in comparison with that of the physical timescale. In order to reduce the runtime of ramp cases, a variable time-step approach was adopted. These time-steps were refined around the ramp section for more accuracy.

A minimum magnitude of 10⁻⁶ for all equations residual levels were achieved for all steady-state simulations. Furthermore, a sizeable parameter was selected for upper limit time-step number in ramp cases, and mean temperature at the exit was analyzed to monitor convergence with respect to time. A new stoppage criterion was defined to terminate the simulations if the deviation in the last 1000 time-steps is lower than 0.5 K [42]. These proved to be adequate convergence criteria for almost all cases.

RLR is a non-dimensional quantity, calculated for the Nusselt number response of each obstacle. A graphical representation for the RLR definition is shown in Fig. 3(b).

Equation (10) was derived to represent a staggered arrangement whereas Eq. (11) was extracted particularly for square obstacles. Both correlations can be used for a large range of properties and Reynolds numbers. The close similarity of the simulated results and the correlations confirms the validity of the numerical simulations.

The computational mesh was systematically refined to find the ideal balance between the relative accuracy and low computational costs. A total of seven grids were created between a range of 2600–2,570,000 cells. Early investigations on a single Reynolds number showed that the Nusselt number shows minimal variance for grid ranges between 225,000–2,570,000 cells. However, further refinement in the grid sizes did not show any significant change in the results (as shown in Fig. 4). Therefore, it is safe to argue that the results are grid independent using the 613,000 cell grid.

Figure 5 illustrates the spatial response of the ten obstacles pore-scale model to a steady inlet flow velocity. This figure displays snapshots of three configurations where the wall temperature alters at intervals of a 100 K for their respective converged solutions. The inlet flow enters at ambient temperature from the left-hand side and exits toward the right. As the inlet flow travels downstream, the constant temperature around the circular obstacles influences the flow field around and behind the obstacles, creating wake regions. This effect is readily observed in the flow field where there is a larger temperature difference between the inlet flow and the circular obstacles and is particularly visible when the wall temperature is at its lowest. Such effects decay after the third obstacle as the flow travels downstream almost appearing to have reached thermal equilibrium. However, the diffusion effects are the strongest throughout the porous configuration when the obstacle temperature is the largest. This can be seen when the wall temperature is 700 K, and the flow field appears to reach very close to thermal equilibrium further downstream at the tenth obstacle. The flow field behavior and temperature change is monitored throughout the system by the surface-averaged Nusselt Number around each circular obstacle. The analyses discussed in this section aim to understand the transient response by observing the response lag ratio (RLR) of the pore-scale model subject to an unsteady ramp input.

A comparison of Nusselt numbers is made for the inlet flow subject to a ramp between inlet temperature and inlet velocity in Fig. 6(a). As the velocity ramp is introduced, the first obstacle responds gradually albeit the largest increase in the Nusselt number is also reported. The opposite is observed at obstacles closer to the outlet with quicker response time and lesser augmentation in the Nusselt number. This is expected to be due to the thermal boundary conditions set for the circular obstacles. As the time-dependent ramp travels throughout the model, the viscous effects continue to decay with the rate of heat transfer plummeting as it reaches the outlet. A similar trend is observed when the system is subject to an inlet temperature ramp. To achieve constant convective heat transfer, the system requires less time when the flow is closer to the outlet when compared to the inlet. As the inlet temperature rises, the local temperature difference between the fluid and solid surfaces decreases. This limits the rate of heat transfer and reduces the Nusselt number when compared with the steady-state value. It is noteworthy that the heat transfer rate per obstacle and response time both change when the system is subject to a velocity ramp whereas only response time differs when the system is subject to a temperature ramp. Figure 6(b) displays the response lag ratio for each obstacle. The lag augments as the fluid travels downstream of the model where the viscous effects resist the inertial forces of the fluid flow. As the Reynolds number increases, the viscous effects are also amplified, inflating the overall lag of the system. This phenomenon is further analyzed in Fig. 7, visualizing the spatiotemporal flow field snapshots equally split over the respective ramp duration. The time-dependent inlet temperature is subject to an increase of 30% of its base value. In Fig. 7(a), it can be seen before the introduction of the ramp, and the flow field of the obstacles closer to the exit of the system exhibits higher temperature readings, almost matching the temperature of the surface of the obstacles (700 K). As the inlet ramp in temperature is introduced, the overall flow field temperature surges, forcing the higher temperatures to move further upstream. This is because the convective heat transfer diminishes as the gap between the temperature difference narrows, decreasing the Nusselt number.

However, the opposite is observed in Fig. 8 where the temporal evolutions for the velocity field are introduced by a 10-s ramp input. Thus, the increase in velocity and constant temperature at the circular obstacles provide a streamlined contour pattern where convective heat transfer is increasing. Also, as the Reynolds number increases further with the ramp so does the inertial effects of the flow, combined with the thermal diffusion occurring at the circular obstacles, the wake region expands gradually with the passage of time. It is apparent from Fig. 8 that there is a greater wake flow when compared to Fig. 9, which is lower in the Reynolds number. This can also be observed between Figs. 9(a) and 9(e) where the wake flow contours surge in velocity as the Reynolds number augments. A similar pattern is also observed in the general flow between the two configurations, where a more streamline contour pattern is observed as the flow moves toward the outlet for a high Reynolds number and a more scattered wider flow pattern for the lower Reynolds number.

An identical configuration to Fig. 8 is presented in Fig. 10 where the working fluid has been switched from air to hydrogen. A similar streamlined contour flow field is observed to that of air; however, there is a drastic increase in the system's velocity albeit the same Reynolds number. This is primarily due to the fluid properties. Hydrogen inhibits a lower density than air, making the fluid lighter therefore more responsive to the identical Reynolds number. This trend is clearly observed during the general flow of the model and within the wake region as the ramp input is implemented. It is understood for a high-density working fluid with identical configuration in the pore-scale model, and the system will exhibit a lower fluid velocity, therefore making it less sensitive to inlet disturbances in comparison with a low-density fluid.

Although the fluid properties of hydrogen make it more sensitive to inlet disturbances, it also appears to be a lot more responsive to the thermal diffusion and viscous effects incurred around the obstacles when compared with air. This is observed in the response lag ratio of hydrogen in Fig. 11(a), where all respective porosities report a higher response lag ratio when compared with its counterpart, air. A monotonic increase in the response lag ratio is observed for air across all three porosities of the system. However, hydrogen exhibits a response lag ratio further downstream for larger porosities. The change in porosity amongst different working fluids is further analyzed in Fig. 11(b). The common trend changed where the largest porosity outputs a higher response lag ratio for high-density fluids (carbon dioxide and air), and the opposite is noted for hydrogen when compared with the lowest porosity. As the Reynolds number is increased with decreasing porosity, the RLR decreases rather than increasing for air and carbon dioxide due to their fluid properties. With hydrogen being more sensitive and less dense, the change in the Reynolds number with porosity has little or no effect in the overall RLR.

The RLR is further discussed under the variable parameters of amplitude and ramp duration for the different working fluids. Figure 12(a) shows that the increase in amplitude augments the inlet velocity and therefore increasing forced convective heat transfer and viscous effects around the obstacles, resulting in a higher response lag ratio. Fluid properties continue to impact the degree of RLR reported, with hydrogen (lower density fluid) showcasing larger increases in RLR as the amplitude increases. However, for air (higher density fluid), although there is an increase in RLR, it is minimal. A general trend is also observed in Fig. 12(b), as the ramp duration is prolonged amongst the working fluid the RLR is noted to decrease. The lower RLR is reported as the gradual increase in velocity due to longer ramp duration decays the viscous effects and allows more time for heat to transfer from the circular obstacles to fluid with ease. Contrary to the previous trend reported in Fig. 12(a), Fig. 12(b) showcases higher RLR for air than hydrogen in general. Yet, it is also more responsive to the increase in ramp duration allowing the system to adapt to a higher density fluid resulting in a lower RLR than hydrogen for a 30-s ramp duration. The RLR being the response lag ratio is representative of the attenuation effect of the porous medium upon an input disturbance. Such attenuation is generally more intense when the disturbance includes strong gradients. For a ramp disturbance, strong gradients imply higher amplitude and shorter duration. That is why in Fig. 12 increase of Ramp amplitude results in higher values of RLR while increases in ramp duration reduces that.

A pore-scale analysis was conducted to study the effects of forced convection in porous media subject to a ramp disturbance superimposed on the entrance of a reticulated porous structure. Comprehending the transient response of porous media has grown in importance with its widespread use in numerous emerging applications. This study allows the prediction of thermal response of the pertinent system subject to time-dependent disturbances at the inlet. The porous model was made up of ten consecutive cylindrical obstacles with set isothermal boundary conditions and was subject to ramp disturbance for the flow velocity and temperature at the inlet. The RLR and transient behavior at each obstacle was examined in detail by monitoring the spatiotemporal evolution of the flow and temperature fields as well as the surface-averaged Nusselt number.

They key findings are as follows:

• An imposed ramp temperature disturbance at the inlet has an identical thermal response throughout the system, where an inlet velocity ramp input shows a varying response for each cylindrical obstacle. This makes a pore-scale analysis essential to monitor the transient thermal response of a porous structure.

• An increase in the Reynolds number has a minuscule rise on the RLR.

• An increase in porosity decreases the RLR. Yet, when the Reynolds number is decreased for a structure with a large porosity, RLR increases for high-density fluids only (air and carbon dioxide).

• An increase in the amplitude increases the RLR.

• An increase in the ramp duration decreases the RLR, more strongly for high-density fluids.

N. Karimi and B. Yadollahi acknowledge the financial support of Engineering and Physical Science Research Council (Grant No. EP/N020472/1).

a =

amplitude = max · u − min · u

k =

thermal conductivity (W K⁻¹ m⁻¹)

p =

pressure (Pa)

t =

time (s)

u =

flow velocity (m s⁻¹)

D =

obstacle diameter (m)

H =

height (m)

T =

temperature (K)

c_p =

specific heat capacity (J K⁻¹ kg⁻¹)

h_o =

external heat loss coefficient (W m⁻² K⁻¹)

q″=

heat flux (W m⁻²)

Nu=

Nusselt number (−)

Pr=

Prandtl number (−)

Re=

Reynolds number (−)

ɛ =

porosity

μ =

dynamic viscosity (kg m⁻¹ s⁻¹)

ρ =

density (kg m⁻³)

ψ =

temperature/velocity

f =

fluid

in=

inlet

L =

based on obstacle diameter size

ramp=

function of ramp

ref=

reference

x, y, z =

Cartesian coordinates