1 Introduction

Open-cell foams which are also known as solid sponges benefit from advantageous morphological properties such as a continuous solid network, the possibility for lateral flow, high porosity and large specific surface area (Zhu et al. 2018). These structural features make open-cell foams superior compared to other structures such as honeycombs and pellets in terms of low pressure drop, high thermal conductivity and excellent solid–gas heat and mass transfer. This has attracted the attention of many investigations in a variety of applications since the last decade (Smorygo et al. 2017). For example, open-cell foams have been investigated by experimental and numerical studies for utilization as heat exchangers (Sertkaya et al. 2012; Kuruneru et al. 2018), heat sinks (Al-Athel et al. 2017), in acoustic and thermo-acoustic applications (Álvarez-Láinez et al. 2014; Napolitano et al. 2017), as volumetric solar energy receivers (Du et al. 2018; Avila-Marin et al. 2018) and as monolithic catalyst supports (Frey et al. 2016; Della Torre et al. 2016; Ulpts et al. 2016, 2017). In industry and literature, open-cell foams are classified based on their pore density expressed as the number of pores per inches (PPI number). These porous structures contain an interconnected solid network of struts which are connected via pore windows. Porosity, pore, window and strut diameter, shown in Fig. 1, are the most important properties influencing the performance of the open-cell foams. A number of studies experimentally investigated the impact of these morphological characteristics on fluid flow parameters and thermal performance, such as pressure loss (Kumar and Topin 2017; Della Torre et al. 2014; Inayat et al. 2016), effective thermal conductivity (Bracconi et al. 2018; Ranut et al. 2015) and heat transfer coefficient (Xia et al. 2017; Dietrich 2013).

Fig. 1
figure 1

Microscopic image of a 20 PPI Al2O3 foam indicating the area-equivalent strut diameter (ds), window diameter (dw) and a single-pore diameter (dp)

Full size image

1.1 Computational fluid dynamics

Since the geometries of many open-cell foams are irregular and complex, in a large fraction of computational fluid dynamics (CFD) simulation studies, idealized or simplified geometries of the samples such as Kelvin cells and Laguerre–Voronoi unit cell models are used to reduce computational time and difficulties with the generation of the computational network (Kumar and Topin 2017; Lucci et al. 2017; Nie et al. 2017; Pabst et al. 2018). Recently, several studies utilized micro-computed tomography (µCT) imaging techniques and subsequent image processing to obtain real digital structures of samples (Meinicke et al. 2016; Zafari et al. 2015; Fan et al. 2016; Ranut et al. 2014). Since meshing the entire and realistic structure of open-cell foams is a demanding procedure and in order to reduce computational time, the majority of µCT-based and idealized geometry-based CFD studies restricted the calculation to a representative elementary volume (REV) of the porous structure.

In a CFD simulation study, Lucci et al. (2017) compared the real geometry of open-cell foams from reconstructed µCT images to their equivalent idealized geometries obtained from the Kelvin cell model. The modeled Kelvin cell structures matched the µCT scans in terms of porosity and specific surface area. The generated ideal geometry predicted a 15–30% larger pore diameter, lower pressure drop and a higher mass transfer coefficient compared to the real structure of the sponge. Della Torre et al. (2014) analyzed the effect of morphological parameters on pressure drop along the ideal and real microstructure of the foams in different flow regimes. The results showed that for different flow regimes, pressure loss along the foam increases quadratically with Reynolds number, which reinforces the assumption that a laminar modeling approach can predict pressure drop with adequate precision even at turbulent flow regimes. In the other studies Meinicke et al. (2016) and Jin et al. (2015), it was further mentioned that microscopic turbulence does not occur in this type of porous structures in the low range of Reynolds numbers. In fact, turbulence is limited to the single microscopic level cells. Thus, at the macroscopic level, it is not necessary to consider turbulent models for flow studies at low or medium Reynolds numbers.

Zafari et al. (2015) investigated incompressible laminar flow and heat transfer within highly porous metal foams (85–90% porosity) and reported a more than 70% increase in velocity compared to the inlet velocities of the fluid (velocity at the inlet of pipe) when the flow reached the foam for different Reynolds numbers. Similarly, Fan et al. (2016) examined the heat transfer and incompressible laminar fluid flow within SiC foams using µCT-based finite element CFD simulations and observed a roughly 60% enhancement of fluid velocities within the foams compared to their corresponding inlet velocities for a wide range of Reynolds numbers. This indicates that incompressible flows within sponges experience velocity enhancement depending on the porosity of the structure. Additionally, it shows that the flow pattern relies more on morphological properties than on the inlet velocity. Full-field analysis of transport phenomena within open-cell foams can be obtained from both numerical methods and noninvasive measurement techniques (Ulpts et al. 2016). Among the available methods, CFD simulation is one of the most cost-effective options for the study of transport phenomena through open-cell foams. In addition, the reliability of such approach gives an insight into the performance of the open-cell foams in different applications. In fact, CFD is an essential tool for precisely elucidating transport phenomena which are often arduous to measure experimentally.

1.2 Experimental investigations

Nuclear magnetic resonance (NMR) and microparticle image velocimetry (µPIV) are two well-known noninvasive methods which have been applied for measurement of flow within open-cell foams (Meinicke et al. 2016; Onstad et al. 2011). The majority of published numerical investigations on the fluid flow within foams were validated against experimental data and literature correlations by comparing pressure loss along the sponge. In contrast, such nondestructive measurement techniques provide the opportunity to investigate transport phenomena at the pore scale level and validate the quality of numerical studies against experimentally determined local velocity and temperature profiles.

Meinicke et al. (2016) utilized CFD simulations and µPIV measurements to investigate different liquid flow regimes within REV of transparent SiO2 glass sponges. The obtained velocity fields and profiles from CFD simulations and in situ µPIV measurements were compared. Surprisingly, only the CFD simulations with the least refinement in computational networks (against common practice in CFD simulations) showed sufficient agreement with µPIV measurements and only at low superficial velocities.

Nuclear magnetic resonance velocimetry (MRV) is a powerful toolkit for the precise characterization of gas velocity also in opaque systems which do not provide optical access. MRV gives proper spatially resolved velocity fields allowing comparison with numerical simulations. However, MRV of gases is challenging due to their low density, short transversal relaxation time and high diffusivity. Therefore, special care should be taken to achieve a sufficient signal-to-noise ratio (SNR) in gas imaging. In recent years, some studies have been conducted to measure thermally polarized gases in porous media, in particular in heterogeneous catalysts (Newling et al. 2004; Ramskill et al. 2016; Koptyug et al. 2000, 2001; Koptyung et al. 2009). Nevertheless, the MRV methods needed to be optimized to make it a routine method for a full-field analysis of gas flow in heterogeneous catalysts. Accordingly, in this study an approach is suggested to measure velocity fields of thermally polarized gases in monolithic structures with 3D spatial resolution and appropriate SNR. The proposed MRV method is based on the adjustment of repetition time (TR), excitation flip angle (FA) and echo time (TE) considering the transverse relaxation time (T2) and the diffusion coefficient (D) of the desired gas. Furthermore, RF-phase cycling was used to suppress unwanted coherence pathways and avoid systematic phase and thus velocity errors.

MRV techniques have been applied for flow investigation of staggered herringbone mixers (Wiese et al. 2018), porous filters (Huang et al. 2017), serpentine mixing chips (Paulsen et al. 2010), blood vessels (Edelhoff et al. 2013) and the cylindrical coquette flow (Hanlon et al. 1998). In addition to flow measurements, NMR-based techniques provide the opportunity to measure temperature and species concentration during thermal and catalytic processes (Ulpts et al. 2015, 2016, 2017). Onstad et al. (2011) performed full-field measurements of flow through highly porous open-cell foams by MRV to better understand the flow behavior inside the foams. Although the results differed in the inlet velocity increase within foams (around 30%) compared to reported CFD simulations (Zafari et al. 2015; Fan et al. 2016), similar trends for velocity profiles at different Reynolds numbers in the laminar flow regime were obtained.

Full-field measurement of flow within non-transparent complex geometries such as open-cell foams, which are required to validate simulations of fluid flow motion based on a real and whole foam structure, is considered a demanding task (Della Torre et al. 2014). Furthermore, measurements of gas flow by MRV are more demanding compared to liquids.

1.3 Comparison of simulations and experiments

The combination of numerical models and experimental measurement techniques provides the opportunity to assess the quality of models for predicting transport phenomena within open-cell foams that are based on µCT data. In addition, the vast potential of NMR techniques for operando measurements of temperature and species concentration profiles within open-cell foams gives the opportunity to evaluate CFD simulation results in terms of temperature fields and local concentration maps in the future. The development of verified CFD models considering real catalytic sponges will enable us to investigate and optimize the performance of monolithic fixed-bed reactors. The irregular and random structure of the open-cell foams makes precise matching of obtained velocity fields from MRV and CFD cumbersome. This is difficult in both non-transparent and transparent foams which forced studies in the past to reduce the region of comparison even from REV to just a few pores (Ranut et al. 2015; Meinicke et al. 2016; Zafari et al. 2015). Such an approach, however, has not yet provided a sufficient match between simulation and experimental data. Here, we hypothesize that transport phenomena within open-cell foams are highly dependent on the interconnected local morphology. We therefore think that it is insufficient to select only a few pores of interest for a CFD simulation without considering the neighboring pores.

We thus claim that reliable simulations of incompressible gas flow through open porous foams can be obtained if the entire foam is included in the simulations. For demonstrating this, we investigate gas transport in a tube filled with three 20 PPI Al2O3 open-cell foams arranged in a row using simulations and experiments. The CFD simulations are based on µCT images of all three foams. From simulation data, the pressure loss of the system is examined and compared against correlations available in the literature. Experimental flow mapping data are obtained by MRV techniques. Finally, the obtained velocity fields from CFD and MRV are aligned and compared qualitatively and quantitatively.

2 Materials and methods

2.1 Reconstruction of foam structure

The μCT images were taken using molybdenum as target material, angular resolution of 0.25°, acceleration voltage of 90 kV, cathode current of 170 μA and exposure time of 0.4 s. The binarized 2D black-and-white sequence of CT images (voxel size of roughly 0.0265 mm) was processed via image processing and CAD software to reconstruct the real geometry of foams. For binarization of CT images, a threshold value should be determined to differentiate between solid and fluid regions. Plotting a histogram of the gray values of the entire CT stack gave two very distinct non-overlapping peaks (for the two phases, solid and fluid). Hence, selection of a threshold value was straightforward by selecting the value that lies in the center between both peaks (2750). After image thresholding, the cross sections of the obtained 3D structures were compared against the CT images to examine the validity of the chosen threshold value (see Fig. 2).

Fig. 2
figure 2

CT images (a) and cross sections of the reconstructed foam (b) and corresponding meshing of a selected part of the CT image (c and d)

Full size image

The triangulated structure of the foams had 73 million surface cells and was exported as an STL (stereolithography) file. The quality of triangulated STL file was improved using the meshing software MeshLab (https://www.meshlab.net/). Afterward, the geometry of reactor was generated with a CAD software (FreeCAD, https://www.freecadweb.org/) and the reconstructed foam structure was inserted into the reactor at an axial coordinate of 150 mm. The procedure used to reconstruct the foam geometry is shown in Fig. 3.

Fig. 3
figure 3

Procedure for preparation of foam reactor geometry

Full size image

Several methods were utilized to determine structural features of foams. Porosity, window and pore size distribution were obtained via volume imaging techniques. The porosity was calculated by dividing the open pore space by the total space (open porosity space and strut space) of the foams. In order to calculate the open porosity, a cylindrical volume at the center of foams was considered as REV. The radius of the REV increased, and the averaged open porosity of corresponding REVs was obtained for each radius until the obtained porosities reached a stable value. Afterward, window and pore size distribution were obtained by analyzing the determined REV with the chosen radius. In order to determine the density of the foams, diameter, height and mass of foams were measured using a standard caliper and a laboratory balance, respectively. The specific surface areas of foams were calculated via the ideal geometry modeler of Inayat et al. (2011). Table 1 contains all measured morphological properties of foams.

Table 1 Morphological properties of utilized 20 PPI open-cell foams for simulation and experimental velocimetry measurements
Full size table

The commercial meshing software cfMesh was employed to prepare the CFD grid and to generate a Cartesian mesh consisting of predominantly hexahedral cells with polyhedra in the transition regions between the cells of different sizes in the fluid domain. This meshing software is based on the surface triangulated geometry of the foam reactor which was imported as an STL file.

2.2 CFD simulation model

Flow within porous media can be categorized into four different regimes based on the pore Reynolds numbers (Dybbs and Edwards 1984):

$${\text{Re}}_{p} = \frac{{\rho_{{\text{f}}} u_{\sup } d_{p} }}{{\mu_{{\text{f}}} }}.$$
(1)

Here, usup is the superficial velocity (the volumetric gas flow at the pipe inlet divided by cross-sectional area of the pipe), ρf is the fluid density, \(\mu_{{\text{f}}}\) is the dynamic viscosity of fluid and dp is the pore diameter. For Rep < 1, the flow is controlled by viscous forces and the pressure drop depends linearly on velocity. The Darcy–Forchheimer regime, in which pressure loss is nonlinearly correlated with velocity, occurs in the range of 1–10 < Rep < 150. This regime is dominated by inertial forces, and the pressure drop can be written as

$$- \frac{\partial p}{{\partial z}} = \frac{{\mu_{{\text{f}}} }}{{k_{1} }}u_{\sup } + \frac{{\rho_{{\text{f}}} }}{{k_{2} }}u_{\sup }^{2} .$$
(2)

Here, k1 and k2 are the Darcian and non-Darcian permeabilities, p is the pressure and z is the axial coordinate. At higher Reynolds numbers, 150 < Rep < 300, which is called post-Forchheimer regime, laminar wake wavering occurs and for Rep > 300 fully turbulent flow appears.

In this work, CFD simulations of steady-state flow through three 20 PPI Al2O3 foams arranged in a row are performed for different flow rates in the Darcy–Forchheimer regime. For case studies, the structure of the sponge was adopted to the geometry of an NMR compatible model reactor with 25 mm inner diameter. Since in all simulated case studies, the Mach number is smaller than 0.2, the assumption of incompressible flow is acceptable for simulations. We used methane as probing gas since it is easily available and NMR measurements of the CO2 methanation reaction are planned to be done in near future. The governing steady-state equations of continuity and momentum for incompressible and Newtonian flow can be written as (Versteeg and Malalasekera 2007)

$$\nabla \cdot {\varvec{u}} = 0,$$
(3)
$$\rho_{{\text{f}}} \left( {{\varvec{u}} \cdot \nabla } \right){\varvec{u}} - \mu_{{\text{f}}} \nabla^{2} {\varvec{u}} = \nabla p.$$
(4)

Here, u denotes the velocity vector. To solve Eqs. (3) and (4), the finite volume-based CFD software OpenFOAM 4.0 was utilized. OpenFOAM is based on the finite volume method (FVM) that discretizes the conservation equations on a fixed mesh. The SIMPLEC (semi-implicit method for pressure-linked equations-consistent) algorithm was applied for pressure–velocity coupling, and the simpleFoam solver in the OpenFOAM library was chosen for running laminar steady-state simulations (Versteeg and Malalasekera 2007; Report et al. 2016). The absolute tolerance for solving the system of equations was set to 10–6. For the approximation of the convection cycle, the upwind method was used.

The boundary conditions were as follows: At the inlet, a fixed superficial velocity in the z direction, u = (0, 0, uz), was assumed. At the outlet, a fixed pressure, p = 0, was prescribed. The tube walls and the sponge surface were prescribed with a no-slip boundary condition (u = 0). In order to ensure a fully developed flow before reaching the porous structure, the adequate entrance length Le was calculated and considered in the model by (Spurk and Aksel 2010):

$$L_{{\text{e}}} = 0.056 {\text{Re}} D_{{{\text{in}}}} ,$$
(5)
$${\text{Re}} = \frac{{\rho_{{\text{f}}} u_{{\text{z}}} D_{{{\text{in}}}} }}{{\mu_{{\text{f}}} }}.$$
(6)

Here, Din is the diameter of the reactor and uz is the inlet velocity. CFD simulations were carried out for inlet velocities of 0.017, 0.034 and 0.051 ms−1 (0.5–1.5 standard liter per minute, SL min−1) which belong to the Darcy–Forchheimer in laminar flow regime. The calculated Reynolds number based on the diameter of the reactor (25 mm) was 76.49 for the inlet velocity of 0.051 m s−1, and its corresponding pore Reynolds number (Rep) based on the average pore diameter is 10.44. Physical properties of methane as operating fluid were calculated at a temperature and pressure of 293 K and 1 bar. The mass density ρf and dynamic viscosity µf were 0.6596 kg m−3 and 10.9945 × 10–6 Pa s, respectively (Touloukian and Vargaftik 2014).

To ensure that the simulations are grid independent, a grid analysis was carried out for different meshing densities, which differed in their cell numbers. For this purpose, the pressure loss caused by the foam’s structure at an initial velocity of 0.051 ms−1 (1.5 SL min−1) was compared for the different grids to determine the critical number of cells at which the pressure loss is independent of the grid. A mesh containing 41.1 million cells was determined to be sufficient to achieve grid independence (~ 10 million cells per sponge). The results of the grid independency test is shown in the Supporting Documents (see Figure S1). To better resolve the large gradients occurring in the foam’s structure, a mesh refinement was applied which started from 10 mm before the sponges and continued until 10 mm after the sponges. Simulations were parallelized on 12 processors (Intel® Xeon(R) E5-2650 v4 2.20 GHz CPU speed), and each simulation took around 12 h to converge for the chosen mesh density (41.1 million cells).

In order to assess the advantage of using the demanding µCT-based CFD simulations, which consider the real structure of solid sponge, we compared µCT-based CFD simulations with homogenous (empirical) calculations based on the Darcy–Forchheimer pressure drop equation (Eq. (2)). The latter represents the effect of the foam on the flow during calculations. By fitting of the Darcy–Forchheimer equation to the obtained pressure losses from µCT-based CFD simulations for different superficial velocities, it was possible to determine the corresponding values of k1 and k2 which are 8.5413 × 10–8 m2 and 1.7301 × 10–3 m for 20 PPI sponges, respectively. This means that both µCT-based and homogeneous models are finally based on µCT images. Practically, the homogenous CFD simulations were performed using porousSimpleFoam within OpenFOAM’s CFD package and the tool blockMesh was used for the creation of the computational network.

2.3 MRV measurements

A spin echo (SE) and phase contrast (PC)-based MRV pulse sequence was implemented. In the employed non-slice-selective 3D MRI sequence, spatial resolution was achieved by one read and two phase encoding gradients. Velocity encoding was performed using a four-step Hadamard scheme.

All measurements were performed on a horizontal 7 T scanner (BioSpec 70/20 USR, Bruker BioSpin MRI, Ettlingen, Germany), which is equipped with a 114-mm bore gradient system (B-GA 12S2), maximum gradient strength of 440 mT m−1 and a maximum slew rate of 3440 mT m−1 ms−1. ParaVision 5.1 was used as the user interface to acquire the measurements. The post-processing of the data was performed using a self-developed MATLAB script (R2014a MathWorks, Natick, MA, USA). A horizontal 72-mm bore birdcage 1H quadrature transceiver RF coil (Bruker BioSpin MRI, Ettlingen, Germany) was used in all measurements.

A glass reactor was used (ID 25 mm, length 80 cm) for the measurements. The utilized heat shrink tubing helped to preserve the samples in the same position of µCT images and also to hold them tight in the reactor. The dimension of the monolithic cylinder after attaching the sponges was as follows: diameter 25 mm and length 68 mm. Methane flow was regulated by mass flow controllers (F-201CV, Bronkhorst High-Tech, Ruurlo, The Netherlands; FMA-2618-A, Omega Engineering, Norwalk, CT, USA). The experimental setup is illustrated in the Supporting Documents (see Figure S2). The gas reactor was placed in the RF coil. Prior to each measurement, the reactor was purged. All the measurements were performed at ambient temperature and pressure of 1 bar. Various volumetric flow rates were applied, i.e., 0.5, 1 and 1.5 SL min−1. The highest obtained velocity in NMR was predicted to be 0.3 m s−1 according to the theory for a flow of 1.5 SL min−1. The parameters of the PC MRV pulse sequence were optimized as described in detail in Huang (2017), primarily to avoid SNR losses that are caused by relaxation and diffusion processes. The measurement parameters were as follows in PC-MRV. Field of view (FOV) was set to 96 × 64 × 64 mm3 in a matrix 120 × 80 × 80 giving a resolution of 0.8 × 0.8 × 0.8 mm3. According to the Hadamard velocity encoding approach, for each volumetric flow rate, four measurement steps were applied. A velocity encoding range [− VENC, + VENC] with VENC values = 300, 250 and 150 mm s−1 was used according to the expected maximal velocity for flow rates of 1.5, 1 and 0.5 SL min−1, respectively. Flow encoding duration and flow encoding delay were 0.37 ms and 1.69 ms, respectively. The measurements were averaged 8 or 16 times for each flow rate to improve SNR. Using a FA of 127°, TE of 3.15 ms and TR of 50 ms, the total measurement time for measurement was about three hours. The data were processed using homemade MATLAB code (R2017a MathWorks, Natick, MA, USA). Binary mask was used based on Otsu thresholding algorithm after the reconstruction corresponding to the four steps in the utilized Hadamard. Pixels included in the mask were considered as methane. Commercial open-cell foams were provided by Hoffmann ceramic GmbH (Breitscheid, Germany). Three 20 PPI Al2O3 foams were utilized to examine flow within porous media. In order to fix the position of foams, they were wrapped with a polyethylene heat shrink tubing with meltable adhesive liner from Conrad Electronic SE Company in the axial direction. This helps to assure μCT images, and NMR data are performed at the same orientation for each pore of foam.

3 Results and discussion

3.1 CFD model validation

To guarantee that the simulated flow results are physically correct, they were validated against pressure loss correlations for foams at flow rates from 0.5 to 20 SL min−1 (Fig. 4). Three correlations of Dietrich et al. (2009), Lacroix et al. (2007) and Fourie et al. (2002) were employed. Since most published experimental data on pressure loss were acquired at porosities of above 0.9 and air as operating fluid, no experimental data sets with approximately equal sponge properties could be found. By inserting morphological properties of sponges such as porosity ε, pore diameter dP, window diameter dw and specific surface area SV into the correlations, a pressure loss estimate was obtained. The CFD simulation results show good agreement with the literature correlations, especially in the range of the inlet velocities applied in this study (below 0.1 m s−1).

Fig. 4
figure 4

Comparison of pressure loss obtained by CFD simulation (symbols) foams with literature correlations data (lines)

Full size image

Average axial and radial velocity profiles were obtained from the CFD simulations. In order to obtain average radial velocity profiles of the axially directed flow at different positions, several planes perpendicular to the direction of flow were considered at different axial positions z. These planes were segmented into several rings with equal area; an average velocity can be calculated for each individual ring segment (Fig. 5). Here, 100 equal area segments were used for the following analysis. The reason for choosing rings with equal area instead of the rings with equal diameter (or thickness) is as follows: In the case of equal-diameter rings, the inner rings contain only few cells, so the average is quite sensitive to local high or low porosities. Also, average axial velocity profiles were obtained by averaging the velocity over the whole plane area at different axial positions z.

Fig. 5
figure 5

Schematic view of the adopted segmentation method for calculating radial velocity profiles. Each section has equal area

Full size image

3.2 Comparison of micro-CT-based and homogenous models

We compared the average axial velocity profiles obtained from heterogeneous CFD simulations with profiles obtained using the homogeneous models for different flow rates from 1.5 to 20 SL min−1 (Figure S3 in Supporting Documents). Both the µCT-based and the homogeneous model show the expected increase in the velocity in the porous section of reactor. It should be noted that throughout the manuscript, we plot averaged velocity profiles (Figs. 8, 10) that show the axial component of the velocity (uz) and thus do not represent the magnitude of velocity.

To calculate the averaged velocity profile over whole slices or rings, the integral of the velocity’s axial component (uz) over the corresponding cross-sectional area was calculated using ParaView (integrate variable filter). The average value of the axial component of the velocity can then be obtained by dividing the obtained flow rate from ParaView by the corresponding cross-sectional area. In case of CFD simulations, the utilized cross-sectional area contains only the fluid part or open area of the foam which was obtained from CT images.

As expected, the local fluctuations in the velocity profiles due to the local heterogeneity cannot be observed in the homogenous model. As transport phenomena within solid sponges are highly dependent on local morphological characteristics, µCT-based models naturally contain vastly more information on flow and thus can be advantageous in deducing design guidelines. In addition to the velocity increase at the beginning of the porous section of the foam, we can observe local fluctuations in the two gaps between the middle sponge and its surrounding sponges. These fluctuations of velocity are due to the narrow non-porous region in the gap. Due to the two gaps between each two adjacent sponges, the velocity shows local minima when the flow passes from one sponge to the next.

Average radial profiles of uz (as obtained from the local averaging method described in Fig. 5) within the foam obtained from CFD were also compared at different axial locations (Fig. 6). The figure additionally shows the velocity profile obtained from the homogeneous model (lines) right in front of the foam (− 3 mm) and inside the foam. The three panels in the figure represent profiles of the axial velocity as a function of the normalized radial position of the reactor for the first, middle and the last sponge, respectively. Right in front of the sponge, the flow pattern changes from a parabolic laminar flow profile due to the sponge’s porous structure. The boundary layer thickness of flow through the sponge is easily identified from the plots. Since porous media properties for all three sponges are the same, the obtained average radial velocity profiles from the homogeneous model are identical at all axial locations inside the sponge. The red markers represent the radial profile of uz for 3 mm before the beginning of porous structure. In all three subsections, the blue markers belong to radial profile of uz of 3 mm after the beginning of each sponge, the green markers represent radial profile of uz in the center of each sponge and the gray ones show radial profiles of uz in 3 mm before the end of each sponge. Fluctuations of the velocity as a function of r are due to the local porosity and window size. Due to the irregularity of the structure, every plane perpendicular to the flow direction at each axial location has its unique radial velocity profile. Furthermore, every ring in each plane has its unique velocity value which is related to the local porosities. The top section of Fig. 6 illustrates the considered positions of the slices in the µCT-based model.

Fig. 6
figure 6

The average radial profiles of uz for the first sponge (left), second (middle) and third sponge (right) from the homogeneous (dashed lines) and µCT-based models (markers): 3 mm before the beginning of the first sponge (red markers), 3 mm after the beginning of each sponge (blue markers), at the center of each sponge (green markers), and 3 mm before the end of each sponge (gray markers). For all plots, the flow rate was 1.5 SL min−1. The boundary layer is indicated by the shaded region

Full size image

The obtained velocity profiles from µCT-based and homogeneous models for the axial position of 3 mm in front of the first sponges in Fig. 6 (red markers titled with − 3 mm in the left panel) showed small deviation. This deviation explains the effect of the sponges on the flow filed within pipe and before them which cannot be predicted by homogenous model. The difference between homogeneous and µCT-based models increased for the axial positions which are located within sponges. Furthermore, average radial profiles of uz of homogeneous and µCT-based models showed higher deviation in comparison with average axial profiles of uz (Figure S3 in Supporting Documents). As it is mentioned in Fig. 5, the averaging area for radial profiles of uz (rings with the same area) is smaller than the ones for average axial velocity profiles (the whole cross section). This indicates that homogenous results have better agreement with µCT-based results when they are averaged on the larger region of interest. Nevertheless, the effect of local morphology cannot be distinguished with the homogenous model in terms of average radial and axial profiles of the axial velocity uz.

3.3 Comparison of CFD and MRV

The velocity fields, obtained via CFD at different flow rates, were compared to MRV measurements. The computational cells and MRV voxels have an edge length of 0.1 and 0.8 mm, respectively. Besides the advantages over homogeneous models, this resolution gap between CFD results and MRV data further underlines the demand for comparing velocity pattern of different methods. Using grid networks with less refinement and thus a computational cell size that matches the MRV voxel size would, according to the mesh independency test, not deliver accurate results as the mesh is too coarse for an accurate solution of the momentum equations. Furthermore, determination of the exact assessment position and orientation alignment of CFD and MRV data is a demanding task. It should be noticed that experimental measurement errors are possible during MRV measurements. The flow behavior in foams is highly dependent on the extremely complex morphological characteristics. Therefore, as a first step, it is essential to prove the capability of MRV to measure gas flow and assess its quality for a simple system. For this, we investigated flow of methane inside an empty reactor (ID 25 mm) at flow rates of 0.5 and 1.5 SL min−1. The CFD and MRV results were compared in terms of average axial velocity; both methods agree reasonably well with a maximum deviation of ± 10% (see Supporting Documents, Figure S5). This justifies the ability to use MRV as a reliable tool to assess CFD simulations and also for studying gas flow within porous structures.

Three different flow rates (0.5, 1.0 and 1.5 SL min−1) were chosen to investigate the gas flow of methane within the three 20 PPI Al2O3 sponges. Figure 7 shows the aligned (the alignment method is demonstrated in Fig. 12) velocity fields for the flow rate of 1.5 SL min−1 from CFD simulation (top row) and MRV measurement (bottom row) for random cross-sectional slices at different axial positions of the middle foam. By comparing the obtained velocity fields, similar flow pattern can be observed and distinguished from both methods. Naturally, in CFD simulations only the fluid part of the foam is considered as the solid part of the sponge is not included in the simulation. Therefore, the net cross-sectional area of fluid can be calculated at each axial position. In the case of the MRV fields, the maxima and minima of the velocity are less pronounced compared to the CFD velocity fields. Also, the non-fluid part of each slice has a nonzero value of velocity as many voxels at the fluid–solid interface belong to both regions. Consequently, for quantitative comparisons, a threshold filter was applied to obtain the open area of each cross-sectional slice. To find a suitable threshold value, the averaged axial velocity profile and the resulting porosity over the corresponding surface were calculated and compared with the CFD values to find the best match, which resulted in a velocity threshold value of 15 percent of the inlet velocity. Figure S4 in the Supporting Documents illustrates the comparison of open porosity (open area divided by the area of whole cross section) from MRV and CFD measurements for a flow rate of 1.5 SL min−1. It can be concluded that the obtained gas velocity fields from MRV measurements cannot predict the foam’s open porosity as precise as µCT imaging. However, MRV offers the opportunity to investigate velocity patterns of gas flow within complex opaque structures which cannot be offered by other noninvasive techniques.

Fig. 7
figure 7

Simulated and measured velocity fields of gas flow (1.5 SL min−1) within the foam reactor via CFD and MRV measurements within radial cross sections in different axial positions

Full size image

Figure 8a shows the average axial velocity profiles within the porous section of the reactor obtained via CFD and MRV. The results show remarkable agreement with an accuracy of ± 15%. The middle porous section of reactor (middle sponge) shows a better agreement in both the overall trend and the local values. Further, in the entrance and terminal regions of the porous structure (the first and third sponge), a considerable gap between CFD and velocimetry results was observed. The 68 mm total length of the utilized sample causes increasing deviations at the edges of the field of view. In fact, these regions are far away from magnet iso-center in the NMR machine. For the axial position far from the iso-center of the magnet, lower SNR and the effect of nonlinearity in gradients can be observed. Thus, the uncertainty increases with increasing distance from the center which originates from MRV measurements as the CFD simulations use µCT images with constant resolution and mesh density. The uncertainty of MRV data is fully discussed in the Supporting Documents (section S2).

Fig. 8
figure 8

Comparison of average axial velocity profile of methane within a foam reactor packed with a three 20 PPI Al2O3 sponges in a row and for different flow rates of 0.5, 1.0 and 1.5 SL min−1b within the middle sponge. The empty and solid markers represent the obtained results from CFD simulations and velocimetry measurements, respectively

Full size image

Considering only the middle sponge (Fig. 8b), the deviation decreases to ± 10%. Investigating only the middle sponge separately is justified as the three sponges are loosely stacked and then taped. The µCT measurements are performed in a different building than the MRV measurements. Thus, it cannot be guaranteed that the sponges did not move slightly relative to each other. We thus only aligned the middle segments for comparison. The quantitative agreement between CFD and MRV results underlines the reliability of CFD simulations. This also underpins the capability of NMR-based measurements for in situ temperature and concentration measurements since NMR temperature and concentration mapping data showed good agreement with 2D pseudo-heterogeneous modeling results in previous studies (Ulpts et al. 2017).

Figure 9 shows parity plots of the CFD and MRV comparison in terms of average axial velocity profiles. These comparisons verify the agreement between simulation and experimental results with ± 10% accuracy in major parts of the reactor. By limiting the field of view from three sponges (Fig. 9a–c) to that sponge which is nearest to the magnet iso-center (Fig. 9d–f), these deviations decreased significantly. The points outside of the ± 15% gray bands in Fig. 9a–c belong to the inlet and terminal regions of the foams.

Fig. 9
figure 9

The comparison of CFD and MRV results in terms of average axial velocities for three different flow rates of 0.5 (top row), 1 (middle row) and 1.5 (bottom row) SL min−1 and for different field of views of three sponges in row (left panels) and middle sponge only (right panels). The gray bands indicate a deviation of ± 15% (left panels) and ± 10% (right panels)

Full size image

Table 2 shows the linear equation coefficients (y = mx + b) of the best-fit lines in Fig. 9 and the respective R2 values. The good quality of agreement confirms the importance of considering the whole structure of the ceramic foams in simulations. In order to prove this hypothesis, new simulations were carried out on a REV. For this, we simulated gas with a flow rate of 1 and 1.5 SL min−1 in a restricted part of the middle sponge with a diameter of 8 mm. The adopted boundary conditions for the inlet and outlet were the same as for the full system described in Sect. 2.2. A symmetry boundary condition was chosen for the wall, and the number of computational cells was 11.2 million. Simulation results were compared with MRV data of their equivalent region of interest in terms of average axial velocity profile. Results showed noticeable deviation caused by neglecting neighboring pores of the REV in simulation. Figures S6 and S7 in the Supporting Documents illustrate the comparison of average axial velocity profiles and parity plots of CFD and MRV data from the REV of the foam. In addition to the different trends of the obtained average axial velocity profiles from CFD and MRV, the results show a higher deviation (within ± 20%) in comparison with the full-field system.

Table 2 Fit parameters and R2 values of parity plots (Fig. 9)
Full size table

In addition to axial velocity profiles, CFD and MRV data were compared in terms of radial profiles of the axial velocity, uz (as obtained from the equal-area segmentation method, Fig. 5). We considered 20 rings with equal area over the slices of the velocity fields obtained from MRV and CFD data. CFD and MRV radial profiles for different flow rate and at different sites in the direction of flow are shown in Fig. 10. As it was mentioned before, the MRV data have a lower resolution compared to CFD data (larger voxel size compared to mesh cell size). It means that the CFD data are more sensitive to changes in the local morphology. Consequently, the obtained radial velocity profiles from the CFD simulation show more distinctive fluctuations. In other words, maxima and minima of velocity are much more pronounced in the CFD data. This is also obvious in the obtained velocity field from CFD and MRV (Fig. 7).

Fig. 10
figure 10

Average radial velocity profiles from CFD and MRV data for different flow rates of 0.5, 1.0 and 1.5 SL min−1 and for different positions in the direction of flow

Full size image

In order to quantify the illustrated qualitative agreement between CFD and MRV in Fig. 7, voxel-wise comparison of the velocity fields was carried out. The velocity fields from MRV measurements are structured and voxelized. Due to non-uniform distribution of points in the simulation mesh, the velocity fields from CFD data are unstructured and irregular, but they have a ca. 8 times higher resolution compared to the MRV data. In order to make these two data sets comparable, the CFD data were mapped on the regular grid given by the MRV results. Afterward, the structured CFD data were convolved to the same resolution as the MRV data. In order to fix axial and angular deviation between CFD and MRV data, each MRV velocity field was compared to the CFD velocity fields, while it was being rotated from 0 to 360 degree. In each step, the structural similarity (similarity index) between the MRV slice and the CFD slice was calculated. A clear maximum in the similarity index as a function of rotational angle indicates that the correct rotational angle has been found. Figure 11 shows the voxelized CFD and MRV data (for the same slices as in Fig. 7). Here, again, the top velocity fields in each set belong to the CFD data and the bottom ones belong to the MRV data. Figure 12a shows, as an example, the variation of the similarity index against the rotational angle of the MRV velocity field for slice 4 (z = 39.9 mm). A clear maximum at 70° indicates the ideal rotational angle which corresponds to a similarity index of around 0.7. (The similarity index as a function of the rotational angle for all five slices can be found in Figure S8 in the Supporting Documents.) Figure 12b shows the ideal rotational angles for all z that belong to the middle sponge as well as the corresponding similarity indices. (The hollow square marker indicates slice 4.) The highest possible similarity index is unity for identical images. The averaged angular deviation is 70.86 degree and fluctuates only by three degrees. The similarity index is ca. 0.7 for all slices.

Fig. 11
figure 11

The voxelized CFD (top) and MRV (bottom) velocity fields corresponding to Fig. 7

Full size image
Fig. 12
figure 12

Similarity index vs. rotational angle for slice 4 at z = 39.9 mm (a) and optimum rotation angle (as extracted from the maximum) and similarity index of MRV velocity fields against their corresponding CFD velocity fields (b)

Full size image

4 Conclusions

In this work, gas flow within open-cell foams was investigated numerically by CFD simulation, and the results were compared to MRV measurements. In a first step of this investigation of gas flow pattern in ceramic open-cell foams, mesh-independent CFD calculations were validated against conventional pressure drop correlations. Further, in order to demonstrate the predictive power of the calculated results that are based on µCT data, we compared CFD results with MRV velocimetry measurement data of gas flow with respect to axial and radial profiles of the velocity in axial direction. It became evident that such a comparison requires spatially resolved information about the original morphology of the foam. Homogeneous simulations, however, were found to be insufficient to calculate local and average fluctuations of velocities. The matching showed a good agreement with deviations below 10 percent in most cases. Further, the obtained velocity fields from CFD and MRV showed good qualitative and also quantitative agreement when they were compared voxel by voxel. These results underpin the reliability of both methods. Such a good agreement, however, could not be obtained if only a few pores were considered in CFD simulations (as shown by comparing the MRV data against an REV simulation). It can be concluded that the whole structure of open-cell foams utilized in experimental measurements should be considered also in CFD simulations. Implications of this finding can be expected for CFD-based simulations of other quantities within open-cell foams such as temperature fields and concentration maps in monolithic catalysts.

Our findings reveal that both methods, MRV and CFD simulations, can be used for providing reliable gas flow velocity patterns within non-transparent open-cell foams. They also underpin that the reliability of such CFD simulations reduces substantially if only a section of the foam is considered in the calculations.