Flow and Heat Transfer Characteristics in Porous Media with Explicit Structure: A Multi-Physical Field Coupling Study

1  Introduction

Porous media with complex internal structures are widely encountered in a broad range of engineering applications, including thermal management systems, energy conversion devices, and chemical processing units such as catalytic reactors, heat exchangers, and filtration systems. The coupled transport of momentum and heat within such media plays a crucial role in determining system efficiency and operational stability. In particular, under conditions involving forced convection and strong thermal gradients, the interaction between flow resistance and heat transfer becomes highly nonlinear and sensitive to internal structural characteristics [1–3].

Conventional modeling approaches for porous media typically rely on volume-averaged formulations, in which the detailed internal structure is replaced by effective macroscopic parameters such as permeability, porosity, and thermal conductivity. Classical models based on Darcy’s law and its extensions, including the Brinkman and Forchheimer corrections, have been extensively applied due to their simplicity and computational efficiency [1,4]. However, these homogenized approaches inherently neglect pore-scale geometric details and are therefore unable to accurately capture local flow heterogeneity, preferential pathways, and temperature non-uniformity [5,6]. More recently, hybrid pore-network-continuum (PNM-continuum) frameworks have been developed to bridge pore-scale resolution and continuum-scale efficiency by coupling explicit pore network models with Darcy-scale representations of microporous regions [7,8]. While these methods significantly reduce computational cost and enable multiscale transport simulations, their applicability to densely packed, low-porosity structured media has not been fully established. This limitation becomes particularly pronounced in structured or layered porous systems, where anisotropic geometries and local constrictions strongly influence macroscopic transport behavior.

To overcome the limitations of homogenized models, increasing research efforts have been devoted to resolving the internal structure of porous media through high-fidelity numerical simulations. Pore-scale computational fluid dynamics (CFD) methods, including direct numerical simulation (DNS) and lattice Boltzmann methods, have demonstrated strong capability in capturing detailed transport phenomena such as local velocity fluctuations, vortex dynamics, and thermal gradients in both ordered and randomly heterogeneous porous structures [9–11]. Recent pore-scale studies have further revealed that vortex formation in pore-throat expansion zones and the resulting convective enhancement are strongly dependent on porosity, with a sharp increase in the Peclet number observed when porosity exceeds 0.7 [11]. Furthermore, permeability has been identified as a critical parameter governing forced convective heat transfer in rock media with a porosity of approximately 0.24, where higher Nusselt numbers are consistently observed at lower Darcy numbers [12]. However, fully hybrid explicit-structure approaches remain computationally expensive and impractical for engineering-scale simulations involving large arrays of discrete porous elements. As an alternative, hybrid modeling strategies that combine explicit structural representation with porous media theory have gained increasing attention in recent years [7,8,13]. These approaches allow for partial resolution of geometric features while retaining reasonable computational cost, thereby providing a balance between physical fidelity and efficiency.

Another layer of complexity arises from the presence of multicomponent substances within porous structures, such as water, glycerin, and other volatile compounds. The transport behavior of these components is highly sensitive to local temperature gradients and flow field variations, particularly under heating conditions where evaporation or phase-change processes may occur [14–16]. Even in the absence of explicit phase-change modeling, the strong coupling between thermal and flow fields can lead to significant spatial variations in temperature distribution. Homogenized models often fail to capture these localized effects, highlighting the necessity of incorporating structural information into the modeling framework [17,18].

In parallel with these numerical developments, experimental investigations have provided valuable insights into flow resistance and heat transfer characteristics across diverse porous structures. Pressure drop measurements over anisotropic 3D-printed lattice substrates have demonstrated that streamwise-preferential permeability can effectively modulate drag response over a wide Reynolds number range (Re = 500–4000), with the wall-normal permeability identified as a critical parameter governing the onset of flow instabilities [19]. Studies on sintered metal porous plates have combined experimental flow resistance data with numerical simulations to validate Darcy-Forchheimer model predictions [20]. Meanwhile, methodologies for systematically evaluating Forchheimer resistance coefficients for permeable screens and wire mesh panels have been proposed, enabling more accurate macroscopic representations of complex thin porous structures [21]. In the context of wire woven mesh porous media, recent REV-scale simulations incorporating local thermal non-equilibrium (LTNE) have revealed the trade-off between enhanced heat transfer and increased pressure drop across different stacking configurations and porosity ranges [13]. The emergence of topology optimization and microstructure-informed modeling techniques has further enabled more accurate representation of complex porous geometries and their associated transport behavior [7,22]. A quantitative comparison of representative recent studies is provided in Table 1.

Despite these advancements, several critical research gaps remain. First, as summarized in Table 1, the majority of existing pore-scale and hybrid-scale studies have focused on media with relatively high porosity, typically exceeding 0.24, while structured systems at porosities of 0.1 to 0.3 composed of densely packed discrete elements remain largely unexplored. In such systems, local geometric features, including narrow flow passages and abrupt cross-sectional variations, may induce flow channeling, jet-like acceleration, and non-uniform heat transfer that cannot be adequately captured by either fully homogenized continuum models or conventional pore-network-continuum frameworks [6,8]. Second, few studies have systematically compared the relative importance of porosity, flow rate, and heating power on both flow resistance and thermal distribution within a single modeling framework. Third, the validity of the local thermal equilibrium (LTE) assumption has been examined primarily for high-porosity metal foams and packed beds [23], whereas its applicability to low-porosity, multi-layered sheet arrays with internal heat sources remains to be quantitatively assessed. Finally, the long-term structural integrity of explicit-structure porous media under repeated thermal cycling, which is inherent to the target application involving periodic heating and cooling, has not been addressed in existing literature, representing a significant knowledge gap that warrants future investigation.

To address these challenges, the present study develops a three-dimensional explicit structure–porous media hybrid model for multi-layered porous sheet arrays. By explicitly reconstructing the geometric arrangement of the porous sheets and incorporating Darcy–Forchheimer resistance terms, the proposed model captures both macroscopic flow pathways and microscopic resistance effects. The coupled flow and heat transfer processes are systematically investigated under varying porosity, inlet velocity, and heating power conditions. In addition, experimental measurements are conducted to validate the numerical model and evaluate its predictive accuracy.

The primary objective of this work is to elucidate the mechanisms governing flow resistance and temperature distribution in structured porous media with explicit geometric representation. Particular emphasis is placed on identifying the transition between viscous-dominated and inertia-influenced flow regimes, as well as clarifying the role of structural heterogeneity in shaping local and global thermal behavior. While the present study focuses on steady-state flow and heat transfer under fixed structural configurations, future work should extend the proposed framework to account for cyclic thermal loading and potential structural degradation effects, which are relevant to the long-term performance of porous media in practical thermal management applications. The findings of this study are expected to provide useful insights for the design and optimization of thermal systems involving complex porous structures, and to contribute to the development of more accurate and physically representative modeling approaches.

2  Control Equations and Boundary Conditions

To investigate the coupled flow and heat transfer characteristics within the explicit-structure porous segment, a hybrid modeling strategy is adopted. The open channels between the sheets are treated as free-flow regions governed by the Navier–Stokes equations, while the porous sheets are modeled using an equivalent porous media formulation to account for their internal resistance effects. The working fluid is assumed to be air and treated as an incompressible Newtonian fluid. This assumption is justified by the low Mach number condition (Ma ≪ 0.3) and the relatively small pressure variations within the computational domain, under which density changes can be reasonably neglected. Additionally, gravitational effects are ignored due to the small characteristic length scale and the dominance of viscous forces. The governing equations for the free-flow region are expressed as:

∇⋅u=0(1)

ρ(u⋅∇)u=−∇p+μ∇2u(2)

ρcp(u⋅∇T)=∇⋅(k∇T)+Q˙(3)

where u is the velocity vector, p is pressure, ρ is density, μ is dynamic viscosity, cp is specific heat capacity, k is thermal conductivity, and Q˙ is the volumetric heat source; notably, Q˙ is set to 0 W/m3 throughout all computational domains in this study, as the thermal energy is introduced exclusively via the boundary heat source (Qb) applied to the central element surface, which is detailed in the subsequent section.

Within the porous sheets, the flow is significantly influenced by the solid matrix. To capture both viscous and inertial resistance effects, the Darcy-Forchheimer model is employed [4,24]. The Brinkman term (μ/ε λ2 u) is retained in the momentum equation to account for viscous shear effects at the porous-free flow interface, which is essential in the present hybrid formulation where explicit porous sheets are adjacent to open channel regions. At the low porosity and low Reynolds number conditions investigated, the Brinkman contribution is small relative to the Darcy resistance but ensures smooth stress continuity across the porous-free flow interface [25]. The governing equations for the porous region are written as:

∇⋅u=0(4)

ρε(u⋅∇)(uε)=−∇p+με∇2u−μαu−Cfρ|u|u(5)

ρcp,eff(u⋅∇T)=∇⋅(keff∇T)+Q˙(6)

where ε is porosity, α is permeability, and Cf is the inertial resistance coefficient. The term (μ/α)u represents the Darcy viscous resistance, while Cfρ|u|u accounts for the Forchheimer inertial loss. This formulation allows the model to capture the transition from viscous-dominated flow at low velocities to inertia-influenced behavior at higher velocities [25]. In the present study, the permeability and inertial resistance coefficient are set to α = 1 × 10−12 m2 and Cf = 0.55 for all three investigated porosity values, namely ε of 0.1, 0.2, and 0.3.

For the energy equation, the Local Thermal Equilibrium (LTE) assumption is adopted. The LTE assumption is considered reasonable under the present low-velocity and steady-state conditions, where the characteristic length scale is small. The Biot number (Bi) is evaluated as:

Bi=hsf⋅Lcks(7)

where hsf is the solid-fluid interfacial heat transfer coefficient, Lc is the characteristic length, and ks is the solid thermal conductivity. According to the authoritative heat transfer criteria outlined by Yang and Tao [26], the convective heat transfer coefficient for air under typical convection regimes spans from 20 to 100 W/(m2·K). Given the relatively low-velocity conditions investigated in the present study, an intermediate average value of hsf ≈ 62 W/(m2·K) is reasonably adopted herein. With Lc = 0.23 (corresponding to the sheet thickness) and ks = 0.21 W/(m·K), the calculated Biot number is Bi ≈ 0.068. Even evaluating at the upper bound of the standard range, the Biot number remains around 0.11, which is still on the order of 0.1 and thus firmly supports the validity of the core assumption. According to established heat transfer criteria for porous media [4,27], a Biot number satisfying Bi < 0.1 indicates that the internal thermal conduction resistance within the solid phase is negligible compared to the convective thermal resistance at the solid-fluid interface. Therefore, the LTE assumption serves as a justified engineering approximation under the investigated steady-state, low-velocity conditions.

This analysis provides an approximate validation of thermal equilibrium. A rigorous local thermal non-equilibrium (LTNE) assessment would require explicit evaluation of interfacial heat transfer coefficients, phase temperature differences, and thermal relaxation times, which is beyond the scope of the present study. Nevertheless, the present LTE framework is considered sufficient for capturing the dominant macroscopic thermal behavior under the investigated operating conditions. Under the LTE condition, the effective thermophysical properties are defined based on volume-weighted averages. The solid phase material properties are: thermal conductivity ks = 0.21 W/(m·K), density ρs = 1540 kg/m3, and specific heat capacity cp,s = 2.95 × 103 J/(kg·K). The fluid air properties are obtained from the COMSOL built-in material library, with a thermal conductivity of approximately kf = 0.026 W/(m·K) at ambient conditions.

(ρcp)eff=ερfcp,f+(1−ε)ρscp,s(8)

keff=εkf+(1−ε)ks(9)

This approach significantly reduces computational complexity while maintaining acceptable accuracy for the investigated operating conditions. Rather than claiming the LTE assumption as rigorously proven, it is adopted as a practical engineering approximation suitable for the low-velocity, steady-state conditions examined herein. Under substantially higher heating intensities or transient operating conditions, LTNE effects may become prominent and should be considered in future model refinements.

Regarding boundary conditions, a pressure inlet is specified with a reference pressure of 1 atm, whereas various outlet velocities ranging from 0.1 to 0.4 m/s are imposed to investigate the influence of flow rate on transport characteristics. The outer walls of the computational domain are subjected to convective boundary conditions with a heat transfer coefficient of h = 5 W/(m2·K) and an ambient temperature of T∞ = 293.15 K, expressed as:

−k∂T∂n=h(T−T∞)(10)

where h is the external convective heat transfer coefficient and T∞ is the ambient temperature. The heating effect is applied as a boundary heat source on the surface of the heating element, with the thermal power specified as P (1, 2, or 3 W). This boundary heat source is implemented independently of the convective condition on the outer walls, thereby ensuring physical consistency between the power input and heat dissipation to the surroundings. The initial conditions assume a uniform temperature equal to the ambient temperature (293.15 K) and zero velocity throughout the entire domain.

At the interface between the porous region and the free-flow region, continuity conditions for velocity, pressure, and heat flux are enforced to ensure the physical consistency of the coupled solution.

3  Numerical Solution Method

3.1 Geometry Model

The numerical simulations were carried out using COMSOL Multiphysics, where the coupled governing equations were solved based on the finite element method (FEM) [14]. To accurately capture the structural influence on transport behavior, an explicit structure–porous media hybrid modeling approach was developed.

In this approach, the multi-layered porous segment is represented by an assembly of 129 discrete cuboid sheet elements embedded within a continuous fluid domain. Each sheet element functions as an individual porous medium, and the ensemble of 129 sheets is arranged in a concentric multi-ring configuration to form the cylindrical porous segment. Specifically, the sheets are distributed across multiple radial layers, with the spacing between adjacent sheets determined by a stochastic placement algorithm that ensures non-overlapping conditions while preserving a target solid volume fraction of 71.7 percent. The interstitial regions between the sheets are treated as free-flow channels. This hybrid representation enables the model to simultaneously resolve macroscopic flow pathways at the inter-sheet scale and account for microscopic resistance effects within each porous sheet through the Darcy-Forchheimer source terms.

The geometric configuration is based on the characteristic dimensions of the physical system, which is schematically shown along with its mesh generation in Fig. 1. The porous segment has a length of 12 mm, and each sheet element has a cross-sectional size of 1 mm × 0.23 mm. To mimic realistic structural randomness, a stochastic placement algorithm is employed to distribute the sheet elements within the domain, ensuring non-overlapping conditions and a target solid volume fraction of 71.7%. The inlet velocity range of v = 0.1–0.4 m/s coincides with the operating conditions of the target thermal management device, in which the air flow rate through the porous segment is governed by user-determined puffing parameters. These velocities translate to volume flow rates of approximately 0.2–1.0 L/min through the cylindrical segment cross-section of 7.25 mm diameter, which are representative of the actual device operation.

Figure 1: Structure and mesh diagram of the porous segment. (a) Cross-section of the porous segment; (b) Isometric view of the porous segment.

Compared with fully homogenized models, this explicit representation preserves key geometric features such as preferential flow channels and local constrictions, which are essential for capturing flow channeling effects and non-uniform heat transfer behavior.

3.2 Microstructural Basis and Parameter Selection

To ensure the physical realism of the porous media representation, the microstructure of the sheet material was characterized using scanning electron microscopy (SEM), as shown in Fig. 2 [28]. The observations reveal a compact fibrous structure with irregularly distributed pores and partially collapsed voids, indicating a relatively low effective porosity.

Figure 2: SEM image of the microstructure of the porous sheet. (a) Powdered flakes (b) The front and back of the porous thin slice.

Due to the inherent heterogeneity and the coexistence of open and closed pores, direct measurement of porosity is subject to significant uncertainty. Therefore, a parametric porosity range ε = 0.1–0.3 is adopted in the simulations to cover the plausible physical conditions. This range allows for a sensitivity analysis of structural effects and facilitates comparison with experimental measurements to identify the effective porosity.

The permeability α and inertial resistance coefficient Cf are set to α = 1 × 10−12 m2 and Cf = 0.55. The geometric arrangement of the 129 sheet elements remains unchanged for all investigated cases. The porosity variation, spanning ε from 0.1 to 0.3, refers to the intrinsic porosity of the porous sheet material, which affects the effective permeability and thermal properties within each porous sheet domain. The overall geometric packing configuration and the macroscopic sheet distribution are preserved throughout the simulations to isolate the influence of intrinsic porous material properties on the coupled flow and heat transfer behavior. In the COMSOL implementation, the porosity ε is specified directly as a material parameter within the porous medium domain feature for each sheet element, without any modification to the geometric arrangement or the macroscopic packing configuration.

3.3 Meshing and Grid Independence Study

A hybrid meshing strategy is employed to balance computational efficiency and numerical accuracy. The porous sheet elements are discretized using structured hexahedral meshes generated through a sweep operation, while the fluid regions are meshed using unstructured tetrahedral elements. To accurately resolve near-wall gradients, boundary layer meshes are applied to both the outer wall and the surface of the heating element.

To ensure grid-independent solutions, simulations are performed using three levels of mesh density, with total element numbers of 4.77, 6.10, and 9.80 million cells, respectively. The relative variation in flow resistance between consecutive mesh refinements remains below 2%, while the average wall temperature variation is below 0.25%. These results indicate satisfactory grid-independent behavior for the present simulations.

The results are summarized in Table 2. The variation in flow resistance between the 4.77 and 6.10 million-cell meshes is less than 2%, while the deviation in average wall temperature is within 1 K, corresponding to approximately 0.25% of the mean wall temperature. Further refinement to 9.80 million cells produces negligible changes in both metrics. Considering both computational cost and solution accuracy, the mesh with 4.77 million cells is selected for subsequent simulations.

4  Results and Discussion

4.1 Influencing Factors and Regularities of Flow Resistance

To quantitatively evaluate the impact of the explicit porous structure on flow resistance and to distinguish the viscous Darcy contributions from inertial Forchheimer effects, steady-state simulations were conducted at inlet velocities v ranging from 0.1 to 0.4 m/s for varying porosities (ε = 0.1, 0.2, and 0.3). The resulting pressure drops are compiled in Table 3 and illustrated in Fig. 3. To test the linearity of the flow behavior under the investigated velocity range, a least-squares linear regression forced through the origin was performed for each porosity case:

ΔP=a⋅v(11)

the regression slope a is computed as:

a=∑i=1n(vi⋅ΔPi)∑i=1nvi2(12)

and the coefficient of determination R2 is evaluated as:

R2=1−∑i=1n(ΔPi−a⋅vi)2∑i=1n(ΔPi−ΔP¯)2(13)

where ΔP¯ denotes the mean value of the simulated pressure drops. The regression results, including the fitted slopes and R2 values, are summarized in Table 3 alongside the primary simulation data.

Figure 3: Flow resistance vs. outlet flow velocity under different porosities.

The regression analysis yields the following results. For ε = 0.2 and 0.3, the linear regressions yield ΔP = 0.795v (kPa) and ΔP = 0.796v (kPa), respectively, both with an R2 value exceeding 0.999. The near-unity R2 values confirm that the pressure drop is strictly dominated by the linear Darcy viscous term, whereas the Forchheimer inertial contributions remain negligible within the investigated velocity range for these two porosity levels. For ε = 0.1, the linear regression gives ΔP = 0.727v (kPa) with an R2 of 0.963. This lower R2 value quantitatively indicates the onset of nonlinear inertial effects at higher velocities, which is consistent with the Forchheimer-extended Darcy model where the total pressure drop comprises both a linear viscous term proportional to v and a quadratic inertial term proportional to v2. The transition from viscous-dominated to inertia-influenced behavior for ε = 0.1 becomes prominent at v > 0.3 m/s.

When the inlet velocity is low (i.e., v ≤ 0.2 m/s), the flow resistance decreases as the porosity increases. However, when the velocity further increases (i.e., v ≥ 0.3 m/s), an inversion occurs whereby the flow resistance instead increases with increasing porosity. To elucidate the underlying mechanism of this phenomenon, the cross-section at the center of the porous segment was extracted to observe the magnitude of internal vorticity. The vorticity distributions at flow velocities of 0.2 and 0.4 m/s under different porosities are depicted in Fig. 4.

Figure 4: Velocity magnitude distribution within the YOZ central cross-section under different porosities and inlet velocities. (a) ε = 0.1, v = 0.1 m/s; (b) ε = 0.2, v = 0.1 m/s; (c) ε = 0.3, v = 0.1 m/s; (d) ε = 0.1, v = 0.2 m/s; (e) ε = 0.2, v = 0.2 m/s; (f) ε = 0.3, v = 0.2 m/s; (g) ε = 0.1, v = 0.3 m/s; (h) ε = 0.2, v = 0.3 m/s; (i) ε = 0.3, v = 0.3 m/s; (j) ε = 0.1, v = 0.4 m/s; (k) ε = 0.2, v = 0.4 m/s; (l) ε = 0.3, v = 0.4 m/s.

Fig. 4 illustrates the velocity streamlines and contours on the axial cross-section for three porosity values (ε = 0.1, 0.2, and 0.3) at four nominal velocities ranging from 0.1 to 0.4 m/s. Given that the 129-sheet spatial arrangement and the intrinsic Darcy permeability (K) of the porous matrix are held strictly constant across all cases, the explicit structure modeling uniquely isolates the geometric effect of fluid void fraction (ε) on the local flow fields. A significant channeling effect is observed, where the fluid primarily bypasses the porous sheets to flow through the open gaps between them. These core bypass regions exhibit highly elevated flow velocities, as indicated by the high-velocity zones (red areas) in the contours. The present study primarily provides a qualitative visualization of this channeling effect within the explicit structural domain, whereas the quantitative characterization of flow maldistribution, velocity non-uniformity, and preferential pathways will be investigated in future work. Nonetheless, this channeling effect exerts a direct impact on local heat transfer: high-velocity channels drastically enhance convective cooling in their immediate vicinity, whereas the stagnant or lower-velocity zones sheltered between adjacent sheets suffer from reduced heat removal, potentially inducing local hot spots as evidenced by the temperature stratification discussed in Section 4.2.

Regarding the hydrodynamic characteristics, the pressure drop for ε = 0.1 exhibits a discernible deviation from the linear trend observed in the higher-porosity cases. This is reflected by its lower coefficient of determination (R2 = 0.963) compared with the near-unity values (R2 > 0.999) for ε = 0.2 and 0.3. At low velocities (v ≤ 0.2 m/s), the identical pressure drops across all cases confirm that the system is governed by the same linear viscous regime, dictated by the identical Darcy permeability and fixed sheet layout. However, under high-flow conditions, the lower porosity (ε = 0.1) yields a significantly higher interstitial pore velocity (v/ε), triggering an accelerated onset of microscopic inertial resistance and local flow rejection at the fluid-porous interfaces. This results in a sub-linear macroeconomic deviation, where the fluid is forced into highly streamlined bypass pathways that paradoxically cap the peak local velocity at 6.72 m/s and suppress the overall macroscopic pressure drop at 0.4 m/s. In contrast, for ε = 0.2 and 0.3, the larger void fraction mitigates the interstitial velocity escalation, maintaining an excellent linear relationship (R2 > 0.999) throughout the investigated range. Notably, the unchoked interfacial momentum exchange between the open channels and the high-porosity matrix induces stronger localized jetting within the main gaps. This shifts the peak local velocities to 7.91 and 7.92 m/s for ε = 0.2 and 0.3, respectively, resulting in a higher absolute flow resistance at 0.4 m/s while strictly preserving the macroscopic linearity.

4.2 Influencing Factors and Regularities of Temperature Distribution

The heating power P is incorporated into the thermal model as a boundary heat source (Qb) applied to the surface of the central element. This boundary condition satisfies the relation—n·q = Qb, where n denotes the outward unit normal vector, q represents the conductive heat flux vector, and the boundary heat source intensity is defined as Qb = P/A, with A representing the effective heating surface area of the central element. To investigate the comprehensive effects of heating intensity and structural characteristics on the heat transfer performance of the porous segment, combined operating conditions with heating powers of P = 1, 2, and 3 W and inlet velocities ranging from v = 0.1 to 0.4 m/s were simulated at a fixed porosity of ε = 0.2. Additionally, the influence of different porosities (ε = 0.1, 0.2, and 0.3) on the volumetric average temperature was further examined at a heating power of 1 W and an inlet velocity of 0.1 m/s. The numerical results indicate that the average temperature of the porous segment is governed by the coupled effects of these three factors. Specifically, the average temperature increases significantly with the rise in heating power, while it gradually decreases as the inlet velocity increases. The influence of porosity on the overall average temperature is relatively minor; however, it distinctively alters the local thermal distribution patterns.

As shown in Fig. 5, under identical baseline operating conditions, specifically at a heating power P = 1 W and an inlet velocity v = 0.1 m/s, increasing the porosity ε from 0.1 to 0.3 causes the volumetric average temperature of the porous segment to drop marginally from 400.56 to 399.05 K, representing a variation of less than 0.4%. Conversely, at a fixed ε = 0.2, tripling P from 1 to 3 W drives a dramatic temperature surge from 399.78 to 607.62 K. This contrast indicates that under low-velocity and low-power conditions, the global thermal state is remarkably insensitive to variations in ε. Within the explicit structural domain, changes in ε primarily reshape local fluid flow pathways and thermal flux distributions, whereas the macro-scale heat transfer rate remains fundamentally governed by the broader energy balance between total power input and bulk convective transport.

Figure 5: Variation of volumetric average temperature with inlet velocity under different heating powers (ε = 0.2).

This dominant mechanism is further elucidated by the parametric trends presented in Fig. 5. At any given v, the average temperature scales approximately linearly with P. For instance, at a low velocity v = 0.1 m/s, the temperature climbs from 399.78 to 607.62 K, marking a substantial rise of over 50%. Even when the flow is accelerated to a higher velocity v = 0.4 m/s, a similar increase in P still yields a 14% temperature rise, from 321.48 to 367.33 K. These results confirm that P serves as the primary driver of the system’s thermal state. Concurrently, increasing v effectively mitigates this thermal build-up; under a constant heating load P = 2 W, accelerating v from 0.1 to 0.4 m/s reduces the temperature from 504.99 to 349.11 K. This 30% reduction demonstrates the high efficiency of forced convection in suppressing temperature rise within the solid matrix.

Fig. 6 visualizes these coupled effects through the spatial temperature profiles across the central cross-section of the porous segment. Reading the contours from left to right highlights the velocity-dependent behavior as v increases step-wise from 0.1 to 0.4 m/s under a uniform P. Meanwhile, moving from top to bottom tracks the thermal evolution as P scales up from 1 to 3 W at a constant v. At low inlet velocities, the prolonged fluid residence time between adjacent sheets facilitates thorough solid-to-fluid heat exchange, allowing the local domains to closely approach a state of local thermal equilibrium (LTE). In contrast, under high-velocity conditions, the convective transport timescale becomes significantly shorter than the thermal diffusion timescale. The fluid sweeps through the intra-sheet gaps too rapidly to achieve thermal saturation, resulting in a lower global temperature field but creating a highly intensified, steep temperature gradient concentrated immediately around the heating element.

Figure 6: Temperature distribution of the central cross-section under different inlet velocities and heating powers (ε = 0.2).

Fig. 7 illustrates the radial temperature profiles across the porous segment under various combinations of inlet velocity v and heating power P. As observed, the temperature monotonically decreases with increasing radial distance r. At a constant v, elevating P induces a global upward shift of the thermal profiles, reaffirming that the overall thermal magnitude of the system scales directly with the power input. Conversely, under a uniform P, higher values of v drive a distinct downward shift of the curves, demonstrating that intensified convective heat transfer effectively cools the solid matrix. Notably, at a high velocity of v = 0.4 m/s, the localized thermal gradient adjacent to the heating core is significantly steeper than that observed at v = 0.1 m/s. This behavior perfectly aligns with the previously discussed phenomenon where high-speed flows restrict heat penetration and drastically intensify the temperature gradient in the immediate vicinity of the heating element.

Figure 7: Radial temperature distribution curves under different velocities and heating powers (porosity = 0.2).

To further analyze the heat transfer characteristics of individual sheets, a statistical analysis was conducted under baseline operating conditions, namely a porosity ε = 0.2, an inlet velocity v = 0.3 m/s, and a heating power P = 3 W. Each sheet within the porous segment was numbered sequentially to track its average temperature, as illustrated in Fig. 8. Proceeding radially inward from the outer ring to the inner ring, the average temperature of the sheets exhibits a gradual upward trend, reflecting the dominant conductive thermal attenuation expanding outward from the central heating source. Notably, a distinct temperature discontinuity emerges between sheet numbers approximately 100 and 110. This localized step-change coincides with a structural region characterized by wider inter-sheet gaps, where intensified local convective heat transfer acts as a dynamic thermal barrier, disrupting downstream heat propagation and inducing the observed temperature mismatch. Although this channeling-induced flow locally mitigates peak temperatures, it simultaneously exacerbates spatial thermal non-uniformity, which can be highly detrimental in applications requiring strict temperature homogeneity.

Figure 8: Volumetric average temperature distribution of each sheet (ε = 0.2, v = 0.3 m/s, P = 3 W). (a) Sheets number distribution; (b) Volumetric average temperature distribution of sheets.

4.3 Experimental Verification

To guarantee experimental reliability, each volumetric flow rate condition was evaluated across three independent experimental runs. The corresponding individual run data, calculated mean values, and standard deviations are summarized in Table 4, with the final flow resistance values represented by the arithmetic mean of these triplicate measurements.

As illustrated in Fig. 9, a micro differential pressure transmitter is integrated with a data acquisition system to monitor the flow resistance across the porous segment in real time. To facilitate a direct correlation with the theoretical framework, the inlet velocity v defined in the preceding numerical simulations was converted into the volumetric flow rate Q via the relationship Q = v·Ac, where the cross-sectional area Ac equals 41.25 mm2 for the cylindrical porous segment with a diameter of 7.25 mm. Accordingly, the investigated velocity range from 0.1 to 0.4 m/s corresponds to a volumetric flow rate spanning from 0.2475 to 0.99 L/min. To enable a direct comparison with the numerical predictions from Fig. 3, both the volumetric flow rate Q and its corresponding inlet velocity v are explicitly mapped on the axes of Fig. 10. The experimental porous segment features an identical geometric configuration to the numerical model, maintaining a length of 12 mm, and is enclosed within a glass sleeve to replicate the specified thermal boundary conditions.

Figure 9: Test bench. (a) Diagram of the test bench; (b) Test bench.

Figure 10: Comparison of simulated and experimental flow resistance. The bottom horizontal axis shows volume flow rate (L/min), and the top horizontal axis shows the corresponding inlet velocity (m/s) for direct comparison with Fig. 3.

Uncertainty Analysis

To guarantee the reliability of the experimental data, a rigorous uncertainty analysis was performed in accordance with the Guide to the Expression of Uncertainty in Measurement (GUM). The experimental system involves two types of uncertainties: Type A uncertainty arising from random errors in repeated measurements, and Type B uncertainty derived from the inherent precision limits of the measuring instruments. The flow resistance is directly characterized by the pressure drop ΔP across the porous sample. The micro differential pressure transmitter has a full scale (FS) of 5000 Pa with an accuracy of ±0.5%  FS, yielding a maximum permissible error of 25 Pa. Assuming a uniform probability distribution, the Type B standard uncertainty of the pressure drop is uB(ΔP)=25/3=14.43 Pa. The suction velocity v is calculated indirectly from the volumetric flow rate Q measured by the mass flow meter and the cross-sectional area of the sample (Ac = 41.25 mm2), defined as v = Q/Ac. The mass flow meter has a full scale of 5 L/min with a precision of ±0.2%  FS, yielding a maximum flow rate error of 0.1 L/min. The corresponding Type B standard uncertainty for the flow rate is uB(Q)=0.1/3=0.0577 L/min.

Each experimental case was repeated three times to mitigate random errors, and the Type A standard uncertainty uA was evaluated from the statistical standard deviation of the repeated trials. By applying the error propagation law, the combined relative standard uncertainties of the measured flow resistance and the suction velocity were determined. Under the tested conditions, the combined relative uncertainty of the volumetric flow rate ranges from 5.83% at the maximum flow rate of 0.99 L/min to 23.31% at the minimum flow rate of 0.2475 L/min. For the experimental pressure drop, the combined relative uncertainty ranges from 4.65% at the maximum pressure drop of 0.31 kPa to 18.04% at the minimum pressure drop of 0.08 kPa. The maximum relative uncertainties occur at the lowest suction velocity of 0.1 m/s due to the low absolute readings relative to the instrument scales. As the suction velocity scales up to 0.4 m/s, the relative uncertainty of the flow resistance drops significantly to below 5%. The overall uncertainty levels remain well within acceptable limits for multi-physics and fluid dynamics investigations, confirming the high fidelity of the experimental data and validating its close agreement with the numerical simulation results.

As illustrated in Fig. 10, the measured flow resistance, characterized by the pressure drop ΔP, exhibits exceptional agreement with the numerical predictions across the entire evaluated operating envelope. Specifically, at the maximum volumetric flow rate Q = 0.99 L/min, which corresponds to a suction velocity v = 0.4 m/s, the experimental ΔP of 0.31 kPa deviates from the numerical prediction of 0.322 kPa by a minor −3.7%. Conversely, at the minimum Q = 0.2475 L/min, where v scales down to 0.1 m/s, the measured ΔP of 0.08 kPa shows a modest deviation of +4.8% against the predicted value of 0.076 kPa. Both bounding operating conditions are well encapsulated within the uncertainty envelopes established by the GUM-based uncertainty analysis, thereby demonstrating the robust predictive fidelity of the numerical framework under varied hydrodynamic loads. For thorough documentation, the complete triplicate measurement dataset is archived in Table 4, while a comprehensive data cross-comparison is tabulated in Table 5.

In addition, the simulated flow resistance under various porosities ε was cross-referenced with experimental measurements to systematically evaluate the fidelity of the selected structural domain. The comparative analysis indicates that the numerical flow resistance obtained at ε = 0.2 and 0.3 demonstrates robust alignment with the experimental data. Conversely, the case featuring a lower porosity of ε = 0.1 tends to noticeably underestimate the measured pressure drop ΔP. This discrepancy suggests that while the exact microscopic porosity of the individual sheets eludes high-precision direct characterization, the macroscopic effective porosity governing the bulk transport behavior of the porous matrix is highly likely to reside within the bounded interval of 0.2 ≤ ε ≤ 0.3 under the investigated operating envelope.

It must be acknowledged that the current experimental validation is primarily constrained to hydrodynamic flow resistance characteristics. Fully capturing the internal thermal fields would necessitate the embedding of intra-matrix micro-sensors within the porous domain. However, such an implementation presents severe diagnostic challenges due to the stringent spatial constraints of the assembly, which features a miniature diameter of 7.25 mm and hosts 129 densely packed micro-sheets under a localized high-flux heating configuration. Consequently, the thermal distributions presented herein should currently be interpreted as physically informed numerical predictions rather than fully empirical thermal measurements. Future endeavors will focus on incorporating spatially resolved, non-intrusive temperature diagnostic techniques to rigorously assess the predictive accuracy of the conjugating thermal transport model and the localized validity of the local thermal equilibrium (LTE) hypothesis.

5  Conclusion

In this study, a three-dimensional hybrid framework coupling explicit structural representation with porous media theory was developed to systematically investigate the conjugated flow and heat transfer characteristics within multi-layered porous sheet arrays. The model integrates geometric reconstruction with macroscopic transport principles and is rigorously validated against experimental measurements. The primary findings are summarized as follows:

(1) The flow resistance exhibits a predominantly linear dependence on the inlet velocity v under low-flow conditions, indicating viscous-dominated Darcy behavior. However, within low-porosity domains, a distinct transition toward nonlinear Forchheimer behavior emerges at elevated velocities owing to the growing contribution of inertial effects, a phenomenon structurally driven by contracted flow cross-sections and intensified localized acceleration.

(2) The explicit structural representation successfully captures a pronounced channeling effect, wherein the fluid preferentially bypasses narrow constrictions to traverse larger inter-sheet gaps, thereby creating localized high-velocity pathways. These preferential flow channels dominate macroscopic energy dissipation and primarily govern the observed deviations from classical linear Darcy transport laws.

(3) The global thermal state of the porous segment is predominantly governed by the heating power P, exhibiting a near-linear scaling with increasing power input. Conversely, accelerating v significantly suppresses the temperature rise, an effect directly attributed to the intensification of convective heat removal timescales over conductive transport.

(4) Under the baseline low-power and low-velocity operating regime, specifically at P = 1 W and v = 0.1 m/s, the impact of structural porosity ε on the core volumetric average temperature remains marginal across the evaluated domain of 0.1 ≤ ε ≤ 0.3. Nevertheless, under elevated thermal loads or intensified flow rates, the role of ε is expected to become highly pronounced and warrants comprehensive future exploration. Crucially, variations in ε significantly alter the localized flow field topologies, thereby inducing substantial spatial thermal non-uniformity, particularly within regions characterized by structural irregularities.

(5) The numerical predictions demonstrate robust alignment with the experimental flow resistance data. Across the entire evaluated volumetric flow rate envelope spanning from 0.2475 to 0.99 L/min, which corresponds to suction velocities of 0.1 to 0.4 m/s, the relative deviation between the simulated and measured pressure drops remains strictly bounded between −3.7% and +4.8%. This close agreement further confirms that the macroscopic effective porosity of the physical system resides within the range of 0.2 ≤ ε ≤ 0.3, validating the predictive fidelity of the hybrid modeling framework. Although direct thermal validation was circumvented due to diagnostic limitations, the tightly coupled mathematical formulation of the governing fluid-thermal equations provides a physically consistent foundation for the predicted internal temperature fields. Future experimental campaigns will integrate internal thermal sensing to comprehensively validate the heat transport component.

In summary, this work underscores the critical necessity of explicitly resolving microstructural topologies within porous domains to ensure the high-fidelity prediction of coupled thermal-fluid physics. The developed hybrid framework serves as a potent diagnostic and optimization tool for advanced thermal management architectures utilizing intricate porous configurations. Nonetheless, several boundaries of the present study should be acknowledged: (i) the internal heat transfer profiles require dedicated empirical validation via future thermal measurements; (ii) the local thermal equilibrium (LTE) hypothesis, while quantitatively justified under the present steady-state conditions by a low Biot number Bi = 0.068, which is well below the standard threshold of 0.1, may necessitate multi-scale revisions under transient or high-flux scenarios; (iii) the influence of ε on the spatial temperature distribution has been explored only within a restricted parametric subset and should be systematically generalized using a full factorial design matrix; and (iv) the long-term structural stability of the multi-layered assembly under cyclic thermal fatigue represents a vital avenue for downstream exploration. As this research constitutes the foundational systematic investigation into this novel multi-layered porous sheet architecture, subsequent iterations of this modeling framework will expand to encompass mass transport kinetics, cyclic thermal loading hysteresis, and structural degradation mechanics.

Acknowledgement: Not applicable.

Funding Statement: The authors received no specific funding for this study.

Author Contributions: Conceptualization, Kai Luo, Haibing Chen and Wei Tang; methodology, Kai Luo, Shaohua Bi and Weifeng He; software, Kai Luo; validation, Kai Luo, Yifei Xie and Kun Chen; formal analysis, Kai Luo and Jirong Zhang; investigation, Kai Luo, Yifei Xie and Kun Chen; resources, Haibing Chen and Wei Tang; data curation, Kai Luo and Kun Chen; writing—original draft preparation, Kai Luo; writing—review and editing, Haibing Chen, Wei Tang and Shaohua Bi; visualization, Kai Luo; supervision, Haibing Chen and Wei Tang; project administration, Wei Tang; funding acquisition, Haibing Chen and Wei Tang. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data will be made available on request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Zheng KC, Wen Z, Wang ZS, Lou GF, Liu XL, Wu WF. Review on forced convection heat transfer in porous media. Acta Phys Sin. 2012;61(1):014401. doi:10.7498/aps.61.014401. [Google Scholar] [CrossRef]

2. Alabdaly IK. Heat transfer and flow behavior in porous media: a review. Anbar J Eng Sci. 2025;16(2):66–86. doi:10.37649/ajes-25-022. [Google Scholar] [CrossRef]

3. Jin K, Zhao J, Yao G, Xu D, Wen D. A fully coupled multiscale phase-change model at the porous interface for transpiration cooling: coupling dynamics pore-scale networks to continuum-scale free flow. Int J Heat Mass Transf. 2025;236(10):126396. doi:10.1016/j.ijheatmasstransfer.2024.126396. [Google Scholar] [CrossRef]

4. Nield DA, Bejan A. Convection in porous media. 5th ed. Cham, Switzerland: Springer International Publishing; 2017. doi:10.1007/978-3-319-49562-0. [Google Scholar] [CrossRef]

5. Li Q, Wang J, Li J, Shi J. Fundamental numerical analysis of a porous micro-combustor filled with alumina spheres: pore-scale vs. volume-averaged models. Appl Sci. 2021;11(16):7496. doi:10.3390/app11167496. [Google Scholar] [CrossRef]

6. Chu X, Yang G, Pandey S, Weigand B. Direct numerical simulation of convective heat transfer in porous media. Int J Heat Mass Transf. 2019;133(1):11–20. doi:10.1016/j.ijheatmasstransfer.2018.11.172. [Google Scholar] [CrossRef]

7. Zhang L, Guo B, Qin C, Xiong Y. A hybrid pore-network-continuum modeling framework for flow and transport in 3D digital images of porous media. Adv Water Resour. 2024;190(2):104753. doi:10.1016/j.advwatres.2024.104753. [Google Scholar] [CrossRef]

8. García-Salaberri PA, Zenyuk IV. A general-purpose tool for modeling multifunctional thin porous media (POREnetfrom pore network to effective property tensors. Heliyon. 2024;10(4):e26253. doi:10.1016/j.heliyon.2024.e26253. [Google Scholar] [PubMed] [CrossRef]

9. Simonov OA, Erina YY, Ponomarev AA. Review of modern models of porous media for numerical simulation of fluid flows. Heliyon. 2023;9(12):e22292. doi:10.1016/j.heliyon.2023.e22292. [Google Scholar] [PubMed] [CrossRef]

10. Meinicke S, Dubil K, Wetzel T, Dietrich B. Characterization of heat transfer in consolidated, highly porous media using a hybrid-scale CFD approach. Int J Heat Mass Transf. 2020;149:119201. doi:10.1016/j.ijheatmasstransfer.2019.119201. [Google Scholar] [CrossRef]

11. Ma R, Wang Z, Wang L, Zhang Z, Li S, Liu B. Pore-scale flow and conjugate heat transfer in heterogeneous porous Media: vortex dynamics and thermal management implications. Therm Sci Eng Prog. 2025;68(23–24):104378. doi:10.1016/j.tsep.2025.104378. [Google Scholar] [CrossRef]

12. Yang M, Wang D, Dong Z, Zhang Z, Li Y, Du D. Insight into the permeability effect on forced convective heat transfer characteristics in porous media based on the pore-scale numerical study. Int Commun Heat Mass Transf. 2025;165(1):109004. doi:10.1016/j.icheatmasstransfer.2025.109004. [Google Scholar] [CrossRef]

13. Trilok G, Gnanasekaran N, Mobedi M. Reconciliation of wire woven mesh porous media in controlling the desired heat transfer and pressure drop. Int J Numer Meth Heat Fluid Flow. 2025;35(3):981–1011. doi:10.1108/hff-01-2024-0014. [Google Scholar] [CrossRef]

14. Zhang B, Huang Y, Cui P, Wang Z, Ding D, Pan Z, et al. Enhanced evaporation of ternary mixtures in porous medium with microcolumn configuration. Front Heat Mass Transf. 2024;22(4):997–1016. doi:10.32604/fhmt.2024.053592. [Google Scholar] [CrossRef]

15. Ding H, Wang J, Zhang Y, Yan H. Analytical model for transient transport of organic contaminants in non-isothermal fractured porous media. Comput Geotech. 2025;187(4):107458. doi:10.1016/j.compgeo.2025.107458. [Google Scholar] [CrossRef]

16. Rehman MAU, Farooq MA, Hassan AM. Thermal radiation effects on 2D stagnation point flow of a heated stretchable sheet with variable viscosity and MHD in a porous medium. Front Heat Mass Transf. 2024;22(1):263–86. doi:10.32604/fhmt.2023.044587. [Google Scholar] [CrossRef]

17. He Y, Chen Y, Du S, Li D, Gao Y. Optimization studies of transpiration cooling using porous medium with gradually-changed structure. Front Heat Mass Transf. 2020;15(1):1–13. doi:10.5098/hmt.15.19. [Google Scholar] [CrossRef]

18. Wang S, Wang H, Zhang Y, Zhang Z, Jia L. Experimental study on the flow boiling of R134a in sintered porous microchannels. Front Heat Mass Transf. 2025;23(6):1721–40. doi:10.32604/fhmt.2025.073226. [Google Scholar] [CrossRef]

19. Vijay S, Luhar M. Pressure drop measurements over anisotropic porous substrates in channel flow. Exp Fluids. 2024;65(9):135. doi:10.1007/s00348-024-03873-2. [Google Scholar] [CrossRef]

20. Lin S, Shi F, Feng W. Experimental and numerical simulation research on fluid pressure drop inside sintered metal porous plates. J Phys Conf Ser. 2024;2775(1):012022. doi:10.1088/1742-6596/2775/1/012022. [Google Scholar] [CrossRef]

21. Marykovskiy Y, Pomaranzi G, Schito P, Zasso A. A method to evaluate forchheimer resistance coefficients for permeable screens and air louvers modelled as a porous medium. Fluids. 2024;9(7):147. doi:10.3390/fluids9070147. [Google Scholar] [CrossRef]

22. Zhang Z, Zhang X, Duan J, Liu Y, Yan G, Sun M, et al. Applicability of forced convection heat transfer topology optimization methods at the cell scale and heat exchanger core scale. Appl Therm Eng. 2025;281(2):128709. doi:10.1016/j.applthermaleng.2025.128709. [Google Scholar] [CrossRef]

23. Al-Sumaily GF, Al Ezzi A, Dhahad HA, Thompson MC, Yusaf T. Legitimacy of the local thermal equilibrium hypothesis in porous media: a comprehensive review. Energies. 2021;14(23):8114. doi:10.3390/en14238114. [Google Scholar] [CrossRef]

24. Andrade JS Jr, Costa UM, Almeida MP, Makse HA, Stanley HE. Inertial effects on fluid flow through disordered porous media. Phys Rev Lett. 1999;82(26):5249–52. doi:10.1103/physrevlett.82.5249. [Google Scholar] [CrossRef]

25. Sivanesapillai R, Steeb H, Hartmaier A. Transition of effective hydraulic properties from low to high Reynolds number flow in porous media. Geophys Res Lett. 2014;41(14):4920–8. doi:10.1002/2014gl060232. [Google Scholar] [CrossRef]

26. Yang SM, Tao WQ. Heat transfer. 4th ed. Beijing, China: Higher Education Press; 2006. (In Chinese). [Google Scholar]

27. Kaviany M. Principles of heat transfer in porous media. 2nd ed. Berlin/Heidelberg, Germany: Springer; 2012. doi:10.1007/978-1-4612-4254-3. [Google Scholar] [CrossRef]

28. Kerdudou P, Chéné JB, Jacqus G, Perrot C, Berger S, Leroy P. A semi-empirical approach to link macroscopic parameters to microstructure of fibrous materials. In: Proceedings of the 44th International Congress and Exposition on Noise Control Engineering (Inter-Noise 2015); 2015 Aug 9–12; San Francisco, CA, USA. p. 5468–79. [Google Scholar]