Pore-Scale Numerical Investigation of Pressure Drop Behaviour Across Open-Cell Metal Foams

The development and validation of a grid-based pore-scale numerical modelling methodology applied to five different commercial metal foam samples is described. The 3-D digital representation of the foam geometry was obtained by the use of X-ray microcomputer tomography scans, and macroscopic properties such as porosity, specific surface and pore size distribution are directly calculated from tomographic data. Pressure drop measurements were performed on all the samples under a wide range of flow velocities, with focus on the turbulent flow regime. Airflow pore-scale simulations were carried out solving the continuity and Navier–Stokes equations using a commercial finite volume code. The feasibility of using Reynolds-averaged Navier–Stokes models to account for the turbulence within the pore space was evaluated. Macroscopic transport quantities are calculated from the pore-scale simulations by averaging. Permeability and Forchheimer coefficient values are obtained from the pressure gradient data for both experiments and simulations and used for validation. Results have shown that viscous losses are practically negligible under the conditions investigated and pressure losses are dominated by inertial effects. Simulations performed on samples with varying thickness in the flow direction showed the pressure gradient to be affected by the sample thickness. However, as the thickness increased, the pressure gradient tended towards an asymptotic value.


Introduction
Open-cell metal foams have unique properties such as high surface area per unit volume, low density, high stiffness and good energy absorption, making them well suited to a diverge range of engineering applications (Banhart 2001). Pore-scale flow data are of considerable interest for many industrial applications, but obtaining it is far from trivial (Bear and Corapcioglu 2012). The foam's intricate geometry makes experimental acquisition of detailed flow data troublesome, and normally only macroscopic quantities such as pressure drop can be measured, drawing many researchers to employ numerical pore-scale simulations in order to overcome such limitations.
There are essentially two classes of numerical approaches employed in pore-scale simulations. The first is network modelling, where a lattice of large pores connected by narrower throats is used to represent the pore space, and the relevant transport equations are computed on top of that network (Blunt 2001;Blunt et al. 2002;Mukherjee et al. 2011). The second is direct simulation, where the flow governing equations are computed directly in the pore space geometry using methods based on first principles. Such type of approach requires an explicit representation of the solid matrix geometry. This paper is concerned with direct pore-scale simulations performed in commercial open-cell metal foams.
Two methods are normally employed to represent the pore-scale geometry: idealized structures or imaging techniques. Both approaches have been used in the past to represent open-cell foam geometries at the pore-scale. The Kelvin cell (Thomson 1887) and Weaire-Phelan structure (Weaire and Phelan 1994) are two idealized structures commonly employed because of their morphological properties akin to foams. Boomsma et al. (2003) investigated laminar water flow through open-cell metal foams using a single Weaire-Phelan unit cell with periodic boundary conditions as the computational domain. Simulation results underpredicted pressure gradient by 25% when compared to experimental data, a fact that was attributed to the lack of wall effects in the computational domain. Krishnan et al. (2006) performed a direct numerical simulation (DNS) of single-phase laminar flow through an idealized periodic bodycentred-cubic (BCC) structure based on the Kelvin cell. Simulated values for pressure gradient and effective thermal conductivity were compared against experiments and semi-empirical models for porosities and showed good agreement for porosities over 0.94. Turbulent flow at the pore-scale was investigated by Hutter et al. (2011), where a large eddy simulation (LES) was carried out in a streamwise (flow direction) periodic structure akin to a metal foam. Simulation results agreed well against experiments but proved to be very sensitive to small changes in geometrical features of the computational model. The use of idealized structures offers a relatively simple way of representing the foam geometry and have the advantage of being periodic in space. On the other hand, they do not have the random disorder present in real foams.
Imaging techniques such as X-ray micro-CT or magnetic resonance imaging (MRI) scans, allied with advances made in 3-D image processing techniques, have made it possible to obtain very realistic digital representations of real porous samples. Additionally, several macroscopic properties of the porous matrix can be directly computed from the image datasets. However, X-ray micro-CT and MRI scans are often time consuming to acquire. Normally, the size of the scanned sample has to be defined in a way that it has statistically representative morphological characteristics. A CT-based image processing technique was developed by Montminy et al. (2004), to generate accurate 3-D digital representation of real polyurethane foam samples. Measurements carried out on the computer generated models were able to determine several morphological features such as strut length distribution, number of pore windows, cell volume distribution and pore size distribution. Maire et al. (2007) employed X-ray tomography scans to characterize several samples of cellular solids.
In a series of papers, Petrasch et al. (2007Petrasch et al. ( , 2008a developed a micro-CT-based method to characterize and simulate reticulate ceramic foams. The porosity, surface-to-volume ratio, minimum representative elementary volume (REV) size were directly computed from the CT scans using robust approaches. Morphological measurements performed directly on the tomographic datasets were validated against experimental measurements of porosity and mean pore diameter. The finite volume method was employed to solve the incompressible Navier-Stokes equations under laminar flow at the pore-scale for a ceramic foam. Simulations results were used to determine permeability, Forchheimer coefficient and the interface heat transfer coefficient. The results derived from the pore-scale simulations were used as a reference for comparison against several analytical macroscopic porous media flow models. It was shown that the analytical models were mostly incapable of achieving agreement with values determined by the pore-scale simulations. Magnico (2009) presented an analysis of hydrodynamic properties of open-cell metal foams using microcomputer tomography and used the lattice Boltzmann method to solve laminar flow through the samples. Results showed that using a REV edge length between 2.5 to 4.5 pore diameters provided the best agreement of simulated transport properties against experiments. Habisreuther et al. (2009) analysed the pressure gradient and tortuosity in opencell foams using the finite volume method. A comparison between a stochastic structure based on the Kelvin unit cell and a MRI-generated foam sample was made. Simulations were run under laminar regime, and pressure gradient results were compared against experimental data from the literature, showing good agreement for the MRI-generated geometry. Significant modifications in the idealized structure were necessary in order achieve the same degree of agreement obtained with the MRI-generated sample.
More recently, Ranut et al. (2014) performed micro-CT-based finite volume single-phase flow and heat transfer simulations on open-cell aluminium foams with three different pore sizes in order to determine pressure drop and heat transfer parameters. Simulation results have shown a certain degree of anisotropy on the porous matrix of the foams. The pressure gradient was strongly affected by varying the Cartesian direction of incoming flow.
One aspect which is commonly ignored in most studies is the foam sample size and its effect on the flow characteristics. Entrance and exit effects were observed in pressure drop measurements performed on several commercial open-cell metal foams by Baril et al. (2008). Results showed the pressure gradient to be dependent on the sample thickness in the direction of flow. Smaller thickness resulted in a higher pressure gradient, which decreased as the sample thickness was increased and tended towards an asymptotic value after a certain critical thickness was reached. Similar results were found in experiments carried out by Dukhan and Patel (2010). Transversal size effects were reported by Dukhan and Ali (2012), though presenting a much smaller magnitude than entrance and exit effects.
De Schampheleire et al. (2016) use direct numerical simulations (DNS) to investigate the relative contributions to the pressure drop from the viscous and form drag terms in a open metal foam. They showed that as the Reynolds number (based on the average filament diameter) increased through the laminar, transitional and turbulent regimes, the contribution from the viscous drag decreased dramatically until it became almost negligible.  Fig. 1 and highlights this effect clearly. The open nature of metal foams means that the individual filaments of the foam are able to act as drag elements and the balance between viscous drag and form drag varies with the Reynolds number. Fig. 1 The contribution of the Darcy (viscous) and Forchheimer (form) terms to the pressure gradient for a metal foam (De Schampheleire et al. 2016) This study aims to address some points commonly overlooked in previous studies. Direct pore-scale simulations are performed on several commercial open-cell metal foams with different pore sizes. X-ray micro-CT scans are employed to obtain accurate digital representations of the foam samples and determine porosity, specific surface and mean pore size directly from the tomographic datasets. The investigated flow range is extended here to the turbulent regime, and the feasibility of using RANS models for modelling the turbulence within the pore space is evaluated. Experimental airflow measurements are performed on the samples and used to extract macroscopic transport properties, namely pressure gradient and Darcian velocity. Computational pore-scale samples representative of the experimental ones are defined, and a commercial finite volume solver is employed to solve incompressible airflow within the pore space. The macroscopic pressure gradient and bulk velocity are calculated from the simulations by averaging and compared against experimental data for validation.
Furthermore, the direct pore-scale simulations are used to show that the size of the REV of an open-cell metal foam cannot be solely defined based on solid matrix morphological properties such as porosity, but rather has to be defined as a function of flow characteristics as well. This is shown by pore-scale simulations on samples with varying thickness in the flow direction. This study intends to show that macroscopic transport properties derived from direct pore-scale simulations can be strongly affected by the computational REV size and that the flow-dependent REV can be orders of magnitude larger than the geometrical one, a fact that is often neglected in many similar studies.

Samples Investigated
Five different commercial foams were investigated: one Recemat ® nickel-chromium foam with 17-23 pores per inch (PPI), two Retimet ® nickel foams with 45 and 80 PPI and two Alantum ® Inconel 625 alloy foams, one with a nominal pore size of 450 µm and a hybrid  Figure 2 shows images obtained using a scanning electron microscope (SEM) of the Inconel 625 alloy 450-µm foam and a pure nickel foam, similar to the Recemat foam, which are manufactured via different routes. The web-like solid matrix made up of many struts surrounded by an interconnected pore space, typical of open-cell foams, is clearly seen. The electrodeposition process employed in the manufacture of the pure nickel foam generates a smooth surface. Conversely, the alloy powder spraying process used on the Inconel foams creates a very rough surface with a granular topology, having a larger surface area per unit volume, which can be advantageous in certain applications.

Experimental Data Acquisition
High-resolution X-ray micro-CT scanning was employed for generating realistic pore-scale representations of real metal foam samples. Electric discharge machining (EDM) was used to cut cylindrical samples from the original foam slabs, ensuring clean cuts and avoiding distortion. A desktop cone-beam micro-CT 40 from SCANCO Medical was employed for scanning of the samples. The instrument presents a trade-off between the scanning resolution and size of the samples. The scan resolution, defined as h, is given in voxels, a volumetric pixel which represents the smallest detectable sub-volume in 3-D Cartesian coordinates. The cylindrical samples were placed inside a 8-mm-diameter sample holder, so that the maximum isotropic digital resolution of h = 12 µm could be achieved. This helps minimize the error in exposure and ensures that the foam micro-structure is well represented. With the intent of obtaining a larger sample, a digital resolution of h = 16 µm was employed for the Recemat NC 1723 sample, which enabled the use of a 12-mm-diameter sample holder.
Pressure drop measurements were performed on all the samples investigated under a wide range of air flow velocities, from laminar to turbulent regimes. The test rig consisted of a manual control valve, pressure regulator and filter, a needle valve, a flow meter and the test section (sample holder). The detailed specifications of the parts are found in Oun and Kennedy (2014). The sample holder was designed to hold samples with a diameter of up to 25 ± 0.2 mm (actual flow diameter of 21.183 mm), and a maximum thickness of 30 mm in the flow direction. The pressure transducers were positioned approximately 25 mm upstream and downstream from the sample. Note that the samples used for the pressure drop measurements were considerably larger (in the transversal direction) than the ones employed in the tomographic scans.

Image Processing
Various image processing routines were written in MATLAB using some of its built-in toolboxes, along with a few community contributed functions (available in the Mathworks website) for rendering 3-D volumes, performing morphological measurements and computing the minimum geometrical representative volume (MGRV) size for each sample.
Tomographic datasets are defined as sets of grey values Ψ i jk uniformly distributed in a Cartesian grid with spacing h. The scanning process generates stacks of N z 8-bit grey-scale images of N x × N y pixels each, corresponding to N x × N y × N z voxels. The 3-D matrix of grey values can be regarded as a smoothed representation of the pore and solid space, respectively, with a smoothing kernel of size h, corresponding to the scan resolution.
The pore and solid spaces of porous media have distinct optical densities, which are associated with local X-ray absorption coefficients, reflecting different phases (solid or fluid) within the material. The process of digitally identifying these phases and partitioning them into disjoint segments is referred to as segmentation. However, it is often necessary to perform image enhancement operations before segmentation. The resulting images obtained from the CT scanning process are not sharp enough and can present significant light scattering and noise. These issues are corrected prior to segmentation by applying image intensity adjustments (to enhance contrast), unsharp masking (sharpens the edges) (Sheppard et al. 2004), and 2-D median filtering (to reduce noise). Due to the nature of the images, it has been noticed that careful selection of the image adjustment parameters is necessary, since it can greatly affect the final outcome of the volume rendering procedure.
The image enhancement is followed by segmentation. The solid-pore interface is determined by defining a grey-value threshold level, Ψ 0 . All pixels in the input image with a value greater or equal than the threshold are replaced with the value 1 for the solid space and value 0 for the pore space, yielding a set of binary values Ψ (bin)ijk . The global grey-value threshold level Ψ 0 is determined using Otsu's method (Otsu 1975). This process yields a continuous representation of the solid-pore interface, represented as an iso-surface, Ψ (x) = Ψ 0 . Therefore, the pore indicator function can be written as and is used to convert the grey-scale matrix into a 0/1 binary matrix. Open-cell metal foam commonly has hollow struts due to their manufacturing process, constituting a number of disjoint pore regions not relevant to fluid flow. Therefore, the largest subset of the pore space is computed by the use of a flood fill operation. This largest subset constitutes the main pore space. The remaining smaller disjoint pores are artificially closed by converting their voxels to 1 (solid), since they are irrelevant for the flow. A second flood fill approach is applied on the solid space, ensuring that no disjoint solid regions remain. Subsequently, a 3-D Gaussian smoothing operation is performed on the volume using a convolution kernel of size Δ k = 3h × 3h × 3h, which prevents a stair-stepped solid-pore interface to be generated. Finally, the marching cubes algorithm (Lorensen and Cline 1987) is applied to generate a triangulated representation of the solid-pore interface, which is converted to a stereo-lithography (.stl) format. A 2-D cross-sectional reconstruction of the metal foam geometry is depicted in Fig. 3a, with the highlighted part showing a hollow strut. The final 3-D rendering of a cylindrical sample with a diameter of 12 mm is shown in Fig. 3b.

Macroscopic Morphological Characterization
The tomographic datasets can be used for direct computation of macroscopic morphological parameters such as porosity, specific surface, mean pore and strut diameter and MGRV size for each sample.

Porosity
The porosity, φ, is determined from the tomographic datasets by counting the number of voxels in the solid phase and dividing by the total number of voxels from the binarized datasets, namely The results of porosity measurements are depicted in Table 1. The total porosity also accounts for the pore space within the solid foam struts, whereas the effective porosity accounts for the main pore space only. In the work of Oun and Kennedy (2014), mercury intrusion porosimetry (MIP) was used to determine the total porosity of the Inconel 450-µm sample, which showed a value of 0.8133, presenting a relative error of less than 1% when compared to the tomographic dataset. This result agrees with data from Lee et al. (2010), which reports a CT-computed porosity value of 0.8214 for the Inconel 450-µm foam. The total porosity value of 0.8899 for the Recemat NC 1723 computed from the tomographic dataset compares well against total porosity values in the range of 0.89-0.90 reported by Medraj et al. (2007), where several samples of the Recemat NC 1723 were analysed. The fraction of unconnected pore space was shown to be generally quite small (<1%) for all foam samples. No porosimetry tests were performed on the other foams. However, since the structure of these foams is similar to the others and they were subjected to the same volume rendering approach, it seems reasonable to assume that results reported in Table 1 are within the same level of accuracy as the other foam samples.

Specific Surface
The specific surface is defined as the solid-pore interface surface area per unit volume and is determined by applying the Cauchy-Crofton theorem (do Carmo 1976), using a MATLAB routine. The total specific surface also accounts for the solid-pore interface area within the hollow struts, whilst the effective surface area accounts for the solid-pore interface area in the main pore space only. As general trend, the specific surface increased with decreasing the pore diameter. Differently from the porosity, the foam samples analysed here showed a large range of specific surface values, as Table 1 shows. The difference between total and effective specific surface values is more significant in samples with larger pore size, where the unconnected pore space inside the struts is better captured by the tomography scans.
Furthermore, values of specific surface are dependent on the tomography resolution. The surface roughness observed on the Inconel foams, as depicted in Fig. 2a, is not entirely captured with the tomography resolution employed here, and hence the generated tomographic datasets have considerably smoother surfaces. No experimental measurements of specific surface were performed. Lee et al. (2010) reported a specific surface value of 7088.8 m −1 for the Inconel 450-µm from CT scans, showing relatively good agreement with the values computed here. Geometrical features smaller than the scan resolution tend not to significantly affect the fluid flow, since they tend to be small compared to the boundary layer thickness.

Mean Pore and Strut Size
The pore space of consolidated media such as the ones that will be covered here can be characterized by means of a pore size distribution. The pore diameter at a given position vector x in the porous domain is defined as the diameter of the largest sphere, D p , that contains x and remains entirely within the pore space, U p , without overlapping the solid matrix, U m (Bear 2013), where d p is the pore diameter. Thus, by attaching a pore diameter to each point of the pore space, a pore size distribution can be defined. The pore size distribution defined by Eq. (3) is determined from the tomographic samples by computing the opening size distribution. Opening is an operation of mathematical morphology, consisting of an erosion followed by a dilation using the same structuring element (SE), thus eliminating all geometrical features smaller than the current SE. Morphological opening was carried out by sequentially inscribing spherical SEs of increasing diameter in the pore space of the tomographic datasets and counting the remaining number of (pore) voxels, until no more voxels are left. The opening size distribution is obtained by relating Table 2 Measurements results for the mean pore and strut diameters along with their normalized coefficient of variation, and minimum MGRV normalized edge length where d SE is the SE diameter and f (d SE ) is the opening size distribution. The mean pore size and the coefficient of variation (C V p ), defined as the standard deviation normalized by the mean pore size, are listed in Table 2. Figure 4b depicts the mean pore diameter as a function of the effective specific surface for all foam samples investigated, showing it to be inversely proportional to the specific surface. The only foam sample that deviates slightly from this trend is the hybrid Inconel 580 + 1200 µm, probably due to the fact that it is composed of two foams with different pore sizes merged together. The mean strut size, d s , is computed in exactly the same way as the mean pore size, but instead of performing the morphological opening operations in the void space, they are performed in the solid matrix using spherical SEs of increasing diameter. Therefore, the mean strut size is equivalent to the mean diameter of the maximum inscribed spherical structuring element. Results for d s and its coefficient of variation (C V s ) are listed in Table 2. The mean strut diameter is on approximately five times smaller than the mean pore size on average.

Determination of the Minimum Geometrical Representative Volume
The minimum geometrical representative volume (MGRV) is defined here in as the smallest sub-volume of the open-cell metal foam for which statistical meaningful local average quantities can be obtained. Therefore, the MGRV is defined as smallest subset of the porous medium for which the a macroscopic morphological quantity, such as the porosity φ, is computed for subsequent larger subsets of the porous medium, until it asymptotically reaches a constant value within a band of ±δ. In that sense, the edge length of a cubic sub-volume of the porous medium is defined as (Petrasch et al. 2008b), where U L * is a cubic subset of the porous domain with edge length L * . This method is identical as the one employed by Petrasch et al. (2007), where a large subset of the tomographic data is initially generated, and the porosity φ is subsequently calculated in sub-volumes of increasing size until φ converges around an average value φ ∞ , satisfying the condition stated by Eq. (5).
The MGRV defines the lower bound for the computational domain size. However, even though increasing the MGRV size would not change macroscopic morphological quantities of a given porous medium, the same cannot be said of its flow transport properties. A flowdependent REV size must be computed based on flow characteristics, which requires flow data and is performed on a later stage in the current work. The effective porosity is the macroscopic morphological quantity analysed here. For the limiting case where the sub-volume cubic edge L MGRV is zero, φ eff is either 0 or 1, since the sampling sub-volume is only a point, which can be located either in the solid or fluid phase. As L MGRV increases, φ eff converges around a mean value within a band of ±0.025φ eff . Figure 4a shows φ eff , as function of the sub-volume cubic edge length normalized by the mean pore size, L MGRV /d p , around six different starting locations for the Recemat sample. φ eff converges for L MGRV /d p ≈ 2. Results for the minimum MGRV normalized cubic edge length, L MGRV /d p , are listed in Table 2 for all samples.

Flow Governing Equations
This section presents the governing equations for the flow of a single-phase fluid within the pore space. The Knudsen number, defined as the ratio of molecular mean free path length to a representative physical length scale, can be used to determine whether the flow phenomena lies under the continuum hypothesis assumption. Values of Knudsen number are expected to be <10 −4 for all cases, using the mean pore size as the representative length scale, meaning that the mean free path length of a molecule is some orders of magnitude smaller than the length scale of the problem. Therefore, under the continuum hypothesis, the flow of an incompressible Newtonian fluid in the pore space can be described by the continuity and Navier-Stokes equations.
For an isentropic gas like air, the compressibility of the fluid can be related to the Mach number, M, as where du and dρ represent differential changes in velocity and density. It is assumed that if M < 0.3, the change in density will be less than 10% and in this case, the assumption that air is an incompressible fluid can be used. In all our simulations, the air velocity, even within the pore space, is less than 100 m/s, which is equal to 0.3M. The conservation equations are derived in an Eulerian frame of reference, for an infinitesimal control volume located within the pore space with all the fluxes across its boundaries set in balance. The conservation of mass in given as where ρ is the fluid density, u is the velocity vector, and S m is a source term which accounts for added mass to the continuous phase. Equation 7 is the general form of the mass conservation equation and is valid for both compressible and incompressible flows. The conservation of momentum in a static reference frame is given by where p is the static pressure, τ is the stress tensor, ρg and F are the gravitational and body force and external body forces terms, respectively. For a Newtonian fluid, the stress tensor is defined as where μ is the fluid dynamic viscosity, I is the unit tensor, and the second term on the right hand side accounts for the effect of volume dilation.

Flow Regimes in Porous Media
Not every flow through porous media is laminar. High-speed flow (high Reynolds number) through porous media can occur and lead to the onset of turbulence within the pore space. This is even more likely if the interstitial fluid is a gas and if the porous medium is coarse (high porosity). The present study is concerned with high-velocity airflow through highly porous open-cell metal foams, and the onset of turbulence within the pore space is very likely. The Reynolds number based on the mean pore diameter, Re p , is commonly employed to characterize flow through highly porous media and is defined as with d p as the mean pore diameter and u D as the Darcian velocity, which is defined as the average superficial velocity, computed as where Q is the volumetric flow rate and A f is the cross-sectional flow area. The literature has distinguished essentially four types of flow regimes in porous media (De Lemos 2012; Pedras and Lemos 2001), which based on Re p , are identified as: (1) Darcy or creeping flow regime (Re p < 1); (2) Forchheimer flow regime (1 < Re p < 150); (3) post-Forchheimer flow regime (150 < Re p < 300) and (4) fully turbulent flow regime (300 < Re p ).

Turbulence Modelling
The investigated flow regime is expected to be predominantly within the turbulent flow range. Hence, a model is required to account for the turbulence within the pore space. Most turbulence models rely on what is referred to as the Reynolds averaging technique. As turbulence is characterized by random fluctuations, the instantaneous Navier-Stokes equations can be decomposed into the mean (ensemble-averaged or time-averaged) and fluctuating components. In this way, the pore velocity components, using Einstein notation, are decomposed as where u i and u i i are the mean and fluctuating velocity components, respectively. Likewise, this decomposition can be applied to other scalar quantities where ϕ denotes any scalar such as pressure, energy or species concentration.
Applying this approach to the instantaneous mass and momentum equations and timeaveraging them (and dropping the over-bar on the mean velocity, u) yields the Reynoldsaveraged Navier-Stokes (RANS) equations, written in Cartesian tensor form as Equations 14 and 15 have the same general form of the instantaneous Navier-Stokes equations, with the velocities and other flow variables now representing time-averaged (or ensemble-averaged) values. Equation 15 includes additional terms representing the effects of turbulence, known as the Reynolds stresses −ρu i u j , which require additional equations to close the problem. Turbulence models arise from the modelling of such terms. In the present work, the re-normalization group (RNG) k-ε and the shear stress transport (SST) k-ω models are employed.

The RNG k-ε Model
The RNG k-ε model is based on the work of Yakhot et al. (1992) and has two transport equations, one for the turbulent kinetic energy, k, and other for the dissipation rate, ε. These two transport equations are given as where μ eff is the effective eddy viscosity, G k accounts for the generation of turbulence kinetic energy due to the mean velocity gradients, G b represents the generation of turbulence kinetic energy due to buoyancy, Y M is the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate, and S k is a source term. The terms α k and α ε are the inverse effective Prandtl numbers for k and ε, respectively. The effective eddy viscosity can be computed from k and ε as where C μ is a constant determined empirically. The term C * 2ε is given by with η = Sk/ε, where S is the modulus of the mean rate-of-strain tensor. The constants C 1ε , C 2ε , β and η 0 are determined from analytical and empirical flow data. Details on the computation of all the extra terms can be found in Versteeg and Malalasekera (2007).

The SST k-ω Model
The SST k-ω model by Menter (1994) is a two-equation model that introduces transport equations for the turbulent kinetic energy, k, and the specific dissipation, ω, which can be thought of as the ratio of ε to k. The SST k-ω transport equations are written as where the newly introduced Γ k and Γ ω represent the effective diffusivity of k and ω, respectively, G ω is the generation of ω, Y ω and Y k account for the dissipation of k and ω due to turbulence, D ω represents the cross-diffusion term, and S ω is a source term. The diffusivity terms Γ k and Γ ω are calculated as where σ k and σ ω are the turbulent Prandtl numbers for k and ω, respectively. The eddy viscosity is then computed as where α * and α 1 are closure coefficients and F 2 is a blending function. Details on the derivation of the extra terms are beyond the scope here and can be found in Menter (1994). The model constants are obtained from analytical and empirical turbulent flow data.

Meshing and Computational Domain Generation
The outcome of the tomographic volume rendering process is a triangulated surface of the foam sample with triangular edge length Δ e = h, exported as a stereo-lithography file to a commercial meshing software, namely Ansys Fluent 15 (ANSYS 2013). The complexity of the foam structures required the use of an unstructured grid generation algorithm. A tetrahedral mesh was firstly created and used as a template for creation of a polyhedral mesh. Curvature-based grid-refinement was employed, i.e., the cell mesh size was refined based on the curvature of the foam surface, where a high level of refinement was achieved due to the intricate features of the geometry. Away from the foam surface, the mesh size increases at a rate of 1.3 until a pre-set maximum value, to ensure a good resolution of the flow boundary layer. Accurate computation of the boundary layer is dependent on the mesh resolution next to the foam surface. The numerical meshes were refined to the point where the first cell next to the surface is within the viscous sub-layer (y + < 5). An illustration of the foam surface mesh and a 2-D cut plane of the pore space mesh is depicted in Fig. 5. The computational domain configuration is defined as a rectangular duct consisting of an undisturbed inlet and exit regions and an intermediate section containing the foam 3-D mesh representation. A Dirichlet boundary condition is specified at the inlet, with an uniform velocity profile normal do the boundary such that u · n = u 0 . A no-slip condition is defined at the solid-fluid interface, u = 0. The side walls of the rectangular duct have free-slip conditions, namely u · n = 0 and n · ∇u = 0. A constant zero gauge pressure is defined over the outlet. This is interpreted as the static pressure of the environment into which the flow exhausts being equal to the atmospheric pressure. The turbulence intensity and length scale at the inlet and outlet are set to 5% and equal to the hydraulic diameter receptively. This computational domain configuration is representative of the central part of the experimental samples (away from the walls), and the thickness of the domains in the flow direction was kept as close as possible to the experimental samples. Figure 6 illustrates the computational domain configuration and its boundary conditions. The optimal domain cross-sectional area orthogonal to the flow was determined through convergence studies on cubic sub-volumes of increasing size, similarly to the MGRV calculations described in Sect. 2.5. Computational samples with volumes ranging from d p × d p × d p to 5d p × 5d p × 5d p were simulated at Re p = 100, and the pressure gradient was computed for each sub-volume. Additionally, mesh independence was investigated by performing simulations on grids of increasing levels of refinement at Re p = 100 and computing the pressure drop at each level. Convergence was determined by comparison against the largest crosssectional area domain and highest level of mesh refinement, respectively.
Analysis of the optimal cross-sectional area showed convergence for pressure gradient values when a volume of 4d p × 4d p × 4d p was employed. Grid independence was achieved when an average mesh length scale of 0.03d p was used, which was valid for all samples. The total mesh count of the domains used for comparison against experiments ranged from Fig. 6 Schematic of the computational domain and its boundary conditions 0.2 − 1.5 × 10 7 polyhedral cells depending on the metal foam type. Larger cell counts were found for foams with smaller pore size, due to a greater surface area and thus requiring a larger number of smaller cells next to the foam surface.

Flow Solving Procedure
The flow governing equations described previously are solved with the finite volume method using the commercial solver Ansys Fluent 15 (ANSYS 2013). Airflow is simulated through static and completely rigid foam samples. The fluid is assumed to be incompressible, and simulations are performed in steady state. The SIMPLE algorithm by Patankar and Spalding (1972) is employed for pressure-velocity coupling. A second-order upwind scheme was used for the spatial discretization, and the linearized equations were solved using algebraic multi-grid acceleration. A detailed description of the numerical methods is beyond the scope here and can be found in Versteeg and Malalasekera (2007). A dimensionless convergence criterion of 10 −6 was set for the equation scaled residuals together with the monitoring of convergence of selected field-variables in order to ensure the solution has reached a steady state.

Pressure Drop Relationships
The pressure gradient is one of the mostly easily measurable macroscopic parameters in applications involving flow through porous media. The pressure gradient is assumed to be unidirectional here and is computed for both the experiments and simulations. At very low fluid velocities, the viscous drag dominates the energy dissipation and the unidirectional pressure gradient behaves according to Darcy's law (Darcy 1856) on a macroscopic level, given by Δp where the quantity K is referred to as the medium permeability, and it is related to the geometry of the solid matrix only. L f is the foam thickness in the flow direction. At higher flow velocities, inertial effects start to become apparent and an additional quadratic term was proposed by Dupuit (1863) and Forchheimer (1901), resulting in the nonlinear Forchheimer where F is the Forchheimer coefficient, which is believed to be fixed for a given class of porous media (Dukhan et al. 2014). The derivation of F is cumbersome, and it is usually obtained from best-fit to experimental data. The pressure gradient, defined as the pressure drop normalized by the foam thickness in the flow direction, Δp/L f , is the parameter chosen for comparison between experiments and pore-scale simulations, since experimental values are available for the samples analysed here. Therefore, the pressure gradient in the pore-scale simulations was calculated by computing the area-averaged static pressure at the domain inlet and outlet, respectively.
The experimental values for Darcian velocity were computed by measuring the volumetric airflow rates upstream the foam samples and dividing by the cross-sectional flow area at the measurement location. In order to be consistent with the experiments, the Darcian velocity is specified at the domain inlet in the pore-scale simulations, using the volumetric airflow rate and inlet cross-sectional area. Because the flow is assumed to be incompressible, and the inlet section upstream the foam sample is a clear channel, the Darcian velocity will have the same value as the prescribed inlet velocity for the pore-scale simulations. The pore velocity, u P , is where φ is the porosity. It should be remembered that this is a macroscale definition of pore velocity, the velocity predicted by the CFD is both a Darcian velocity (in the inlet section, for example) and a pore velocity (in the porous medium). Throughout what follows, when a particular case is referred to as u D = 20 m/s, for example, it is referring to the Darcian (or inlet) velocity.

Comparison to Experiments
Preliminary simulations were carried out in order to verify whether the assumption of a steadystate flow is valid under the flow conditions here. Time-averaged unsteady simulations were compared against steady-state flow solutions, and a deviation of less than 3% in the pressure gradient was found. Therefore, all subsequent pore-scale simulations have assumed steady flow, in order to keep computational costs at a minimum. Simulations were carried out using a high-performance computing (HPC) facility, employing normally 24 Intel Harpertown 3 GHz cores. A steady-state pore-scale simulation took between 30 min to 2 h of computing time to achieve convergence depending on the mesh size. A comparison between the experiments and numerical solutions obtained with the RNG k-ε and SST k-ω turbulence models was made using the Recemat foam sample. Pressure gradient results calculated from the pore-scale simulations are compared against experimental values for validation. Note that the computational domain had a thickness of L f = 10 mm, identical to the experimental sample. The pressure measurements were taken under a range of u D = 2.3-26 m/s, and pore-scale simulations were performed for u D = 5-25 m/s. Pressure gradient results for the Recemat foam are plotted in Fig. 7a as a function of u D . Comparison between experiments and results obtained from pore-scale simulations showed larger devia- Δp/L f (Pa/m) The subscript CFD refers to simulations and exp to experimentally-derived values tions at lower Darcian velocities values, u D < 10 m/s (Re p ≈ 800). Results obtained using the RNG k-ε and SST k-ω turbulence models were very similar. By assuming the pressure gradient curves to behave according to the Forchheimer equation (Eq. 26), permeability and Forchheimer coefficient values can be obtained by fitting a secondorder polynomial to each curve. Hence, validation of results obtained from the pore-scale simulations can be done by comparing values of K and F with the ones calculated from experimental data. Results reported in Table 3 have shown that the RNG k-ε model presented a better agreement than the results obtained with the SST k-ω model. At the present stage, there is no simple explanation for the better agreement shown by RNG k-ε model, because detailed experimental turbulence data within the foam pore space would be necessary in order to fully understand the cause for this. However, all subsequent pore-scale simulations results reported here have been obtained using the RNG k-ε model. Experimental and numerical permeability and Forchheimer coefficients were computed from second-order polynomial curve fits applied to pressure gradient curves for all other samples and are reported in Table 3. The quality of the curve fitting procedure was very good, with the square of correlation factor, R 2 , greater than 0.99 in all cases. Table 3 also shows the investigated Re p range for all samples.
The pressure drop experimental measurements for the Inconel and Retimet foams were carried out for u D = 2.3 − 49.6 m/s, except for the hybrid Inconel 580 + 1200 µm, where All the pore-scale simulations were performed for u D = 5−50 m/s; however, calculations for the Inconel 450 µm and Retimet 80 PPI considered an extended velocity range, u D = 1 − 50 m/s. Pore-scale simulations of the hybrid Inconel 580 + 1200 µm were limited to u D = 26 m/s, since the no pressure drop measurements were carried out above that. The thickness of the Inconel 450-µm sample was identical to the experimental ones, except for the hybrid Inconel 580 + 1200 µm, which presented a small discrepancy of 0.40 mm, due to the size of the tomographic dataset. The Retimet foams computational domains were smaller than ones evaluated experimentally due to the limited size of the tomographic datasets The Inconel and Retimet samples pressure gradient curves obtained with the pore-scale simulations agreed relatively well against experimental values as depicted in Figs. 8 and 9. Greater deviation is observed in the Retimet results, probably because of the discrepancy between the experimental and computational sample thickness. Values of the permeability and Forchheimer coefficient for these two samples are reported in Table 3. Larger pressure gradients are observed in samples with smaller pore sizes, which present a greater flow restriction due to their larger surface area.
Measurements have shown the pore size to be inversely proportional to the specific surface. Similarly, the permeability values show a weak decrease with increasing surface area of the foams, with some samples presenting a small deviation from this general trend. This same trend was observed in pore-scale simulation results by Ranut et al. (2014). Figure 10a depicts the CFD permeability values for all foams as a function of the specific surface area. Permeability is also be a function of porosity, but since the samples investigated here have a relatively narrow porosity range, 0.82-0.89, the specific surface area seems to be a more representative macroscopic parameter for comparison purposes. The Inconel 1200 + 580 µm sample deviates from the general trend, probably because it is composed of two foam layers of different pore sizes merged together, therefore having distinct geometrical features compared to other samples. Figure 10b depicts the numerical Forchheimer coefficient values as a function of the effective specific surface area. Again, the Forchheimer coefficient displays a general increase with increasing effective specific surface of the foams. Higher values of the coefficient have been found to be associated with higher pressure gradients and vice versa. The Forchheimer coefficient obtained for the Retimet 80 PPI foam is similar to the one for the Retimet 45 PPI sample, even though the specific surface of the Retimet 80 PPI is significantly larger. These results may indicate that the Forchheimer coefficient is a function not only of the solid matrix but of the flow characteristics as well, since the pressure gradient curves for these two foams are of similar magnitude.
The large discrepancies between permeability values computed from the experiments and pore-scale simulations can be explained by the physical interpretation of the linear permeability term in the Forchheimer equation. The flow regimes investigated here are well above Darcy's, even for the laminar cases, as shown by the Re p range in Table 3. As a consequence, under these flow regimes, the viscous losses, which are represented by the linear term in Eq. 26, are expected to become negligible. This is in agreement with the findings of De Schampheleire et al. (2016) who state that at high Reynolds numbers the contribution from the Darcy term is negligible, with the bulk of the pressure losses being due to the form drag of the filaments.
Under these flow conditions, the quadratic term in Eq. 26, which represents pressure losses due to inertial effects, is dominant. Thus, for the flow conditions investigated here, the Forchheimer coefficient is most representative parameter. Moreover, this also explains the large Forchheimer coefficient values found. Moreover, since the curve fitting procedure is basically searching a 2-D parameter space, there will inevitably be more than one pair of permeability and Forchheimer coefficients as a solution. If, on the other hand, it is assumed that the linear term in Eq. 26 is zero, then the pressure gradient of the foam samples analysed here can be rewritten as Thus, the pressure gradient curves obtained by both experiments and CFD simulations can be fitted using Eq. 28, and Forchheimer coefficient values can be recomputed. This allows for a more meaningful comparison between results obtained from experiments and the pore-scale simulations. A curve fitting procedure assuming the pressure gradient curves to behave according to Eq. 28 is done, and results are depicted in Table 4, which also shows the normalized RMS deviation between experimental and pore-scale simulation data (using the experiments as reference). Again, the quality of the curve fits is very good, with the square of correlation factor, R 2 , higher than 0.99 for all cases.
The agreement between Forchheimer coefficient values obtained from experimental and pore-scale simulations results is significantly better when the viscous term is neglected. The only sample for which pore-scale simulation results showed a relatively high deviation is the Inconel 580 + 1200 µm. A possible reason for this is due to deviations in the geometry of the tomographic and pressure drop measurements samples, caused probably during the merger of the 580 and 1200 µm foams. Nevertheless, the pore-scale simulations were able to show a reasonable level of agreement with experiments. Furthermore, these results give support to the assumption that pressure losses are caused mainly by inertial effects under turbulent flow conditions. Figure 11a depicts a pressure contour plot on a 2-D cut plane for the Recemat foam, at a flow velocity, u D , of 20 m/s (Re p ≈ 1700). A complex pressure distribution is caused by the intricate foam geometry, with stagnation points seen in front of the foam struts. A pore velocity contour plot on a 2-D cut plane parallel to the flow direction is depicted in Fig. 11b. The flow is separated and accelerated through the pores apertures, showing a preferential flow through large passages, with high-velocity streaks seen through the narrowest gaps. The flow behaviour is similar to the flow past multiple bluff bodies in this case, with a low-velocity region present on the wake of every solid ligament. Simulation results obtained with the other samples showed similar behaviour. The pressure gradient behaviour along the sample thickness is more clearly shown in Fig. 12a, where the area-weighted averaged dimensionless pressure profile was computed in multiple cross sections perpendicular to the flow and is plotted against the streamwise position normalized by the sample thickness, z/L f , for multiple inlet velocities (u 0 as the area-weighted average velocity over the domain inlet). The entrance and exit of the foam sample are depicted by the vertical dotted lines. In all cases, the profile shows an almost linear pressure gradient behaviour, with local fluctuations caused by the intricate foam geometry. At the exit of the sample, the local pressure falls below the far field pressure, with the drop becoming more pronounced for higher flow velocities. This phenomenon is likely to be caused by vortices generated at the wake of the foam. The linear pressure gradient behaviour observed here is in agreement with similar studies (Petrasch et al. 2008a;Akolkar and Petrasch 2012). Figure 12b shows the Darcian velocity profile normalized by u 0 . In this case, the Darcian velocity is computed as the area-weighted average pore velocity in multiple cross sections perpendicular to the flow. As the incompressible fluid enters the foam, its bulk velocity is immediately increased due to the reduction in overall cross-sectional area and continues to fluctuate above unity, proportionally to the cross-sectional area variation along the foam. At the exit, the normalized velocity gradually goes back to unity, and a small exit effect is present at the wake of the foam.

Results with Varying Sample Thickness
The effects of varying the thickness of samples in the flow direction are evaluated numerically for the Inconel 450 µm. The 3-, 4-and 7-mm samples are subsets of the 10-mm one, whereas the 20 mm is two 10 mm combined together. Note that this merging procedure will result in a discontinuity on the structure of the foam at the merging position; however, no statistically significant change was observed in the porosity or specific surface area of the 20-mm sample.
It is clear from Fig. 13 that the pressure gradient decreases with increasing sample thickness, a trend which was observed experimentally in other studies (Dukhan and Patel 2010;Oun and Kennedy 2014). The pressure gradient for the 10-and 20-mm samples almost collapse into a single curve for the velocity range investigated, indicating that the critical sample thickness probably lies between these values. A plot of the pressure gradient as a function of the sample thickness normalized by mean pore size is shown in Fig. 14. The same trend observed in Fig. 13 extends across the entire velocity range, being more pronounced at higher fluid velocities. Apparently, the critical thickness lies between 25 and 50 d p .
These results highlight the importance of adequate sample sizing, which should be applied for both experimental and computational samples. This aspect is commonly neglected in many similar studies published in the literature. Furthermore, these results show that a flow- dependent REV size can be in fact many times larger than the REV size determined solely by macroscopic morphological parameters such as porosity (as shown in Sect. 2.5). This also points out to the fact that both experimental and computational samples investigated in Sect. 4.1 are likely smaller than the flow-dependent REV size (except probably for the Inconel 450 µm). However, it is important to remember that the permeability and Forchheimer coefficient values calculated in Sect. 4.1 were only used for comparison between experiments and simulations and should not be taken as definitive values for the open-cell metal foams analysed here.

Conclusions
This study described a standard grid-based pore-scale CFD approach applied to several distinct commercial open-cell metal foams under a wide range of flow velocities, focusing on high-velocity turbulent flow. Tomographic datasets of each foam were acquired and used for characterization of the five foam specimens and generation of the computational domains used in the pore-scale simulations. The volume rendering approach was validated by comparing porosity values obtained from the tomographic datasets with MIP and literature results, which showed good agreement.
Firstly, the computational domains were designed in order to be representative of the experimental samples used in the pressure drop measurements. The computational sample thickness in the flow direction was kept as close as possible to the corresponding experimental one. Additionally, a sensibility study was carried out in order to find an optimal perpendicular cross-sectional size for the computational domains, so that the pore-scale flow results are not affected. Pressure gradient and Darcian velocity values were computed from the porescale simulations by averaging and compared against experimental data. Permeability and Forchheimer coefficient values were computed from second-order polynomial fits applied to the pressure gradient curves for all cases, assuming the pressure gradient to vary according to the Forchheimer equation.
Comparison of permeability values obtained from experiments and pore-scale simulations showed large discrepancies in the results. This is explained by the fact that the permeability linear term in the Forchheimer equation represents the viscous losses. Since the flow regime investigated here is predominantly in the turbulent regime, viscous losses should be negligible. Additionally, the curve fitting procedure searches a 2-D parameter space, where there might be more than one pair of permeability and Forchheimer coefficient values as a solution.
Therefore, under the flow conditions investigated here, it is reasonable to assume that the pressure losses are caused by inertial effects only, represented by the quadratic term in the Forchheimer equation. Hence, another curve fitting procedure is performed assuming the pressure gradient to vary according to Eq. 28 (which assumes the linear term in the Forchheimer equation to be zero). Comparison between experimental and pore-scale simulations Forchheimer coefficient values calculated with Eq. 28 showed relatively good agreement. These results also support the feasibility of using RANS turbulence models to account for the turbulent within the pore space. Better agreement with experimental data was achieved when the RNG k-ε model was used. Nevertheless, a more detailed study is necessary to evaluate the reasons for this better agreement. Results have shown that permeability is inversely proportional to the specific surface area of the foams, whereas the opposite is true for the Forchheimer coefficient.
A source of deviation between the experiments and simulations could have been caused by a certain degree of anisotropy that is present in the porous matrix of the foams. As shown in the work of Ranut et al. (2014), open-cell metal foams normally have the pores and struts elongated along one particular direction due to their manufacturing process. Simulations results obtained by Ranut et al. (2014) have shown that changing the direction of the incoming flow can cause deviations of up to 20 % on the pressure gradient. Unfortunately, anisotropy in the foam structure and its effect on the flow were not evaluated in the present work, but could explain some of the discrepancies found.
Perhaps one of the most important conclusions from the current work is the fact that the sample size-related effects on the flow were observed by investigating samples of varying thickness in the flow direction. This points out to the fact that the REV size determined solely by morphological properties, such as porosity, might not be sufficient for the REV to provide flow-dependent results. Thus, a flow-dependent REV size has to be determined based on the flow characteristics, which can in turn be determined either experimentally or numerically using pore-scale simulations. Results for the Inconel 450-µm samples with varying thickness have shown the critical thickness to be between 25 and 50 d p , in agreement with experiments by Dukhan and Patel (2010). The critical thickness was defined as the thickness for which pressure gradient results converge asymptotically towards a constant value.
property of the E-BREAK Consortium, and it must not be diffused without the prior written approval of the authors.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.