Drag and heat transfer closures for realistic numerically generated random open-cell solid foams using an immersed boundary method

a r t i c l e i n f o


a b s t r a c t
In this paper, we apply a novel immersed boundary method to simulate pore-scale level fluid flow and convective heat transfer in realistic numerically generated open-cell solid foams in a Cartesian computational domain. Five different periodic foam samples of varying porosities (e ¼ ½0:877; 0:948) are generated by numerically mimicking the actual foam formation process (minimizing surface area). The step-by-step procedure for generating the periodic foam geometries is presented. The specific surface areas of the generated foams of different porosities are compared with real foam geometries showing a reasonable agreement. The Reynolds number (Re) is varied from Re % 0 (creeping flow) to Re % 500, and finally drag and Nusselt correlations have been proposed. A detailed analysis is presented on the local velocity and temperature field for the fluid-solid interaction in a complex cellular porous medium.

Introduction
Open-cell solid foam is a novel material consisting of interconnected ligaments forming a complex network of randomly oriented polyhedrons. Because of its excellent thermal, mechanical and acoustic properties, such foams are extensively used in several industrial applications like heat exchangers, thermal energy absorbers, vaporizers, heat shielding devices, etc. Cellular foams are generally made of metals (e.g. aluminum, nickel, iron-chromium, etc.) or ceramics (e.g. Al 2 O 3 , silicon, carbide, etc.), and the choice of foams depends on the application area. In fixed bed catalytic reactors open-cell solid foams are gaining popularity as catalyst support because of their large surface to volume ratio. It improves the gas-liquid contacting and thereby enhances the heat and/or mass-transfer rate with a minimal increase in pressure drop as compared to other conventional packings.
Due to very complex and random geometrical shapes, most of the reported work in the literature dealing with solid foams is experimental. A comprehensive review on several experimental studies is given by Dietrich (2012). An extensive review of the potential applications of the metal foams as a heat transfer material is given by Ackermann et al. (2014). They have reported a large deviation between the existing literature due to (a) morphological differences among the different commercially available solid foams, and (b) experimental inaccuracies. In contrast, fully resolved numerical simulations at pore-scale level are capable to overcome all experimental uncertainties and to precisely account for specific geometrical details. At the same time arbitrary material properties and flow conditions can be defined and perfectly controlled in numerical simulations. Considering an idealized geometry, many researchers performed pore-scale simulations for opencell foams by resolving the detailed solid structure in a 3D computational domain. Kelvin's unit cell (Tetrakaidecahedron) (Chung et al., 2006;Horneber et al., 2012;Lucci et al., 2014;Bai and Chung, 2011;Das et al., 2017b) and Weaire-Phelan cell (Boomsma et al., 2003) are the most popular geometrical shapes to represent the foam unit cell structure in a periodic domain. One step ahead of the idealized unit cell assumption, several authors (August et al., 2015;Bianchi et al., 2013Bianchi et al., , 2012Bianchi et al., , 2015aZafari et al., 2015;Ranut et al., 2014) presented accurate and assumption free numerical results for real foam structures where a 3D image captured from a micro-CT scan has been converted into a CFD readable geometry. However, those results are only valid for the very particular geometry, mostly presented in a dimensional form and hence those results loose their generality. Authors like August et al. (2015), Wehinger et al. (2016), Kopanidis et al. (2010), etc. tried to create random foam samples numerically to perform heat and fluid flow simulations. However, all of them have considered a small section of the generated foam in a channel with inlet and outlet boundary conditions and hence those results are affected by boundary effects. In Table 1 the numerical studies on real (reconstructed from l À CT images) and realistic (numerically generated) open-cell solid foam geometries have been listed.
The main motivation of the present work is twofold: (i) to generate realistic foam samples and (ii) to derive fluid drag and heat transfer closures for the open-cell solid foams using an immersed boundary method. The realistic random micro-structures of foam are generated with Surface Evolver (Brakke, 1992) simulations which inherently mimic the foam forming process (i.e. minimization of surface energy). The generated random foam samples are periodic in all the three directions that help us to perform simulation in a comparatively smaller computation domain without loosing any generality. To the best of our knowledge, it is for the first time that realistic random periodic structures have been used to perform flow and heat transfer simulations for developing fluid flow and heat transfer closures.
In the accompanying paper (Das et al., submitted for publication) we have presented an accurate immersed boundary method to perform pore-scale simulations for highly complex porous structures in a periodic domain. A very detailed numerical validation and verification have been also presented. The immersed boundary method is a Cartesian grid based computational technique where the complex non-body conforming solid surface is treated as an ''immersed boundary" to impose boundary conditions. The same in-house code is applied to fully resolve the foam geometry on Cartesian grids.
In our previous work (Das et al., 2017b), we have considered an idealized foam geometry (Kelvin's unit cell model with cylindrical strut morphology) in a periodic computational domain. Simulations were performed for a wide range of porosities (e ¼ ½0:638; 0:962]) and Reynolds numbers (Re ¼ ½0 500). Flow was enforced by applying a constant body force for three different A nondimensional drag/pressure drop correlation was proposed for the entire data set. The present work is a natural extension of our previous work. However, we have used a realistic foam geometry. Five different periodic foam samples of varying porosities (e ¼ ½0:877; 0:948) are used where the Reynolds number (Re) is varied from Re % 0 (creeping flow) to Re ¼ 500.

Numerical methodology
For the present fully resolved pore-scale level numerical simulations, the conservation equations for mass, momentum and energy are solved in a 3D Cartesian computational domain. The fol- Table 1 Summary of previous pore-scale level numerical studies for real (reconstructed from l À CT images) and realistic (numerically generated) open-cell solid foams.

Reference
Foam geometry Analysis/study Zafari et al. (2015) l À CT scan of actual foam samples; e ¼ 0:85-0:95 Hydrodynamics and convective heat transfer Della  l À CT scan of actual foam samples; e ¼ 0:79-0:96 Hydrodynamics Bianchi et al. (2013Bianchi et al. ( , 2012 l À CT scan of actual foam samples; e ¼ 0:89-0:95 Hydrodynamics and convective heat transfer Ranut et al. (2014) l À CT scan of actual foam samples; e ¼ 0:92-0:94 Hydrodynamics, convective heat transfer and effective (stagnant) thermal conductivity Iasiello et al. (2014) l À CT scan of actual foam samples; e ¼ 0:894 Hydrodynamics and convective heat transfer Wehinger et al. (2016) Numerically generated by filling voidage between random packings of bidisperse spheres; e ¼ 0:74-0:82 Heat and species transport Bodla et al. (2010) l À CT scan of actual foam samples; e ¼ 0:93 Hydrodynamics, convective heat transfer and effective stagnant thermal conductivity Kopanidis et al. (2010) Numerically generated by surface energy optimization using Surface Elolver (Brakke, 1992); e ¼ 0:97 Hydrodynamics and convective heat transfer August et al. (2015) Numerically generated by filling voidage between random packings of monodisperse spheres; e ¼ 0:85-0:92 Effective stagnant thermal conductivity lowing governing equations are solved for velocity ( u), pressure (p) and temperature (T) respectively: qC p @T @t þ qC p r Á ð u TÞ ¼ r Á ðk rTÞ ð 3Þ where q; l; C p and k are the density, viscosity, specific heat capacity and thermal conductivity of the fluid, respectively. An immersed boundary method is used to enforce boundary conditions at the fluid-solid interface. In our implementation complex solid structures can be imported as a triangulated surface mesh in the CFD domain where boundary conditions at the fluid-solid interface are treated implicitly. The method shows a second order spatial convergence where 20 cells across the diameter produce grid independent results for Reynolds number up to 500. More details of the present IBM implementation can be found in the accompanying paper (Das et al., submitted for publication). All the direct numerical simulations (DNS) are performed in periodic domains using periodic foam geometries. The advantages of using periodic conditions (fully-developed flow) are that there are no entrance length effects and simulations can be performed in a reasonably small domain. For temperature field, we have used Prandtl number (Pr ¼ lC p =k) equal to unity, which assures equal boundary layer thickness (and as a result same grid requirements) for the velocity and temperature field. The temperature gradient inside the solid is neglected and hence a constant wall temperature is applied at the fluid-solid interface. In the periodic domain the flow is sustained by a source term in the momentum equation. For moderate to high Peclet number (Pe) flows (negligible axial conduction), with the solid surface is at constant temperature, T S , the non-dimensional temperature field (H) repeats its exact morphology in the stream-wise flow direction (here we consider along the x-direction), i.e. Hðx; y; zÞ ¼ Hðx þ L; y; zÞ; where H ¼ In the above equation T b j x represents bulk or cup-average temperature at a constant x-plane. At low Peclet number, due to finite axial conduction we cannot use the scaled-periodic temperature boundary condition (Das et al., submitted for publication). For such cases a Danckwerts' boundary condition (Danckwerts, 1953) is applied across the stream-wise flow direction where the velocity field is also periodic. It assumes that at low Pe, for thermally fully developed flows, where the reference temperature, T Ref can be considered as a developing far stream temperature field. It has no influence on the nondimensional quantities and in the present case we have used T Ref ¼ 0.

Reconstruction of open-cell solid foam samples
Open-cell solid foams are generally manufactured in a casting process where in a molten metal pool high pressure gas bubbles are induced. It finally creates the open structures after solidification. Another very common way is to use molds made by random packings of polymer spheres and in the casting process the interstitial space between the spheres is filled with molten material. After solidification of the metal the spheres are burned and removed from the system. The final geometrical shape of the foam depends on several factors, however at the molten stage i.e. before complete solidification, due to surface tension, a system with minimal surface area is obtained. In computer simulation the random foam structures are generated similarly by minimizing the surface area (Kraynik, 2006;Kraynik et al., 2003Kraynik et al., , 2004Jang et al., 2008). The numerical procedure of generating realistic open-cell solid foams is shown in Fig. 1 and discussed in detail here.
The first step of the procedure is the generation of a random close packing (RCP) consisting of mono-disperse spheres (Fig. 1a). These RCPs are packings one would obtain in reality by pouring spheres of an equal size in a vase and subsequently vibrating the vase. Since our aim is to generate spatially periodic random mono-disperse foams, it is logical to start with generating a spatially periodic RCP. The techniques mentioned above do not result in such a periodic RCP. The method used in this study is based on the Lubachevsky-Stillinger (LS) algorithm (Lubachevsky and Stillinger, 1990). The method starts with the insertion of a predefined number of small spheres randomly into a periodic box or unit cell. The diameter of the spheres is chosen sufficiently small to ensure that inserted spheres do not overlap with others. After the insertion of all spheres is finished, a growth process is initiated, i.e. the spheres grow with a predefined growth rate until a pair of spheres touches each other. When spheres touch, they are moved by an elastic collision as observed in regular event based molecular dynamics simulations. This growth and collision process continues until the desired packing density is achieved. Care must be taken when picking a growth rate. RCPs are only formed when the growth rate is chosen high enough, i.e. the spheres growth rate is so high that they get ''jammed" before they have time to settle in a crystalline lattice. In the present case we intend to create random foams of $ 27 cells (restricting the computational expenses for flow simulations), hence we have used 27 spheres in a periodic domain. A growth rate of 0:9 is used to obtain a sufficiently high randomness and finally it produces a periodic packing of e ¼ 0:62. The next step is to compute a spatially periodic Voronoi tessellation ( Fig. 1b) based on the centers of the spheres that make up the previously generated RCP unit cell. A Voronoi tessellation is a partitioning of space into Voronoi cells. These Voronoi cells are constructed around predefined seeds, in this case the centers of the spheres in the RCP. The partitioning is done in such a way that all possible points enclosed within a certain Voronoi cell are closer to the corresponding seed than to any other seed. As a result, the closest pair of seeds corresponds to two adjacent cells in the Voronoi tessellation and the line joining those two seeds is always perpendicular to the corresponding face of the Voronoi cell. There are several open source software packages available for the computation of Voronoi tessellations based on the coordinates of the seeds. The open source software library Voro++, developed by Rycroft (Rycroft, 2009), was used in the current work.
The Voronoi structure is used as an input to create a dry foam (Fig. 1c). The dry foam, otherwise known as soap froth, is the structure that essentially would remain when the liquid volume fraction of a liquid foam is decreased to zero. In this dry limit the ligaments reduce to curved edges and the films to simple surfaces of minimal area, all together forming a system of polyhedran which represents the skeleton of the closed-cell foam. Empirical work on soap froth (dry foam) by Plateau (Plateau, 1873) suggests that under a static equilibrium condition it must satisfy Plateau's Laws: (a) constant mean curvature everywhere along the same piece of film, (b) three films meet at an edge under equal dihedral angles of 120 and (c) four edges meet at a vertex under equal tetrahedral angles of 109:47 . At first sight the Voronoi tessellated unit cell, already looks like a piece of dry foam, however the tessellation is not in mechanical equilibrium. Furthermore, the volumes of the Voronoi cells are not equal as typically observed in a mono-disperse piece of dry foam. To generate a stable mono-disperse dry foam, the Voronoi cells are relaxed to attain mechanical equilibrium (minimal surface energy). An open source software package Surface Evolver (Brakke, 1992) is used for this purpose. The Surface Evolver is specially designed to compute the equilibrium shapes of liquid surfaces under the influence of surface tension or other body forces. Next to this, the user can also work with periodic domains and has several features to handle complex topological transitions as observed in real dry foam formation. The iterative relaxation process adjusts the vertex positions, turns straight edges into curved edges, and flat films into curved films in such a way that the total surface area is minimized. The process inherently ensures that Plateau's laws are satisfied. After surface relaxation, the surface energy can be reduced even further by annealing the relaxed mono-disperse dry foam. During the annealing process the dry foam unit cell is subjected to large deformations by means of shearing the unit cell. Subsequent relaxation of the deformed dry foam unit cell as usual, induces topological transitions. Performing several deformation-relaxation cycles in all Cartesian directions will ultimately result in a stable realistic mono-disperse dry foam with minimal surface energy. The annealing process was continued until no significant decrease in surface energy was observed anymore. Starting with a certain Voronoi tessellation results in a unique random mono-disperse dry foam. The interested reader may refer to the intensive work by Kraynik et al. (2003Kraynik et al. ( , 2004) on formation of dry foams.
Wet foamsare generated from the dry foam through a second stage in the Surface Evolver. All the trivalent edges of the periodic dry foam unit cell are replaced with triangular ligaments and the vertices at which those edges intersect are replaced with octahedra (Fig. 1d). Half of each octahedron's faces are omitted and the resulting holes are connected with the ends of the triangular ligaments. All together, this system of connected ligaments and octahedra forms a single liquid containing body, a network of Plateau borders (Plateau, 1873), closely resembling the structure observed in most real open-cell solid foams. As no mass can leave the periodic domain during further processing, the final porosity of the open-cell solid foams for flow simulations will depend on the total volume of the liquid body inserted in this step. Using a higher (lower) liquid fraction results in thicker (thinner) Plateau borders, and in open-cell solid foams of lower (higher) porosity. The same random dry foam can be used as a precursor to generate a range of realistic open-cell solid foams with different porosities.
Finally, the wet foam unit cells have to be turned into solid foam unit cells. Again Surface Evolver is used to minimize the surface area of the wet foam (with ligaments of triangular crosssection). The idea here is to minimize the surface area of the wet foam while keeping the Plateau edges fixated. If the edges of the foam skeleton (i.e. Plateau borders) are not kept fixed, in the surface energy minimization process, the wet foam will try to attain a spherical lump of liquid. Unfortunately, the complete fixation of the Plateau edges is not possible due to numerical constraints (Brakke, 1992). Hence, simulations are performed with a small relaxation time and as a result the basic cell skeleton evolves to some degree in the process. In the present study, a total of five different foam samples of porosities e ¼ 0:877; 0:891; 0:910; 0:929 and 0:948 are generated with more or less the same skeleton structure and in the subsequent section they are denoted as P1 to P5, respectively.  rial and production method. In general, ligaments of circular crosssections (e % 0:85-0:9) and convex-triangular cross-sections (e % 0:9-0:94) are commonly observed in real open-cell foams, whereas for highly porous metallic foams (e % 0:97) the crosssection of the ligaments are concave-triangular (Bhattacharya et al., 2002). In the cross-sectional view of the present foam samples ( Fig. 2) similar geometrical shapes can be observed. Also, from the cross-sectional view of the foam samples one can observe that the structures are periodic.
The specific surface area (S v ), i.e. the available solid surface area (A S ) per unit volume (V) of the foam is an important parameter for heat and mass transfer. In the literature several empirical (based on actual measurements) and semi-empirical models (based on idealized foam geometry) are available to compute S v . It is interesting to compare the specific surface area of the present foam samples with real open-cell foams along with the Kelvin cell representation of a foam unit cell. In this analysis the specific surface area (S v ) is non-dimensionalized by the average cell diameter Hence, c can be considered as a shape factor of foams that indicates the morphological difference between two foams of the same porosity. Fig. 3a and b shows foam cells extracted from the generated foam geometries of e ¼ 0:877 and 0:948, respectively. Fig. 3c shows the idealized Kelvin cell of e ¼ 0:95. The cell diameters (d cell ) for random foams are calculated by a trial and error method where a sphere is initially placed inside the cell and then its diameter is increased until it properly touches the ligaments. In this way all the 27 cells are checked and an average d cell is calculated. In Fig. 3a-c the measuring spheres (of diameter d cell ) are also shown. Finally, Fig. 3d shows the variation of c with e where for both the Kelvin cell geometry and the present random foams, c decreases with increase in e (i.e. with decrease in solid volume in the unit cell). The specific surface area is higher for the present foam samples compared to idealized foam samples. The experimental data (Iasiello et al., 2014;Calmidi and Mahajan, 2000;Bianchi et al., 2015b;Kumar et al., 2014;Giani et al., 2005;Joen et al., 2010;Schmierer and Razani, 2006) for c (i.e. S v ) are very scattered, however they lie in close proximity to both the realistic and Kelvin's foams. In real solid foams one may find several defects like, closed holes or lumps of solid which change S v . Also, polydispersity in the size distribution of the foam cells alter the results. Moreover, uncertainties in the measuring technique may result in large deviations.
Random poly-disperse wet foams are probably better in representing real open-cell solid foams than random mono-disperse wet foams. One can reconstruct poly-disperse foam samples by using poly-disperse RCPs. However, it is numerically challenging to create periodic poly-disperse structures. Moreover, due to finite number of foam cells statistical variability arises. In order to obtain good statistically averaged closure terms, for poly-disperse foam geometry requires very large number of foam cells. It will increase the computational cost drastically. For this reason we restrict this work to mono-disperse solid foams.

Flow and heat transfer through open-cell foams
After constructing five realistic foam structures of different porosities (e ¼ 0:877-0:948) the next task is to perform heat and fluid flow simulations. The random foam geometry is stored as a triangulated surface mesh using the STereo Lithography (STL) file format and imported in the periodic CFD domain. For flow through random porous media there exists a large variety of definitions for the length scale, e.g. pore diameter, equivalent spherical diameter, ligament diameter, square root of permeability, etc. In this work we define the length scale/hydraulic diameter based on an equivalent spherical diameter as: where V S and A S represent the solid volume and solid surface area in the domain, both are calculated numerically from volume and surface integrals. The superficial velocity (U) is used to define the Reynolds number (Re ¼ qUd eqv =l). Starting from creeping flow, Re is varied up to $ 500 and in this range in total nine values of Re are considered for foams of five different porosities. The time step size is set on basis of the Courant number, allowing a maximum value of CFL < 0:5. Starting from an initial zero velocity field, the physical time to reach steady state increases with decreasing kinematic vis-cosity (m ¼ l=q), however not in a linear fashion (Hill et al., 2001).
To decrease the total physical time to reach a steady state (or quasisteady state for high Re cases) the value of m should be high and at the same time if we use very high values of m, to maintain a specific Re the velocity should be high, which demands smaller time step size for proper temporal convergence. For this large parameter space it is very important to select the proper fluid properties to optimize the total computational time. We have used l ¼ 0:003 Pa s for all cases, and q ¼ 10 kg=m 3 for Re % 0:1; 1 and 5; q ¼ 150 kg=m 3 for Re % 10; 20 and 50 and q ¼ 1000 kg=m 3 for Re % 100; 200 and 500; and as a result for all the cases the number of time steps required to reach a (quasi) steady state is found 10 5 . A grid of 420 Â 420 Â 420 cells is used which ensures a minimum of $ 18 grid cells across the ligaments. All the heat transfer calculations are performed with C p ¼ 10 J kg À1 K À1 and k ¼ 0:03 W m À1 K À1 , which assures Pr ¼ 1.
At the surface of the foam a constant wall temperature is specified. A temperature-periodic boundary condition (Eq. (4)) is used for Re % 10 to 500 across the periodic faces in the macroscopic flow direction. For Re % 0:1; 1 and 5, a Danckwert type boundary condition (Eq. (5)) across the periodic faces in the flow direction is used. All boundaries parallel to the flow direction (i.e. y and z boundaries) are subjected to a zero heat flux condition. Few additional simulations are performed using the P3 foam sample to check whether the foam sample is isotropic, homogeneous and the coordinate system is aligned along the principal axes of the porous medium. Using P3 foam we have also performed simulations when the macroscopic flow is along the y-and z-direction to the calculate permeability tensor at creeping flow condition. We found that for both creeping flow as-well-as flow at Re % 200, the permeability tensor is symmetric and almost diagonal (off diagonal values are less than 1:5% of the diagonal values). The permeability for all three directions is almost the same ($ 1% variation). This validates that the present numerically generated random foam samples are isotropic as-well-as homogeneous porous media.

Velocity profile
At very low Reynolds number, i.e. for Stokes flow, the pressure force is balanced by the shear stresses induced by the solid surfaces. When the Reynolds number is gradually increased fluid inertia starts contributing: this is known as the transition from the so called Darcy flow regime to the Forchheimer flow regime. The exact limiting value of this transition depends on the structure of the porous medium and definition of Re. In Fig. 4 simulation results are shown at the mid-plane of the P5 foam sample with varying Re. When the flow is inertia dominated, the maximum velocity appears at the zone of minimal flow cross-section area. However, in the viscous dominant Darcy or creeping flow regime the maximum velocity is formed where the flow area is farthest away from the solid walls. Fig. 4a reveals a similar behavior. With increasing Re, the zones with the maximum velocity contours are located further downstream and have a lower strength -this indicates the pronounced effect of boundary layers being formed at the solid walls and the onset of the nonlinear inertial contribution in the macroscopic flow features. With further increase in Re, a jet like feature appears in the parts with the smallest open area available for flow (Fig. 4b and c). This is very similar to the inviscid core formed in channel flows. Dybbs and Edwards (Dybbs and Edwards, 1984) experimentally reported a similar flow feature in porous media made of spherical and cylindrical packings. They termed this jet-like feature the "inertial core" and identified this as the main source for the non-linear relationship between pressure drop and flow rate. The present numerical results suggest that this is indeed the case. For this particular porosity, for Re % 20 or higher we observe recirculation zones. With increasing Re, the volume of the recirculation zones increase which is very similar to flow past a body in an unconfined flow domain. However, due to the presence of another solid object downstream, the strength of the recirculation is less compared to the case where a solid body is placed in an unconfined domain. For this particular case, unsteady behavior is first observed at Re % 100 in the form of a small oscillation of the inertial core. At Re % 200 and Re % 500 the flow becomes completely unsteady, the stable inertial core breaks down and the vigorous fluid motion leads to good mixing throughout the domain (Fig. 4d).
A powerful way to visualize and study a 3D velocity field in a complex porous medium is to analyze the probability density (PD) of the normalized interstitial velocity components. Fig. 5 presents the probability density of velocity components for the P5 foam sample. The macroscopic flow (i.e. U) is along the xdirection. For the present homogeneous isotropic porous media, qualitatively, both PDðu y =UÞ and PDðu z =UÞ plots show very similar behavior. Hence, only one velocity component perpendicular to macroscopic flow direction is shown. By plotting the PD, one can study the variance or spread of the velocity field, and also identify the degree of recirculation.
First, we focus on analyzing the PDðu x =UÞ for varying Re (Fig. 5a). One of the distinguishable features of PDðu x =UÞ is the sharp peak and discontinuity at u x ¼ 0. This behavior has also been observed by Maier et al. (1998) for fluid flow through a packed bed of spheres and follows from the no-slip condition at the solid surface. The effect of Re is quite evident in the normalized velocity distribution. Except at creeping flow (Re % 0:01), for all the cases the presence of recirculation zones can be identified by observing the negative tails in the PDðu x =UÞ plots. The span of the negative tails (i.e. recirculation zone) increases with Re, as expected for any closed system. It indicates a large mixing for high Re flows. At the positive side of the curves, we see a saddle region (0 < u x =U < 0:8) and a clear peak (u=U % 1:5). The saddle region of PDðu x =UÞ indicates a wider distribution of u x and one can correlate this with the formation of boundary layers at the solid walls. One can correlate the peak in the PDðu x =UÞ curve with the inertial core, where the velocity is comparatively high however the spread or the variance of u x is low. With increasing Re the spread of the saddle region decreases indicating a formation of thinner boundary layers. Because the PDðu x =UÞ has a unit integral (area under the curve), with increasing Re the height of the peak also increases. The maximum magnitude of positive u x =U is invariant with Re and for the present very open cellular porous structures we have not observed any long positive velocity tail, often observed for very closed systems (Maier et al., 1998).
Due to the isotropic structure, the PDðu y =UÞ is almost symmetric for the entire range of Re (Fig. 5b) indicating net (nearly) zero macroscopic velocity in the y-direction. The sharp peak and discontinuity at the u y ¼ 0 suggest (near) zero velocity very close to the solid surface. With increasing Re the variance of u y =U also increases to facilitate a higher flow recirculation. The span of both the positive and negative tails for u y increases with Re. The effect of porosity on both PDðu x =UÞ and PDðu y =UÞ is very small, at least in the present porosity range, and hence it is not shown.

Non-dimensional drag force
It is very common to write the pressure loss in a porous medium as the sum of two terms: viscous drag due to the presence of the solid walls, proportional to fluid velocity and viscosity, and form drag, proportional to the density and velocity squared. We start our analysis from the ansatz form of pressure drop (Dp) inside any porous medium of length L: To analyze the resulting drag force we use the viscous-scale to non-dimensionalize the drag or pressure gradient: As a result, Eq. (9) simplifies to a linear function of Re: One of the main advantages of the current non-dimensional form is that the viscous and inertial contribution can be easily considered separately. The most well-known formula for predicting the pressure drop in a randomly packed bed of particles is the Ergun correlation (Ergun and Orning, 1949) which states: where A and B known as Ergun constants. Universal values and/or form of the Ergun constants have been a subject of considerable debate. Several studies on packed beds filled with particles with different shape and size suggested that actually an Ergun-type of pressure drop correlation ought to be determined empirically for each configuration. A large number of studies in the literature only focuses on finding the values of A and B considering the same porosity functions ð1ÀeÞ 2 e 3 and ð1ÀeÞ e 3 for the viscous and inertial contribution. For creeping flow, Ergun adopted the porosity function ð1ÀeÞ 2 e 3 from the Carman (1937) correlation. The Carman-Kozeny correlation is only valid for creeping flow which suggests A ¼ 180. For close packing of spherical and non-spherical particles (e % 0:3-0:5) the porosity dependency can be derived using the analogy of flow through parallel pore conduits (Niven, 2002). However, for highly porous Here, we will first focus on creeping flow. Fig. 6 shows the variation of non-dimensional drag with solid fraction (1 À e) for five different foam samples along with the Ergun equation (Ergun and Orning, 1949), Carman-Kozeny equation (Carman, 1937) and drag correlation derived for idealized foam geometry (Eq. (11)). With increase in the total solids volume fraction (1 À e), the total drag force increases exponentially due to the rapid increase in solid surface area per unit volume of the system. A few additional simulations have been performed using the P3 foam, but with Re % 0:01 and Re % 0:05. The extrapolation of the drag value to Re ¼ 0 is invariant with the present results for Re $ 0:1 and deviation in the nondimensional drag between Re % 0:1 and Re % 1 is only 1%. This suggests that the simulations are performed in the Re ! 0 limit. To fit the simulation data we use a correlation with a similar form as for the idealized foam (Das et al., 2017b), using however a different constant. The fit obtained from nonlinear regression is given by: The maximum deviation of this fit is less than 1:2% where the average deviation is only 0:65%. Fig. 6 shows that both the Ergun and Carman-Kozeny correlations, under-predict the pressure drop for solid foams. However, the present simulation results very closely resemble the drag correlation derived using the Kelvin cell geometry.
With increasing Re, the microscopic inertial forces related to the local acceleration of the fluid particles increase. This inertial contri-bution in the overall pressure drop in porous media is expressed by a U 2 term by almost all authors, and can be justified by the timeintegration of the average 1D macroscopic advection term (u x Á ru x ). Fig. 7 shows the variation of the non-dimensional friction factor (f) with varying Re for different porosities (e). For proper visualization the low to moderate Re cases and the high Re cases are plotted separately. In this figure, the continuous smooth lines are representing a linear-fit. The slopes of the f versus Re curves are constant for the entire range of Re and only a function of porosity. Table 2 contains the values of b for the foam of different porosities obtained from nonlinear regression for the entire range of Re. It shows that for the entire range of porosities, b, and hence f, increase rapidly with the solids volume fraction (1 À e). To find the curve-fit correlation for the variation of b with e we use porosity function ð1ÀeÞ 1:18 e 4:39 derived using the Kelvin cell geometric representation of foam but for a large range of e. However, a different prefactor is used to reduce the error between the curve fit and actual simulation data. The following correlation is finally proposed for the non-dimensional drag force: For the entire range of Re and e, the average error between the proposed correlation (Eq. (13)) and discrete data is less than $ 5%, while the maximum error is $ 13% (only for few cases).
All the reported studies for fluid drag inside the porous media suggest that the non-dimensional drag (i.e. parameters a b in Eq. (9)) is independent of Re irrespective of the geometrical morphology of the solid structure. However, it is debatable whether b is also independent of Re, and only a function of e. Several experimental and numerical studies on moderate and high Re porous media flows show that b is not connected with Re. However, the numerical results by Hill et al. (2001) for ordered arrays of spheres show that the linear dependency of f on Re occurs only at moderate

Table 2
Slopes of the f versus Re curves, i.e. b with e for the present numerically generated foam samples. and high values of Re. They showed that a Re 2 term appears in the non-dimensional drag when the Reynolds number is small but finite, where as f $ Re when Re > 20. In our previous work (Das et al., 2017b) on fully resolved simulation for the idealized foam geometry, we observed similar behavior: f $ Re when Re > 20. At Re ¼ ½0; 20 we found that b is a function of both Re and e. It is very interesting to observe that for the random porous medium f ðReÞ appears linear including at low (finite) and moderate Reynolds numbers. Fig. 7 shows that linear fit very closely resembles the discrete data points at both moderate and high Re cases. Finally, we compare the non-dimensional drag calculated from the present random foam geometries with the drag correlation derived from idealized foam geometries (Das et al., 2017b) in a log-log parity plot (Fig. 8). The maximum deviation between the two different foam samples at same Re and e is less than 15%. It indicates that the idealized Kelvin cell representation for opencell foams very closely approximates the cellular structure of the random foams.

Temperature profile
We consider solid to fluid convective heat transfer in a (quasi)steady flow past the periodic foam structures when the solid surface is subjected to a constant wall temperature. The main aim here is to obtain the heat transfer coefficient at a thermally fully developed condition. The advantages of establishing a thermally developed condition is that it provides statistically more consistent results as there is no entrance length effect. In the direction of the flow the average temperature increases (or decreases) due to fluid heating (or cooling) by the porous structure. However, at the thermally fully developed condition, the non-dimensional temperature field (H) repeats its shape in the streamwise flow direction. Fig. 9 shows the non-dimensional temperature contours for P5 foam at different Reynolds numbers. One can clearly observe the periodic nature of H for Re % 10 and Re % 50. There exists a direct relationship between the velocity profile and convective heat transfer. The non-dimensional temperature contours, qualitatively, follow very similar trends to the velocity contours (compare Fig. 4b and c with Fig. 9a and b). At high Reynolds number, the overall heat transfer enhances because of vigorous flow structures ( Fig. 9c and d). At both Re % 200 and Re % 500 the flow is unsteady. However, the magnitude of fluctuation is very small (maximum variation in H is AE1%). With increasing Re the maximum value of H decreases.
Similarly to the velocity, we plot the probability density of the non-dimensional temperature (H) to analyze the 3D temperature profile. Fig. 10 shows the probability density (PD) of H for the P5 foam sample. All curves start from H ¼ 0 indicating near wall cells that attain the wall temperature. The probability density at H ¼ 0, is higher at low Reynolds number and decreases with increase in Reynolds number. For all the high Re cases, PDðHÞ exhibits a saddle zone, a sharp peak and a positive tail. The saddle zone in the PDðHÞ plot indicates the temperature distribution in the thermal boundary layer which is not pronounced at Re % 10. The maximum value of H in the computational domain always appears away from the solid, in the core section of the flow and decreases with increasing Re. The spread or the variance of H decreases with increasing Re, and since the area under the PDðHÞ curves are constant, the height of the peak increases when the spread becomes thinner. This behavior shows improved mixing (i.e. more homogeneous temperature distribution) with increasing Reynolds number.

Nusselt correlation
In volume-averaged macroscopic models for heat flow through porous media, the important 'unclosed' term is the heat transfer coefficient between the solid and fluid. The accuracy of such averaged (both spatial and temporal) coarse grained models highly depends on the correctness of the heat transfer coefficient. In order to find the heat transfer coefficient (h) we calculate the total heat flow ( _ Q ) from solid surface to the fluid. The solid-fluid temperature Fig. 8. Comparison between the non-dimensional drag force for the random foam geometry (f random ) and the idealized foam geometry (f ideal ).
difference (DT) is the main driving force for the heat transfer and _ Q ¼ h Â DT. In the present study we have defined the temperature differences for h calculation as, where T S is the uniform solid surface temperature and hT f i is the volume-averaged fluid temperature. Other definitions of the temperature difference can also be considered. However, the present definition for DT is more consistent considering its application to volumeaveraged models (Whitaker, 1998). The logarithmic mean temperature difference (LMTD) is very often used in convective heat transfer through porous media (Das et al., 2017a) and in the present system both definitions provide very close values for heat transfer coefficient. The variation of Nusselt number (Nu ¼ hd eqv =k) with Reynolds number (Re) for different foam samples is presented in Fig. 11. The heat transfer increases with Re (due to increased fluid mixing) and with solid fraction (due to increased solid surface area). Historically, extensive experimental and numerical studies have been conducted to study solid-fluid heat transfer processes. Based on theoretical grounds one can expect the Nusselt number to depend on the Reynolds number (flow property) and Prandtl numbers (fluid property). Furthermore, it will depend on a geometrical parameter since the flow field is also strongly influenced by the geometry of the medium. In the current work it is expected that Nu ¼ f ðe; Re; PrÞ. One of the well known forms to express the non-dimensional Nusselt number reads: where the constants a 0 ; a 1 and a 2 are a function of porosity and solid structure. For convective heat transfer at finite Reynolds number the power of the Prandtl number, i.e. a 3 is proposed as 1=3 in most of the studies. Theoretically, it can be derived using boundary layer theory. On the other side, for very low Reynolds number flows, the Nusselt number is often expressed in term of the Peclet number (Pe ¼ Re Â Pr).
We initially tried to fit all our simulation data in the form of a single generalized Nu correlation presented in Eq. 15. However, we experience a large deviation between the fit and the discrete DNS results at the low Re cases. Another well known equation for the fluid-solid convective heat transfer at packed bed and fluidized beds is given by Gunn (1978) as, Nu Gunn ¼ 7 À 10e þ 5e 2 À Á ð1 þ 0:7Re 0:2 Pr 1=3 Þ þ 1:33 À 2:4e þ 1:2e 2 À Á Re 0:7 Pr 1=3 ; for 0 < Re < 10 5 and 0:35 < e < 1. Gunn's correlation covers a large range of Re and e; and the first term indicates conduction from a sphere to stagnant fluid (when Re ! 0 and e ! 1; Nu ¼ 2). Motivated by Gunn's correlation and other correlations (Sun et al., 2015;Tenneti et al., 2013) In the present study we have performed all our simulations at a fixed value of the Prandtl number (Pr ¼ 1) due to numerical constraints at higher grid resolution. Hence, following Gunn (1978), we use Pr 1=3 term to account for the fluid property dependency of Nu. In Fig. 11, the continuous lines represent the Nusselt   correlation (Eq. (16)). The maximum deviation between the correlation and the discrete data points is less than 7%, whereas the average deviation is only 2% for the entire range of Re and e.
At low Pe axial conduction is not negligible and one cannot use periodic boundary condition for nondimensional temperature. For the present steady state thermally fully developed flow, at Re ¼ ½0; 5 we use Danckwerts' boundary condition for energy at the domain boundaries in the streamwise periodic flow direction. On the other hand, at Re ¼ ½10; 500 we use scaled periodic temperature conditions discussed earlier. Fig. 11 shows a smooth monotonic variation of Nu with Re for the entire range, which indicates consistency in the Nu calculation for both high and low Re cases.

Conclusions
Five different periodic foam samples of varying porosities (e ¼ ½0:877; 0:948) have been generated by numerically mimicking the actual foam formation process (minimizing surface area). The step-by-step procedure for generating the periodic foam geometries has been presented. To the best of our knowledge, it is for the first time that realistic random periodic structures have been used to perform flow and heat transfer simulations. The specific surface areas of the generated foams of different porosities have been compared with real foam geometries showing a reasonable agreement.
A versatile sharp interface second order accurate immersed boundary method (Das et al., submitted for publication) has been applied to simulate flow and heat transfer through periodic foam geometry. The Reynolds number (Re) has been varied from Re % 0 (creeping flow) to Re ¼ 500. For the entire range of Re, flow has been solved using periodic boundary conditions across the domain boundaries. We have considered a constant wall temperature at the fluid-solid interface while solving energy equation. We have used Danckwerts' boundary condition at the domain boundaries in the streamwise periodic flow direction when Re ¼ ½0; 5 (dominating axial conduction). On the other hand, for higher Reynolds numbers, scaled periodic temperature conditions has been used that assumes self-similar temperature profile (negligible axial conduction) across the periodic faces. Finally, a drag and Nusselt correlation have been proposed. A detailed analysis is presented on the local velocity and temperature field for the fluid-solid interaction at complex cellular porous medium.