Effects of Microscopic Properties on Macroscopic Thermal Conductivity for Convective Heat Transfer in Porous Materials

Porous materials are widely used in many heat transfer applications. Modeling porous materials at the microscopic level can accurately incorporate the detailed structure and substance parameters and thus provides valuable information for the complex heat transfer processes in such media. In this study, we use the generalized periodic boundary condition for pore-scale simulations of thermal flows in porous materials. A two-dimensional porous model consisting of circular solid domains is considered, and comprehensive simulations are performed to study the influences on macroscopic thermal conductivity from several microscopic system parameters, including the porosity, Reynolds number, and periodic unit aspect ratio and the thermal conductance at the solid–fluid interface. Our results show that, even at the same porosity and Reynolds number, the aspect ratio of the periodic unit and the interfacial thermal conductance can significantly affect the macroscopic thermal behaviors of porous materials. Qualitative analysis is also provided to relate the apparent thermal conductivity to the complex flow and temperature distributions in the microscopic porous structure. The method, findings and discussions presented in this paper could be useful for fundamental studies, material development, and engineering applications of porous thermal flow systems.


Introduction
Fluid flows and the associated heat transfer in porous materials have attracted great interest over the decades for their scientific and practical importance. With the large contact area between the solid matrix and the coolant fluid, heat transfer performance is significantly enhanced in such porous materials. Relevant applications can be found in various industrial areas, including catalytic beds for chemical reactions, water treatment, nanofluids, heat exchangers, heat sinks and thermal energy storage systems [1][2][3][4][5][6][7][8]. Extensive experimental and theoretical investigations have been conducted to study the heat transfer performance in porous materials [9][10][11]. In addition to the high cost and long time duration, experimental studies are not able to reveal the flow and temperature distributions in the complex microscopic porous geometry, while such information is crucial for understanding the mechanism and relationships among various factors involved. On the other hand, theoretical analysis is limited to simple geometry structures and idealized flow situations, and such conditions can hardly be satisfied in realistic systems. Fortunately, with the advances in computer and software technologies, numerical simulations have been proved to be useful for various complex systems, including flows and heat transfer in porous media [12][13][14][15][16]. In the literature, there exist two typical approaches for modeling the flow and heat transfer in porous media: the continuum and the pore-scale approaches. The first method treats the porous materials as continuum media and solves the macroscopic flow and energy equations with apparent parameters [17,18]. Despite the model's simplicity and computational efficiency, this approach relies on the accuracy of those apparent parameters and also cannot incorporate the specific microscopic structure and material thermophysical properties in the porous media. On the other hand, the pore-scale method solves the flow and energy equations with the microscopic structure and properties of the matrix material considered explicitly [19,20]. Pore-scale simulations can provide detailed flow and thermal fields at the microscopic level, and apparent properties can be obtained from these microscopic distributions for macroscopic analysis [21][22][23]. In recent years, the lattice Boltzmann method (LBM) has become especially popular in such simulations, mainly for its convenience in dealing with complex boundary geometry [22][23][24][25].
With the microscopic porous structure modeled explicitly, the physical scale of porescale simulations is limited to small sizes due to the large computation demand. For this reason, pore-scale simulations typically work with a small unit and assume that such identical units are repeated in space. To solve the flow and energy equations in a periodic unit, appropriate boundary conditions are necessary on the side surfaces of this computational unit. Such side surfaces are simply for computational convenience and they are not physical surfaces or boundaries; therefore, definite conditions on such surfaces are not available. To create a thermal gradient, previous studies are usually assigned constant temperatures on two opposite unit surfaces and assumed adiabatic (no heat flux) conditions on other side walls [24][25][26][27]. This treatment has been examined recently and significant concerns have been raised in terms of its validity and accuracy [28]. Few studies treated the unit surface conditions more reasonably by considering the periodic relations of temperature, however, the temperature or heat flux on the fluid-solid interface was fixed at constant values [29][30][31]. Clearly this is physically not true: the thermal condition (temperature and heat flux) at the microscopic fluid-solid interface is determined by the particular local flow and thermal situations and it cannot be specified in advance. The interfacial thermal resistance may also exist at the fluid-solid interface in porous materials [32] or the constituent interface in composite materials [33], and this feature has typically been neglected in previous studies. Furthermore, previous studies usually only considered situations with the flow along the temperature gradient direction [22,25,30], whereas in many experimental setups and industrial applications the temperature gradient is mainly in the orthogonal plane to the fluid flow [34][35][36].
In this paper, we study the convective heat transfer in porous materials using the recent generalized periodic boundary method for thermal flows by Jbeili and Zhang [28]. This boundary method was established with rigorous mathematical justifications and validated by carefully designed numerical tests. Using this method, extensive simulations are conducted based on a two-dimensional (2D) porous model of circular solid particles. Unlike previous studies usually focusing on the porosity and Reynolds number effects, we also investigate the influences of the unit aspect ratio and the interface conductance on the macroscopic thermal performances of porous flows. Our results show that, even with the same porosity and Reynolds number, the unit aspect ratio can greatly affect the apparent thermal conductivity at the microscopic level. It is also interesting to observe that a weaker interfacial conductance can reduce the tortuosity conductivity but increase the dispersion conductivity, and that the dispersion conductivity can respond to the porosity change in a non-monotonic fashion. In addition, qualitative discussions are provided to relate the macroscopic conductivity coefficients to the microscopic flow and thermal fields for a better understanding of the complex nature of porous thermal flows.

Governing Equations and Boundary Conditions for Periodic Unit
For simplicity and clarity, in this study we consider the 2D porous model with circular solid domains as shown in Figure 1a. The fluid flows vertically from the bottom to the top, while the thermal gradient is imposed in the horizontal direction. The governing equations for this thermal flow system include the continuity equation, Equation (1), and the Navier-Stokes equation Equation (2) for the fluid domain Ω 0 : ρ ∂u ∂t + u · ∇u = −∇p + µ∇ 2 u + f; (2) and the energy equation for the fluid (Ω 0 ) and solid (Ω 1 ) domains, respectively [37]. In these equations, ρ represents the fluid density, u is the fluid velocity, p is the pressure, µ is the dynamic viscosity of fluid, f is the body force on the fluid, T is the temperature, and α is the thermal diffusivity. The subscripts f and s for T and α denote the fluid and solid domains respectively. On the fluid-solid interface Γ, the no-slip boundary condition is applied for the fluid flow to incorporate the possible thermal resistance and therefore the temperature discontinuity at the solid-fluid interface Γ, a general conjugate condition for temperature is given as [38][39][40]: Here, q is the heat flux along the local normal direction n, which points from Domain Ω 1 toward Domain Ω 0 ( Figure 1). C is the interface conductance, and K f and K s are the thermal conductivities of the fluid and solid substances, respectively. In the pore-scale approach, only one periodic unit of the porous material is considered in the simulation (Figure 1). To solve the governing equations in this periodic unit, appropriate boundary conditions are required for the side boundaries of the periodic unit. These side boundaries are virtual but not real surfaces, and thus physical constraints on such boundaries are not directly available. For flows in such periodic structures with negligible fluid property changes, it has been well recognized that the flow field would be identical in each periodic unit [21,[41][42][43]: where l and h are the unit length and height (Figure 1), and m and n are arbitrary integers. Based on previous practices in periodic thermal flow simulations [28,[41][42][43][44][45], the following relationship is proposed for the temperature field among periodic units in Figure 1a: T(x + ml, y + nh) = T(x, y) + mlT g .
Here, T g is the global thermal gradient in the horizontal direction, and it can be obtained from the temperature difference T H − T L and the distance L between the locations where these two temperatures T L and T H are imposed (see Figure 1a) as T g = (T H − T L )/L. Assume u(x, y) and T(x, y) are the correct flow and temperature solutions in one periodic unit (0 ≤ x ≤ l and 0 ≤ y ≤ h), it can be readily shown that the velocity from Equation (6) and temperature from Equation (7) satisfy all the governing equations and boundary conditions given in Equations (1)- (5), and thus they are the correct solutions for the unit of ml ≤ x ≤ (m + 1)l and nh ≤ y ≤ (n + 1)h. Therefore, in pore-scale simulations, one can simulate the flow and thermal fields in one periodic unit as shown in Figure 1b, and the solutions can be extended to other units according to Equations (6) and (7). These relations can then be used to establish correct boundary conditions for the side surfaces of the periodic unit. For example, the right boundary x = l for the unit shown in Figure 1b actually is also the left boundary of the next unit on its right. According to Equations (6) and (7) with m = 1 and n = 0, we have the left-right boundary relations for the periodic unit as: u(l, y) = u(0, y), T(l, y) = T(0, y) + lT g .
Similarly, the top-bottom boundary relations for flow and temperature are expressed as: These relations are similar to those used by Kuwahara et al. [21]; however, no mathematical justifications and numerical validations were provided there. In addition to these boundary relations, due to the linearity and homogeneity of the governing equations and boundary conditions, extra anchoring conditions are necessary to ensure the numerical convergence [39,41,42]. For our current system in Figure 1a, the following conditions are adopted in our next simulations: where U 0 is the mean velocity through the porous medium and T 0 is the mean temperature at the left unit boundary. The mean velocity U 0 is related to the Reynolds number Re of the system (to be defined later), while the mean temperature T 0 can be set at an arbitrary value and it has no impact on the result analysis and interpretation.

Numerical Validation of the Periodic Relations among Units
To further verify the flow and temperature relations among units given in Equations (6) and (7), we conduct a direct numerical simulation for the system in Figure 1a. This simulation is also helpful to illustrate our concerns with the periodic unit boundary treatments used in previous studies. The diameter of the solid domains is D = 50 and they are arranged as a square array with l = h = 2D. The porosity of this model medium is thus = π/16. The rectangular simulation domain has a length of L = 20D in the horizontal direction and a hight of H = 10D in the vertical direction, and thus the simulation domain consists 10 × 5 = 50 identical periodic units as shown in Figure 1b u(0, y) = u(L, y), u(x, 0) = u(x, H), and T(x, 0) = T(x, H). Please note that the temperature relation Equation (7) is not involved in this validation simulation. A body force f in the upward direction is utilized to generate the flow, and its magnitude is adjusted to obtain the Reynolds number Re = ρU 0 √ lh/µ = 50. Here we use the geometric mean value of the periodic unit length l and height h for the Reynolds number Re, and this choice is made for our convenience in examining the effect of the aspect ratio β = h/l in Section 4.1. The Prandtl number of the fluid is 0.7. The thermal conductivities are set as K f = 0.2 for fluid and K s = 10K f = 2 for solid. A relatively large interface conductance value C = 10 is adopted to minimize the temperature discontinuity across the interface. The lattice Boltzmann method (LBM) with the D2Q9 (2D and with nine lattice velocities) lattice structure is used to solve the flow and thermal fields [44,[46][47][48] for all calculations in this paper, and the counter-extrapolation method [49] is adopted for the conjugate thermal condition on domain interfaces. As in general LBM studies [16,22,31,[48][49][50], all quantities provided above and in the later result presentation are non-dimensional based on LBM simulation units (for example, length in the lattice grid resolution, and time in the simulation time step).
The calculated flow and temperature fields from this direct simulation are displayed in Figure 2. The flow pattern in Figure 2a appears identical in each periodic unit. For the temperature in Figure 2b, it can be seen that the temperature increases in general along the horizontal direction, agreeing with the global thermal gradient generated from the two boundary temperatures T L = 0 at the left and T H = 1 at the right. In Figure 2b, neither the temperature nor the heat flux are constant along the solid-fluid interface (edges of the circular patches), meaning that the constant temperature or flux assumptions for the interfaces in previous studies [22,30,31] are not valid. The wider gaps between isotherm contours in the solid domain imply smaller temperature gradients there, and this is consistent with the larger solid conductivity used in this simulation. With a relatively larger conductance value for the solid-fluid interface, no apparent discontinuity in temperature can be observed across the interfaces. The apparently linear increase in temperature and the similarity in local isotherm contours in Figure 2b appear rational according to the temperature relation Equation (7). Moreover, T ∼ y profiles at six horizontal locations are plotted in Figure 3, one in each panel from the left and with the x positions labeled on top. These horizontal positions are selected at the same relative locations to the nearby solid columns (0.4D to the left of the patch centers). The flat segments on these curves occur in the solid domains, and they are due to the vertical isotherm lines there. The strong variations in these profiles clearly show that in general one should not assign constant temperatures on the periodic unit boundaries, as done in [24][25][26]. According to Equation (7), these temperature profiles should have the same variation features along the vertical direction, except different offset values. This is apparently true by looking at these individual profiles. For a more quantitative confirmation, we then shift each profile by its mean temperature T m , and the six shifted profiles are plotted in the last panel on the right in Figure 3. These six shifted profiles overlap with each other perfectly and we only see one single curve there. These identical shifted profiles precisely confirm that the temperature values at same relative positions in individual periodic units are only different in respective mean temperatures and the variation fashion is identical. Furthermore, the mean temperature values calculated along these six profiles are listed in Table 1. According to Equation (7), these mean temperatures can be related to the mean value at x/D = 0.6 by with x i = 0.6D, 2.6D, 6.6D, 10.6D, 14.6D and 18.6D for the six profiles in Figure 3 (from left to right). The predicted values from this equation are also provided in Table 1 for comparison, and the excellent agreement convincingly confirms the validity of the temperature relation Equation (7) for such porous thermal flows.

Simulation Results and Discussion
The enhanced heat transfer performance in porous flow systems, in principle, benefits from the large solid-fluid contact area and the long, twisted flow path as the fluid passes through the microscopic porous structure. At the macroscopic, practical level, the apparent thermal conductivity is often used to quantify the overall thermal behaviors. The volume averaging analysis [21,[51][52][53] can be applied to obtain the effective conductive tensor from the microscopic flow and thermal fields. Kuwahara et al. [21] further decomposed this conductivity tensor K into three parts where K e is the stagnant conductivity, which is simply the volume average of the fluid and solid conductivities based on the porosity : Please note that the expression for K e in Equation (13) is simply the coefficient for the isotropic part of the conductivity matrix K from the volume average analysis [21,[51][52][53]; and that there are no assumptions involved on the microscopic porous structures. K tor and K dis are called the tortuosity and dispersion conductivity tensors, respectively. When the mean flow direction is along one of the coordinate directions, only diagonal elements are nonzero; and obviously the primary concerns are the diagonal elements in the global thermal gradient direction [21,53,54]. For the setup in Figure 1, we simply use K tor and K dis to denote the xx components of tensors K tor and K dis ; and they can be calculated from the simulated flow and temperature distributions in one periodic unit by [21,53] where c p f is the fluid specific heat, and the volume average temperatureT and the intrinsic average velocityū f are given by: The subscript x denotes the x component of the corresponding vector.
In this section, we apply the boundary conditions in Section 2 to the periodic unit shown in Figure 1b, and examine the effects on the effective macroscopic conductivity coefficients K tor and K dis from several microscopic parameters, including the aspect ratio β = h/l, the Reynolds number Re = ρ √ hlU 0 /µ, the interfacial thermal conductance C, and the porosity = 1 − πD 2 /4hl. We keep the unit area lh = 90, 000 constant in our next simulations, so that the Reynolds number Re is directly proportional to the mean flow velocity U 0 , and the porosity depends purely on the solid diameter D and not affected by the aspect ratio β. The mean temperature at x = 0 is set as T 0 = −0.05 and the global thermal gradient is T g = 0.1/l. Please note that, due to the linearity and homogeneity of the energy equation Equation (3) and the associated thermal boundary conditions, the particular values of T 0 and T g used in the simulations do not affect the calculated conductivities K tor and K dis from Equations (14) and (15). To make this point clear, we consider a periodic unit under two situations: Situation (a) with mean inlet temperature T 0,a and thermal gradient T g,a ; and Situation (b) with mean inlet temperature T 0,b and thermal gradient T g,b . If T a (x, y) is the solution in the periodic unit with Situation (a), the temperature field for Situation  (14) and (15) and considering Γ u x dΓ = 0 for periodic structures, one can find that the K tor and K dis values remain the same for Situations (a) and (b). Other simulation parameters are kept the same as given in Section 3, unless otherwise mentioned. As for the initial state, we start simulations with a linear transition from T 0 at the inlet to T 0 + T g l at the outlet for the temperature, and zero velocity for the flow. The results are taken when steady or quasi-steady states are established. The computer code for this work has been developed based on our programs used in previous publications [28,38,44,49,55,56], and all important elements involved in our simulations, including the LBM algorithms for flow and heat transfer, the no-slip boundary and conjugate interface treatments, as well as the mesh resolution selection, have been validated and confirmed in these previous studies.

Effects of the Aspect Ratio β and Reynolds Number Re
We start with simulations to study the effect of the aspect ratio β on the macroscopic thermal conductivity, which has not been well addressed in previous investigations. Here we set the porosity = 0.85, and accordingly the solid domain diameter is D = 131.11. Following previous studies [29,57,58], two representative Reynolds number values, Re = 50 and 100, are considered in this section. Higher Reynolds numbers are possible for gas flows through porous media [59,60]. Our simulations cover a range of 0.25∼4 for the aspect ratio β; further increasing (>4) or decreasing (<0.25) of the β value will make the gap between the solid surfaces too small for accuracy and stable computations. The flows are always steady for all tested β values at Re = 50 and for β ≤ 1 at Re = 100; however, the flow becomes unsteady for β > 1 and Re = 100. This transition is reasonable, since for a larger β, the gap between the solid surfaces becomes narrower, and the fluid velocity through this gap increases accordingly. This situation is similar to a flow jet coming out from a small opening shooting into a relatively large space. Compared to the condition with the same mean flow velocity (same Reynolds number) but a smaller aspect ratio β, the reduction in flow passage width is less significant (a wider gap between solid surfaces) and the fluid space after the gap is relatively limited. Figure 4 shows the streamlines and temperature distributions for two representative aspect ratios, β = 0.44 and 1.44 at Re = 50 and 100. For steady flows in Figure 4a-c, we see smooth and symmetric streamlines and isotherm contours. On the other hand, when the flow become unsteady in Figure 4d, these lines are distorted and less organized. Both for the steady and unsteady cases, thermal features discussed in Section 3 are noticed as well, including the relatively uniform temperature distributions inside the solid domains (due to the high solid conductivity K s ) and the smooth temperature transition across the solid-fluid interface (due to the large interface conductance C). Moreover, we see significant temperature variations along the unit boundaries for the unsteady system in Figure 4d. The variation fashions are similar on opposite edges. More specific, the temperature variations along the top and bottom boundaries appear identical, whereas the temperature is low on the left and high on the right. All these observations are consistent with the thermal boundary conditions in Equations (8) and (9).
For unsteady flow situations, the mean flow velocity U 0 and Reynolds number Re both vary with time, and it is difficult to achieve the time-averaged Reynolds number exactly at 100. For such situations, we accept the simulation results when the Re variation is limited in the range of 95 ∼ 105. The unsteady flow certainly affects the thermal field, and accordingly, the tortuosity and dispersion conductivities K tor and K dis calculated from Equations (14) and (15) also vary with time. Figure 5 displays the temporal oscillations of Re, K tor and K dis for Re = 100 and β = 1.44. The oscillation periods for K tor and K dis appear to be identical, as twice that for Re. In our next analysis, we use the average values over a variation period for the tortuosity and dispersion conductivities K tor and K dis . Figure 6 plots the apparent conductivities K tor and K dis changing with the aspect ratio β and Reynolds number Re. Clearly, even under the same porosity = 0.85, the macroscopic conductivities can be greatly affected by the aspect ratio β, and the influence can be enhanced with a higher Reynolds number Re. The tortuosity conductivity K tor decreases with β in Figure 6a, which can be explained by looking at the K tor definition Equation (14) and the temperature distributions in Figure 4. Equation (14) can be interpreted as a weighted sum of the fluid temperature over the solid-fluid interface Γ, with n x , the x-component of the normal vector n, as the weighting factor. For the system considered here, the left surface has n x > 0 and the right surface has n x < 0. On each semicircular surface, the portion around the horizontal centerline (the surface is approximately aligned in the vertical direction and thus it has a larger |n x | value) plays a more determinant role than the portions near the unit edges (the surface is approximately aligned in the horizontal and thus |n x | ≈ 0). Therefore, by comparing the temperatures on the parts around the horizontal centerline of the two semicircular surfaces in Figure 4, we can have a qualitative understanding of the K tor dependence on β: For a large β, the two semicircular surfaces are closer and the temperature difference on the closest portions becomes smaller, and this smaller temperature difference results in a smaller tortuosity conductivity, as observed in Figure 6a. The Reynolds number Re appears to be less influential on K tor , except for the highly unsteady case with β = 4. The weak influence of Re on K tor is consistent with that observed in Ref. [21]. Unlike the tortuosity conductivity K tor , the dispersion conductivity K dis increases by orders with β in Figure 6b. Similarly, we attempt to interpret this observation by examining the K dis definition Equation (15) and flow and thermal fields in Figure 4. Clearly, Equation (15) represents the disorder level of the flow and thermal distributions in the periodic unit, since the integrand is simply the product of the temperature fluctuation T −T and the flow fluctuation (u −ū f ) x . Therefore, the more disordered the flow and temperature distributions in the periodic unit, the larger K dis value will be obtained. This analysis is consistent with our general understanding of the thermal dispersion process and it explains the increasing trend of K dis with β and Re in Figure 6b as well. The similarity in flow and thermal fields between Re = 50 and 100 for β ≤ 1 steady systems suggests that K dis does not change much with Re there; however, for β > 1, the Re = 100 systems become unsteady and the flow and temperature distributions become significantly disordered, resulting in an abrupt increase in K dis . We also notice that only the x-component of the flow fluctuation is involved in Equation (15). For the steady flows at β ≤ 1, the nonzero x fluid velocity is limited in the small circulation areas near the four unit corners, whereas for large β values, the circulation region is large. The high Reynolds number Re = 100 is also helpful in increasing the circulation size, and it can even make the x velocity component nonzero over almost the entire fluid domain for the unsteady flows with β > 1. This is the reason that we see the K dis increase with β and Re by orders for β > 1 in Figure 6b.   Figure 6. Plots of K tor (a) and K dis (b) changing with the aspect ratio β at two Reynolds numbers Re = 50 (black squares) and 100 (blue circles).

Effect of the Interfacial Conductance C
Next we shift our attention to the influence of the interfacial thermal conductance C on the apparent conductivities. Again this is an important topic [32,33], however, it has not been investigated adequately. In this part, we fix the porosity = 0.85 and the Reynolds number Re = 50 for simplicity, and test three interfacial conductance values C = 10 (virtually no thermal resistance at the solid-fluid interface as observed in Figures 2b and 4), 5 × 10 −4 , and 5 × 10 −5 . Figure 7 displays the flow and thermal distributions at these three conductance values with two aspect ratios β = 0.44 and 1.44; and Figure 8 collects the calculated K tor and K dis conductivities in our simulations. The analysis in the previous section on the relationships between macroscopic conductivities and microscopic flow and thermal situations can still be applied here. Please note the change in interfacial conductance C does not affect the flow field and thus we can focus on the temperature response to different C values in Figure 7. As conductance C decreases, the solid domains become more insulated, resulting in a nearly constant temperature in each solid patch (i.e., the same colors and no isotherm contours in solid domains). With the solid part being insulated and thus its high-conductivity influence reduced, the fluid temperature around the solid domains exhibits a faster change along the thermal gradient direction. For example, the fluid temperatures near the solid surfaces at the narrowest gap location in Figure 7a are −0.0471 (left) vs. 0.0473 (right) for C = 10, −0.0364 vs. 0.0366 for C = 5 × 10 −4 , and −0.0251 vs. 0.0254 for C = 5 × 10 −5 . As discussed in the previous section, the large temperature difference at C = 10 introduces a larger tortuosity conductivity K tor , and vise versa, as shown in Figure 8a. On the other hand, the dispersion conductivity K dis increases as we reduce the interfacial conductance C, and the C influence becomes negligible for large aspect ratios. This change should be attributed to the temperature change in circulation areas; however, due to the complexity in flow structures and temperature changes with C in Figure 7, we are not able to provide a direct qualitative explanation here. Overall, the influence of the interface conductance C on the apparent conductivities is less dramatic in Figure 8 than that from the unit aspect ratio β in Figure 6, but it cannot be neglected for accurate analysis.

Effect of the Porosity
Unlike the effects of aspect ratio β and the interfacial conductance C, extensive studies have been conducted for the influence of the porosity on the macroscopic thermal behaviors of porous materials [21,61,62]. As mentioned in Section 1, previous pore-scale simulations usually set isothermal conditions on a pair of unit boundaries to create the thermal gradient and applied adiabatic conditions on other unit side surfaces [24,27]. Such artificial, unphysical treatments interfere with the natural heat transfer processes among microscopic periodic units, and the results from such boundary methods could be inaccurate or misleading [28]. In this section, we use the generalized periodic condition Equation (7) to examine the porosity effect on macroscopic conductivities K tor and K dis . Our next simulations cover a range of porosity of = 0.5∼0.99 for aspect ratios β = 0.69 and 1, and = 0.65∼0.99 for and β = 1.44. The Reynolds number is set at Re = 50 and the interfacial conductance is fixed at C = 10. Further reducing will increase the solid cylinder diameter D and thus make the gap between the solid surfaces very narrow ( Figure 9). As discussed in Section 4.1, such a narrow gap may turn the flow unsteady for the large aspect ratio β = 1.44, and also a finer spacial resolution is required to accurately capture the flow and thermal variations in the gap regions.  Results from this set of simulations are collected in Figure 9 for the flow and temperature distributions and Figure 10 for the calculated conductivities K tor and K dis at different porosity values. Please note that these figures are plotted in terms of the solid fraction 1 − as in Ref. [21]. With the solid diameter D increase (from left to right in Figure 9), a pair of circulation vortices are developed in the wake region above the solid cylinder, and the vortex size grows gradually till it completely fills the vertical space between two cylinders. Meanwhile, the original relatively organized temperature field is gradually distorted by the increasing size of the solid domain. Such changes in the flow and thermal patterns therefore affect the macroscopic thermal behaviors as characterized by the tortuosity (K tor ) and dispersion (K dis ) conductivities ( Figure 10). The continuous increase in K tor with the solid fraction 1 − is mainly because of the larger area (length in our 2D system) of the solid-fluid interface Γ (see Equation (14)). This trend is similar to that observed in Ref. [21]; however, the K tor magnitude is smaller here. This can be attributed to the different solid domain shapes: For the same porosity, compared to the circular shape in this work, the square solid domain shape in Ref. [21] yields a longer interface length, and more profoundly, half of the interface has the local normal direction aligned the x direction. All these are favorable for a larger K tor according to its definition in Equation (14). As for the dispersion conductivity K dis in Figure 10b, in general, K dis increases with the solid fraction 1 − ; however, unlike the monotonic growth in Ref. [21], local maximum and minimum states are observed here. Due to the complexity in flow and temperature fields as well as their involvement in the K dis calculation, it is difficult to provide detailed insights and mechanism for the K dis variations. Nevertheless, this could be an interesting topic to explore in the future. The general increasing trend can be qualitatively related to the increasing size of the circulation region, which is the major contributor to K dis via the large x velocity component in this region. A similar analysis can be applied for the gentle variation and slow recovery of K dis for 1 − = 0.2∼0.5 at β = 0.69 in Figure 10b, since the flow pattern remains almost unchanged in these systems (see the subplots for 1 − = 0.2 and 0.35 at β = 0.69 in Figure 9). Finally, for the three curves with different aspect ratios, we see K tor is smaller but K dis is larger for a higher aspect ratio β. This agrees well with our findings and discussions for the aspect ratio effect on conductivity in Section 4.1 (see Figure 6).

Summary
In this study, we have first justified and validated the generalized periodic condition [21,28] for pore-scale simulations of thermal flows in porous media. Extensive simulations have then been carried out using a 2D porous model to investigate the influence of several microscopic parameters on macroscopic thermal conductivity. Among them, the effects of the aspect ratio of the periodic unit and the thermal conductance at the pore surface have not been addressed adequately in previous studies. Our results show that these microscopic properties can dramatically change the flow and thermal fields in the microscopic porous structure, and affect the apparent thermal performances of the porous materials at the macroscopic level. Therefore, these microscopic factors need to be considered carefully for more accurate and reliable simulation results, which are crucial for both fundamental research and practical applications. In addition, thorough discussions are attempted to qualitatively explore the relationship between the apparent conductivity at the macroscopic level and the complex thermal flow situations in the microscopic porous structure, and our analysis and findings could be helpful for a better understanding of the underlying thermal processes.
We are aware that several serious limitations exist for this study and one should not over interpret the results obtained from a specific system. The simple 2D porous model and the relatively low Reynolds numbers considered here may appear less realistic for practical systems; however, they are helpful for understanding the micro-macro relations and fundamental mechanisms involved in the complex thermal flow processes in porous materials. We have fixed the solid and fluid substance properties (conductivity, diffusivity, and the Prandtl number), which can certainly affect the apparent thermal behaviors as well. The 2D unit geometry adopted in this work is also symmetric along both the flow and thermal gradient directions; and the anisotropic effects could be an interesting topic for future research. The periodic conditions Equations (8) and (9) utilized in these calculations require the flow and thermal fields to be fully developed and thus they are not applicable to the entrance and development regions [63,64]. Moreover, in our simulations, we have assumed that the material properties are constant in the flow and heat transfer processes and thus steady or quasi-steady states can be established. In some situations, the microscopic porous structure may undergo a dynamic change due to swelling and erosion [65]; and caution must be taken for applying results from this study to such systems. Nevertheless, the boundary method, simulation results, and analysis discussions can be useful for the research and applications of porous thermal flows. The generalized periodic boundary condition, although presented in 2D, can be readily applied to three-dimensional pore-scale thermal flow simulations.