CO2 methanation: on the modeling of reacting laminar flows in structured Ni/MgAl2O4 catalysts

Journal of Physics: Conference Series PAPER • OPEN ACCESS CO2 methanation: on the modeling of reacting laminar flows in structured Ni/MgAl2O4 catalysts To cite this article: J A Medina Méndez et al 2022 J. Phys.: Conf. Ser. 2367 012015 View the article online for updates and enhancements. You may also like Absolute OH density and gas temperature measurements by laser induced fluorescence in a microsecond pulsed discharge generated in a conductive NaCl solution Jianan Wang, Marien Simeni Simeni, Mingzhe Rong et al. - Generation and dissipation of runaway electrons in ASDEX Upgrade experiments G. Pautasso, M. Dibon, M. Dunne et al. - Modeling of Ion Crossover in Vanadium Redox Flow Batteries: A ComputationallyEfficient Lumped Parameter Approach for Extended Cycling Philipp A. Boettcher, Ertan Agar, C. R. Dennison et al. - This content was downloaded from IP address 150.214.182.230 on 15/02/2023 at 15:06 Content from this work may be used under the terms of theCreativeCommonsAttribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Published under licence by IOP Publishing Ltd XXV Fluid Mechanics Conference (FMC 2022) Journal of Physics: Conference Series 2367 (2022) 012015 IOP Publishing doi:10.1088/1742-6596/2367/1/012015 1 CO2methanation: on the modeling of reacting laminar flows in structured Ni/MgAl2O4catalysts J A Medina M´endez1, M Gonz´alez-Casta˜no2, F M Baena-Moreno2,3 and H Arellano-Garc´ıa2 1Numerical Fluid and Gas Dynamics, BTU Cottbus-Senftenberg, Siemens-Halske-Ring 15a, 03046 Cottbus, Germany 2Process and Plant Technology, BTU Cottbus-Senftenberg, Burger Chaussee 2 LG 4/1, 03044 Cottbus, Germany 3Chemical and Environmental Engineering Department, Technical School of Engineering, University of Seville, 41092 Sevilla, Spain E-mail: [email protected] Abstract. This work presents an overview of issues for the modeling of laminar flows in monolith catalysts. Both 0-D and 3-D models are evaluated for a parallel channel structured honeycomb catalyst (PC-HC), and a gyroid 3-D printed structured catalyst (G-3D). At the 0-D homogeneous reactor modeling level, the analysis is focused on the effect of the bulk porosity, as well as the model choice to represent Nusselt number effects. Results show the better suitability of a long tube Nusselt number model for the representation of the maximum temperature achievable in the 0-D homogeneous reactor, as well as a modest effect of the porosity on the catalyst CO2conversion. A more detailed insight on heat transfer and the core reaction zone inside the monolith can be obtained at the 3-D homogeneous reactor modeling level. 1. Introduction Fluid dynamic treatment of structured catalysts is generally conceived within the framework of porous media flows. For homogeneous media, the most widely used closure model for the upscaling of the momentum balance using volume-averaging theory (VAT) is the DarcyForchheimer equation [1]. In the experimental field, the use of VAT generally leads to confusion regarding the nature of the flow. Specifically, the use of a variety of Reynolds numbers for flow characterization, as well as deviations from the classical Darcy’s law in the observed regimes, have led to some confusion regarding the laminar or turbulent nature of the flow, see [1]. In this work, we focus on the macroscopic (volume-averaged) evaluation of the flow over the porous media, and its corresponding bulk flow Reynolds number Re. Within this volume-averaged point of view, structured catalysts generally operate at laminar flow regimes due to the very small cross-sectional areas of their associated monolith channels. It is clear that, on one hand, a structured catalyst is better represented in terms of a heterogeneous model, which is faithful to the dynamics of both the solid and the working fluid phases. Homogeneous reactor models, on the other hand, are still widely used for industrial applications due to its simplicity and highly efficient use of computational resources. Despite the intrinsic loss of accuracy associated to homogeneous reactor models, the use of heterogeneous reactor models may not always be justified [2]. XXV Fluid Mechanics Conference (FMC 2022) Journal of Physics: Conference Series 2367 (2022) 012015 IOP Publishing doi:10.1088/1742-6596/2367/1/012015 2 Within the widespread field of applications for structured catalysts, one application with current strategic importance is that of emissions control. The CO or CO2methanation has some fundamental significance, since it is associated both to natural gas production, as well as CO and CO2emissions removal [3]. For the optimization of the process of natural gas production, the development of selective, reliable and efficient catalysts is required. Although structured catalysts have been received with initial skepticism due to their poor performance in comparison to pellet beds, as well as their virtual adiabatic character, specific requirements for structured catalysts can be met with dedicated designs [4]. One example of such dedicated design is the gyroid 3-D printed catalyst presented by [5]. The latter uses Ruthenium-promoted magnesiumaluminate supported Nickel as a catalyst material, 0.5%Ru-15%Ni/MgAl2O4. Preliminary evaluations of this catalyst have already demonstrated a reliable and highly efficient CO2 conversion. However, the optimization of the design is complicated due to the multitude of factors playing a role in the catalytic efficiency: from the chemical kinetics taking place, all the way to the catalyst-enhanced flow drag, and the associated heat and mass transfer. We focus on the evaluation of the recently introduced Ru-Ni/MgAl2O4structured catalysts [5], from the numerical fluid-dynamic point of view. This paper is structured as follows. In section 2, we briefly detail the relevant aspects associated with a zero dimensional homogeneous model for the numerical fluid-dynamic evaluation of structured catalysts. Next, in section 3, we explain some aspects concerning the three dimensional homogeneous modeling. The results for the evaluation of a parallel channel structured honeycomb catalyst (PC-HC), and a gyroid 3-D printed structured catalyst (G-3D), as in [5], are shown in section 4. Finally, closing remarks and an outlook of potential future work is shown in section 5. 2. Zero-dimensional structured catalyst model An intended 0-D homogeneous model for structured catalysts follows from the model developed by [6]. The velocity Ubis assumed constant and homogeneous in the 0-D reactor. The model relates streamwise advancement of plug-flow-like reactor models to time, by the spatial-totemporal transformation dt= dz/Ub, where dzis a differential of the reactor streamwise coordinate, and dtis a differential of time. The reactor time-scale tτis then obtained as L/Ub, where Lis the monolith length. The chemical species (k-th mass fractions) Ykand temperature Tequations are formulated for one reactor volume, in the form of a balance, and are then expressed per unit volume, in order to obtain the corresponding ordinary differential equations (ODEs) for the homogeneous fluid phase, ρdYk dt= ˙ωkWk(1) ρcp dT dt=X jhρcat ˙ωR,j −∆HR,0i−4 DUA(T−Tw) (2) Here, is the bulk porosity of the catalyst, ˙ωkis a molar-based species reaction rate and Wk is the k-th species molecular weight. Also, ρis the fluid density, ρcat is the catalyst density (catalyst mass per reactor volume), and cpthe specific heat capacity at constant pressure. ˙ωR,j is the j-th reaction net molar rate of progress per unit mass of catalyst, i.e., in kmol/(kgcat·s), specified by the used chemical kinetics mechanism. This is related to the species reaction rate by ˙ωk=ρcat Pjνjk ˙ωR,j, being νjk the stoichiometric coefficient of species kin reaction j. Also, ∆HR,0is the reaction enthalpy at reference temperature conditions, e.g., 298 K. Furthermore, Dis the monolith diameter, UAis a global convective heat transfer coefficient, and Twis a constant reactor cooling wall temperature. Equations (1) and (2) can be integrated in time using an explicit Euler method. XXV Fluid Mechanics Conference (FMC 2022) Journal of Physics: Conference Series 2367 (2022) 012015 IOP Publishing doi:10.1088/1742-6596/2367/1/012015 3 2.1. Chemical kinetics Chemical kinetics are modeled in this work by the methanation reaction mechanism developed by [7]. This is a reduced chemistry mechanism for Ni/MgAl2O4-spinel catalysts, with 3 reactions and 5 chemical species (H2, H2O, CH4, CO, CO2). N2is added as an inert specie. This mechanism was also used by [6] for the study of CO2methanation reactors. It is noted that this mechanism may not incorporate the effects of the added Ru in the intended catalyst to evaluate, i.e., the Ru-Ni/MgAl2O4catalyst from [5]. Both the reaction mechanism and the thermodynamic data for the species contained in the mechanism are supplied in the form of a CTI file to the open source library Cantera [8], in order to handle the reaction rate calculations. The methanation mechanism developed by [7] also contains the presence of so-called DEN inhibition terms. These limit the reaction rates described by the rate coefficients supplied in [7]. The inhibition terms are calculated in this work in an explicit way, in addition to the Cantera-calculated reaction rates, and are applied as multiplicative factors to such rates afterwards. 2.2. Convective heat transfer The global heat transfer coefficient UAis calculated in [6] based on a Nusselt number (Nu) correlation for the heat conduction at the wall, and a radial effective transport specialized for fixed-bed reactors. More generally, UAis related to the bulk Nusselt number as UA=Nu(t)λ Ω(3) Here, λis the thermal conductivity of the fluid, Ω is a characteristic length, and Nu(t) is a time-dependent Nusselt number. We consider the analytically derived solutions to the Nusselt problem obtained by [9] and [10]. In [9], asymptotic analytical solutions are presented for short reactors with uniform inlet flow, as well as fully developed laminar parabolic inlet profiles. Short reactors are defined at negligible aspect ratio, or conversely, at very large transverse Peclet number PeT>> 1, where PeT=UbΩ02/(LDm), being Ω0=Dh/2 the hydraulic radius in the short reactor case, and Dm a mass-diffusion coefficient. For PeT>> 1, the Nusselt scaling considering flat uniform inlet velocity profiles in a (constant pressure) plug-flow reactor (PFR), is then given by [9], previous consideration of the space-to-time analogy discussed in this work, Nu(t)=2sPeT Lef (4) Note that for this case, Ω = Dh, the hydraulic diameter of the monolith, is the correct characteristic length associated to Nu or UA, see [9]. Also, Lefis the fluid Lewis number in equation (4), where Lef=λ/(ρcpDm). The time-dependency for Nu in equation (4) arises from the temperature dependency of the fluid properties. We also evaluate a long tube Nu model. In [10], a solution to the Nu problem for long reactors is obtained (PeT<< 1). In the case of long reactors, Ω0=Dh/4, while Ω = Ω0is the correct characteristic length to associate with Nu, see [10]. We focus again on the flat uniform inlet velocity profile, for a constant operating pressure and fixed wall temperature reactor, i.e., a developing flow consideration for which Nu has now a streamwise dependency Nu = Nu(z), see [10], Nu(z) = Nu∞+0.44Pr−1 6PeTL z 1+0.314qPeTL z (5) The streamwise-dependency can be exchanged by the time-dependency in our model. In equation (5), Pr = cpµ/λ is the Prandtl number, with µbeing the fluid dynamic viscosity. XXV Fluid Mechanics Conference (FMC 2022) Journal of Physics: Conference Series 2367 (2022) 012015 IOP Publishing doi:10.1088/1742-6596/2367/1/012015 4 Also, Nu∞= 2.966 is the asymptotic Nusselt number value for a honeycomb-shaped channel cavity in [10]. The hydraulic diameter, Dh, is calculated as a function of the monolith total volume, and the interfacial fluid and catalyst area, which corresponds to the volume analog of the traditional definition based on the surface area and the wet perimeter. Dhcan then be estimated as Dh=4 4 D+ (1 −)4S πD2L (6) where Sis the supplied catalyst surface area. It is noted that the estimation of Dhalso allows the estimation of the number of channels Nchan in the monolith, assuming that there is one channel per unit area by Dh, that is, Nchan =D2/D2 h. Figure 1 shows results for a PFR evaluation case using reference data from [6]. The PFR model used by [6] is compared to the short and long 0-D reactor models explained before. Figure 1 shows that the model of [6] is more dynamically similar to a short reactor Nu model. In general, the short tube Nu model minimizes the convective heat transfer, thus allowing larger maximum process temperatures and reaction product yields for the same initial or inlet conditions. That is, the short tube Nu model imposes the dominance of fast reactions, the faster ignition-like behavior in figure 1. Figure 1. Maximum temperature in the PFR for a simulated reactor timescale tτ, obtained as a function of the feed (initial or inlet) temperature Tfeed; see input parameters in [6]. The pseudohomogeneous reactor model data from [6] is shown with a continuous line (——) and markers (•). The short tube Nu model described in this work is shown with a dashed line (– – –), while the long tube Nu model is shown with a dotted line (······). 3. Three-dimensional structured catalyst model We now detail some considerations for a 3-D homogeneous structured catalyst model. For the 3-D model, we have resorted to the use of the open source library OpenFoam (v8) [11]. Specifically, we resort to the use of the reactingFoam solver, which solves a system of equations for the conservation of chemical species (mass fractions), mass (in the form of density), momentum, and total energy (formulated in terms of enthalpy). Assuming homogeneous media, i.e., uniform porosity, these equations are adapted from the homogeneous model in [2], that is, ∂(ρYk) ∂t +∇·(ρV Yk) = −∇·(ρV kYk) + ˙ωkWk(7) ∂ρ ∂t +∇·(ρV ) = 0 (8) ∂(ρV ) ∂t +∇·(ρV ◦V) = −∇P+∇·τ+S(9) ∂(ρfsh) ∂t +∇·(ρV h) + ∂(ρeK) ∂t +∇·(ρV eK)−∂P ∂t =∇·qeff −ρcat X j ˙ωR,j∆HR,0(10) XXV Fluid Mechanics Conference (FMC 2022) Journal of Physics: Conference Series 2367 (2022) 012015 IOP Publishing doi:10.1088/1742-6596/2367/1/012015 5 Here, ∇is the traditional Nabla operator, with ·and ◦symbolizing the scalar and dyadic products. Additionally, Vis the velocity field, while Vkis a Fickian diffusion velocity field for species k. Likewise, Pis the thermodynamic pressure, Tis the temperature, τis the shear stress tensor, eK=|V|2/2 is the specific kinetic energy, and qeff =λeff ∇Tis an effective heat flux, where λeff is an effective thermal conductivity. This is assumed to have a uniform value in axial direction λeff,z, and a different uniform value in radial direction λeff,r. The values λeff,z and λeff,r considered for this work are those mentioned in [5]; see [4] for related details. Also, ρfs is a fluid-solid weighted density, similar to the treatment in [2], i.e., ρfs = [ρcp+ (1 −)ρscp,s]/cp, where ρsand cp,sare the density and specific heat capacities of the solid matrix. The enforcement of the zero-net species mass-diffusion flux is achieved by the offset of errors into inert N2[12], i.e., YN2= 1 −Pk6=N2Yk. The chemical kinetics mechanism used for the 3-D homogeneous reactor model is the same one considered for the 0-D homogeneous reactor model. Unlike in the 0-D model, Cantera is not used for the processing of the reaction rates. Instead, OpenFoam itself is used. It is noted that the calculation of the equilibrium constant in OpenFoam differs from the calculation method used in Cantera, thus yielding different reaction rates (see the online OpenFoam v8 documentation). The solution of equations (7-10) considers uniform initial conditions. The initial composition considered is that of a pure H2flow at the specified reactor feed temperature. This is similar to the reactor heating situation for catalyst pre-treatment in [7]. Likewise, zero Neumann boundary conditions (BCs) are considered for the chemical species at the reactor walls and reactor outlet, while fixed Dirichlet BCs are considered for the inlet gas composition. The velocity field considers no-slip BCs at the reactor walls, and inlet-outlet BCs. Finally, the temperature field considers zero Neumann BCs at the reactor outlet, fixed inlet gas temperature Dirichlet BCs, and fixed reactor (cooling) wall Dirichlet BCs, Tw. 3.1. Porous media modeling The modeling of the porous media enters the momentum equation (9) via the additional drag term S. Following the assumptions for homogeneous porous media, Sis modeled following the Darcy-Forchheimer equation [1]. In practical terms, this is, S=−µ Kuˆez−F √Kρ|u|uˆez(11) Here, uis the streamwise component of the velocity field, whereas ˆezis the canonical streamwise vector, and µis the dynamic viscosity. The permeability of the media Kis estimated as in [13], using the hydraulic radius model and the Carman-Kozeny equation, K=3d2 36kK(1 −)2(12) where dis an equivalent spherical diameter obtained for the porous media, according to the specific surface area of the catalyst based on the solid volume, see [13]. Furthermore, kKis the Kozeny constant, approximated as 5. Finally, Fin equation (11) is obtained using the modified Ergun relation [13], such that F=1.8 √1805(13) 4. Case configuration and results Dimensional parameters corresponding to the PC-HC and the G-3D structured catalyst reactors introduced in [5] are shown in table 1. In tables 2 and 3, we detail the operating and feed conditions for both the 0-D homogeneous reactor simulations, and the 3-D homogeneous reactor simulations. These are also taken from [5]. XXV Fluid Mechanics Conference (FMC 2022) Journal of Physics: Conference Series 2367 (2022) 012015 IOP Publishing doi:10.1088/1742-6596/2367/1/012015 6 Table 1. Dimensional parameters for PC-HC and G-3D structured catalyst reactors. Parameter Dimension Value PC-HC Value G-3D Monolith diameter (D) cm 1.6 1.6 Monolith or reactor length (L) cm 3.0 3.0 Catalyst surface area (S) cm2540.0 504.0 Porosity () - 0.88 0.79 Estimated hydraulic diameter (Dh) mm 2.66 1.58 Effective radial thermal conductivity (λeff,r) W/(m K) 1.23 2.17 Effective axial thermal conductivity (λeff,z) W/(m K) 2.31 3.91 Solid matrix specific heat capacity (cp,s) J/(kg K) 466.0 466.0 Solid matrix density (ρs) kg/m37850.0 7850.0 Table 2. Operating conditions evaluated in the simulations. Case Condition Dimension Value All cases Weight hourly space velocity (WHSV) m3/(kgcat s) 0.0167 Gas volume flow rate (Q) m3/s 3.3333 ×10−6 Catalyst mass (mcat =Q/WHSV) kg 2.0×10−4 Catalyst mass per reactor volume (ρcat) kg/m31989.44 Reactor (constant) pressure kPa 101.325 Table 3. Feed compositions evaluated in the simulations. Case Condition Case I (0-D) Feed stream molar composition: 15% CO2+ 60% H2+ 25% N2 Inlet and reactor wall temperature: 473.15 −673.15 K Case II (0-D) Feed stream molar composition: 8% CO2+ 32% H2+ 12% CH4+ 48% N2 Inlet and reactor wall temperature: 473.15 −673.15 K Case III (3-D) Feed stream molar composition: 15% CO2+ 60% H2+ 25% N2 Inlet and reactor wall temperature: 573.15 K 4.1. Zero-dimensional homogeneous reactor The newly developed PC-HC and G-3D structured catalysts from [5], are evaluated in this section using a 0-D homogeneous reactor long tube Nu model. This evaluation corresponds to the numerical modeling of Cases I and II in table 3, using the operating conditions from table 2. In figure 2(a), the experimental results for CO2conversion for a Case I feed composition are shown and compared to the long 0-D reactor numerical model detailed in this work. The CO2 conversion is defined based on the mole fractions of CO2,XCO2, at the inlet and outlet of the XXV Fluid Mechanics Conference (FMC 2022) Journal of Physics: Conference Series 2367 (2022) 012015 IOP Publishing doi:10.1088/1742-6596/2367/1/012015 7 reactor (or start or end of the simulated reactor time-scale), xCO2=XCO2,inlet or start −XCO2,outlet or end XCO2,inlet or start (14) The numerical results observe a large deviation with respect to the experimental results. First of all, we note that the chemical kinetics developed by [7], which are used in the numerical model, yield a lower CO2conversion in comparison to the actual experimental measurements. The kinetics developed by [7] correspond to Ni/MgAl2O4catalysts, while the catalyst used in [5] is doped with Ru, i.e., Ru-Ni/MgAl2O4. Additionally, it is possible that by neglecting intermediate chemical species in the simplified chemistry, the established chemical equilibrium differs from the real reactor equilibrium. We observe that the real reactor runaway to chemical equilibrium is seen to happen at much higher temperatures in comparison to the numerical 0-D model, even when using the long tube Nu model. Despite the differences between experiments and numerics, the superior performance of the G-3D catalyst is confirmed in this work. The superior performance is attributed to the larger catalyst surface area of the G-3D catalyst, as suggested in [5]. Figure 2(b) shows similar CO2conversion results for a simplified biogas feed composition (Case II). Figure 2. CO2conversion as a function of Tfeed. (a) Feed composition Case I. (b) Feed composition Case II. Refer to tables 1-3 for operating conditions and feed compositions. Experimental data from [5] is shown with dashed lines (– – –) and markers (•), while simulation data is shown with dotted-dashed lines (— ·—). Orange and magenta lines refer to the G-3D and PC-HC catalysts, respectively. The experimental chemical equilibrium, and the equilibrium obtained by the kinetics of [7], are shown with gray lines. 4.2. Three-dimensional homogeneous reactor As an additional numerical tool for the evaluation and analysis of the experimental results in [5], we now perform 3-D homogeneous reactor simulations corresponding to the feed composition Case III (see table 3). Figure 3 shows the inlet, outlet and longitudinal fields for the CO2mass fraction in the simulated PC-HC catalyst, at half of the reactor time-scale. Results for the G-3D catalyst are very similar. There are no appreciable differences between the cases, and both yield CO2conversion values of approximately 90%, based on the cross-sectional averaged XCO2of the reactor inlet and outlet, see equation (14). The ∼90% value is of course larger than that obtained in the 0-D model for the same simulated inlet temperature (573.15 K). The reason for this are the higher flow temperatures in the regions away from the wall and the inlet, according to the temperature and composition BCs, in the 3-D model. XXV Fluid Mechanics Conference (FMC 2022) Journal of Physics: Conference Series 2367 (2022) 012015 IOP Publishing doi:10.1088/1742-6596/2367/1/012015 8 Figure 3. Snapshots of CO2mass fraction field at half of the reactor time-scale (feed composition Case III). (a) Inlet cross-section. (b) Outlet cross-section. (c) Longitudinal section. 5. Discussion and conclusions We presented a framework for the homogeneous modeling of 0-D and 3-D structured catalysts. The modeling framework was applied for the numerical evaluation of the novel structured RuNi/MgAl2O4catalysts introduced in [5]. For that, chemical kinetics developed by [7] were used. Despite the numerical modeling being short on reproducing the experimental data, the numerical results confirm the superior performance of the G-3D catalyst in comparison to the PC-HC. This is solely due to the different porosity and catalyst surface area. The 3-D simulations were not conclusive in this work due to technical difficulties in the handling of the reaction rate calculation by OpenFoam, specifically regarding the calculation of the equilibrium constants, which differs from the method used in Cantera. In future work, we will use the Cantera library in OpenFoam for the calculation of the reaction rates, as a way to verify the 3-D results. The ultimate goal is to have a consistent comparison with the 0-D model, which simultaneously allows the adjustment of the chemical kinetics, and the evaluation of the monolith fluid-dynamic performance based on the porosity, drag, heat transfer, CO2conversion, among other variables. In this way, a reliable numerical tool can be developed for the future optimization of the catalyst design. References [1] Wood B D, He X and Apte S V 2020 Annu. Rev. Fluid Mech. 52 171–203 [2] Liu B, Hayes R E, Yi Y, Mmbaga J, Checkel M D and Zheng M 2007 Comput. Chem. Eng. 31 292–306 [3] Gao F and Goodman D W 2012 Annu. Rev. Phys. Chem. 63 265–286 [4] Tronconi E, Groppi G and Visconti C G 2014 Curr. Opin. Chem. Eng. 555–67 [5] Baena-Moreno F M, Gonz´alez-Casta˜no M, Navarro de Miguel J C, Miah K U M, Ossenbrink R, Odriozola J A and Arellano-Garc´ıa H 2021 ACS Sustainable Chem. Eng. 98198–8206 [6] Schlereth D and Hinrichsen O 2014 Chem. Eng. Res. Des. 92 702–712 [7] Xu J and Froment G F 1989 AIChE J. 35 88–96 [8] Goodwin D G, Speth R L, Moffat H K and Weber B W 2018 Cantera: An objectoriented software toolkit for chemical kinetics, thermodynamics, and transport processes https://www.cantera.org version 2.4.0 [9] Gupta N and Balakotaiah V 2001 Chem. Eng. Sci. 56 4771–4786 [10] Gundlapally S R and Balakotaiah V 2011 Chem. Eng. Sci. 66 1879–1892 [11] Weller H G, Tabor G, Jasak H and Fureby C 1998 Comput. Phys. 12 620–631 [12] Poinsot T and Veynante D 2014 Theoretical and Numerical Combustion (Aquaprint) [13] Kaviany M 2012 Principles of heat transfer in porous media (Springer Science & Business Media)