Morphological Characterization and Effective Thermal Conductivity of Dual-Scale Reticulated Porous Structures

Reticulated porous ceramic (RPC) made of ceria are promising structures used in solar thermochemical redox cycles for splitting CO2 and H2O. They feature dual-scale porosity with mm-size pores for effective radiative heat transfer during reduction and µm-size pores within its struts for enhanced kinetics during oxidation. In this work, the detailed 3D digital representation of the complex dual-scale RPC is obtained using synchrotron submicrometer tomography and X-ray microtomography. Total and open porosity, pore size distribution, mean pore diameter, and specific surface area are extracted from the computer tomography (CT) scans. The 3D digital geometry is then applied in direct pore level simulations (DPLS) of Fourier’s law within the solid and the fluid phases for the accurate determination of the effective thermal conductivity at each porosity scale and combined, and for fluid-to-solid thermal conductivity from 10−5 to 1. Results are compared to predictions by analytical models for structures with a wide range of porosities 0.09–0.9 in both the strut’s µm-scale and bulk’s mm-scale. The morphological properties and effective thermal conductivity determined in this work serve as an input to volume-averaged models for the design and optimization of solar chemical reactors.


Introduction
Foam-type reticulated porous ceramics (RPC) structures are applied in a broad range of physical processes requiring enhanced heat and mass transfer [1,2]. Applications include microelectronics cooling [3], soil dynamics [4,5], catalytic reactors [6], radiant burners [7], tissue engineering [3,8] and volumetric heat exchangers for the conversion of concentrated solar energy [9][10][11]. Of special interest of the latter application is the solar-driven thermochemical redox cycle for splitting CO2 and H2O [12][13][14][15][16], consisting of: (1) a high-temperature endothermic reduction, in which a metal oxide is thermally reduced and oxygen is evolved; and (2) a lower-temperature exothermic oxidation, in which the reduced oxide is re-oxidized with H2O and CO2 to form H2 and CO (syngas), and further processed to liquid hydrocarbon fuels. Ceria-based oxides have emerged as highly attractive redox materials because of the rapid oxygen transport in the bulk [14][15][16][17][18]. Various porous structures made of ceria have been investigated for enhanced reaction rates [18][19][20], including structures with submicron-sized interconnected pores, but these are problematic to retain because of partial sintering at elevated temperatures [19]. Furthermore, their high optical thickness inhibits penetration of concentrated solar radiation, resulting in non-uniform heating and temperature distributions [14]. Most recently, Furler et al. [15] presented a unique and morphologically stable RPC structure featuring dual-scale porosity: mm-size pores with struts containing micron-size pores. The mm-size pores enable volumetric absorption of concentrated solar radiation and effective heat transfer during the reduction step, while the micron-size pores within the struts offer increased specific surface area leading to enhanced reaction kinetics during the oxidation step.
Optimization of solar reactors for thermochemical redox cycles requires computational models of heat transfer and fluid dynamics coupled to the reaction kinetics [14,21]. Since resolving the solar reactor at the pore scale would require tremendous computational demand, volume-averaging theory is often applied for solving the mass, energy, and momentum conservation equations using effective heat and mass transport properties [22][23][24][25]. These can be determined accurately by direct pore-level simulations (DPLS) using the detailed 3D digital geometry of the structure obtained by computer tomography (CT) [26][27][28]. For example, the Monte Carlo ray-tracing method has been applied at the pore level for solving the radiative heat transfer equations and determining the effective extinction coefficient and scattering phase function [29], and the finite volume (FV) technique has been applied at the pore level for solving the Navier-Stokes equations and determining the effective thermal conductivity, permeability, and heat transfer coefficient [30,31].
In this work we apply the tomography-based methodology to investigate RPC structures made of ceria with dual-scale porosity in the mm and µm scales. This structure is schematically depicted in Figure 1 [15]. The total and open porosity, pore size distribution, mean pore diameter and specific surface area are extracted from the CT-scans. The effective thermal conductivity is determined by DPLS for the RPC with non-porous struts and for the RPC with dual-scale porosity. We investigate the effect that the dual-scale porosity has on the morphological properties and on the conduction heat transfer across the RPC, and further compare the results to predictions by analytical models for structures with a wide range of porosities in both the strut's µm-scale and bulk's mm-scale.

Figure 1.
Ceria RPC with dual-scale porosity: mm-size pores for volumetric radiative absorption and effective heat transfer (a) during the reduction step, and struts containing micron-sized pores leading to increased specific surface area (b) for enhanced reaction kinetics during the oxidation step.

RPC Synthesis
The dual-scale RPC structure is manufactured using the Schwartzwalder foam replication method [32]. An organic foam template is coated with multiple slurry layers containing ceria particles and micron-sized carbon grains [15]. The carbon pore former content ranges from 10 to 50 vol%. After firing at high temperatures (>1800 K), the bearing carbon foam and grains are burned und the desired foam-type structure undergoes sintering.

Synchrotron Submicrometer Tomography
Strut samples are scanned using synchrotron submicrometer tomography with a voxel (3D pixel) size of vs = 325 nm and a 0.832 × 0.832 × 0.702 mm 3 field of view. The high-resolution CT is performed at the Swiss Light Source (SLS) of the Paul Scherrer Institute (PSI, Villigen, Switzerland) with the TOMCAT beamline for 40 keV photon energy, 400 µA beam current, a 100 µm thick aluminium filter, 40 µm thick copper filter, a 10 µm thick iron filter, 20× geometrical magnification, 1 s exposure time, 1001 projections. Figure 2 shows exemplary tomograms of strut samples manufactured with various concentrations of pore former ranging from 10 to 50 vol% and their corresponding 3D reconstructions of the pore space within isotropic strut regions. One tomogram contains 2560 × 2560 pixels. Numerous strut samples are scanned to verify reproducibility.  (10, 20, 30, and 50 vol%) and their corresponding 3D digital reconstruction of the void space within isotropic porous strut region.

Micrometer Tomography
A ceria RPC sample with 10 pores per inch (ppi) is scanned by micrometer tomography with a voxel size of vs = 35.7 µm and a 36.56 × 36.56 × 36.56 mm 3 field of view. The low-resolution CT is performed with an unfiltered polychromatic X-ray beam at the Swiss Federal Laboratories for Materials Science and Technology (EMPA, Dübedorf, Switzerland) for 150 keV photon energy, 45 µA beam current, 6.272 s exposure time and 721 projections. Figure 3 shows a tomogram and its corresponding 3D digital reconstruction of the scanned RPC sample. One tomogram contains 1024 × 1024 pixels.

Porosity
For the analysis of the strut structures, isotropic regions within submicrometer tomograms are cropped with the size of 501 × 501 × 501 voxels (0.163 × 0.163 × 0.163 mm 3 ). For the analysis of the RPC structures, micrometer tomograms are cropped to 500 × 500 × 500 voxels (17.85 × 17.85 × 17.85 mm 3 ). The cropped 8 bit tomograms obtained from measurements are pre-processed with a 3D Gaussian blurring filter to remove unwanted image noise derived from the photon sensor. Histograms computed from 3D tomogram stacks show bimodal character representing two grey scale pixel classes, with the threshold found by Otsu's method of intra-class variance minimization [33,34]. Finally, based on the threshold, each pixel is assigned to be either void or solid. Porosity, ε = Vf/V, is defined as the ratio between the void space volume and the total cube volume. It is calculated by counting void and solid voxels of the 3D stack. The representative elementary volume (REV) defines the minimum volume containing a porous zone for which the continuum assumption is valid. It is determined from incrementally growing cubic subvolumes until their calculated porosities convergences within a certain band, ±γ. The conditions for the minimum edge length of the REV are [27]: where VL* is the sample subvolume and L* is the edge length of the sample subvolume. For the RPC with a porosity band of γ = 0.05, lREV ≥ 6.6 mm, which leads to cube structures larger than 186 voxels edge length. For the µm-sized struts with a porosity band of γ = 0.05, lREV ≥ 76.1 µm , which leads to cube structures larger than 235 voxels edge length. Dual-scale porosity, εdual, is calculated from strut-scale porosity, εstrut, and RPC-scale porosity, εRPC, as: (2) εstrut is linearly fitted to the pore former concentration, ϕ, as εstrut = 0.008707 ϕ. RPC structures with varying strut thicknesses are generated by altering the original segmented tomography scans through a dilation process with 3D spherical elements of a certain diameter, d, with Table 1 lists the dilation radius, digital porosity, and mean pore diameter of the original RPC reconstruction and of the digitally altered RPC for 3 increasing strut thicknesses, and the corresponding digital section cut and 3D rendering. As expected, porosity and mean pore diameter of the RPC decrease with increasing strut dilation because the thicker struts consume void space.
Of special interest is the connectivity of the µm-sized pores within the struts. The pore connectivity scales directly with the specific surface area reachable by reacting gases, and thus scales with the fuel production rates [15]. Open porosity, εopen, is defined as the pore space accessible from one of the 6 cube sample surfaces. An iterative routine starts searching from one cube side for connected neighbour void voxels in order to find all pores connected to this side. This reconstruction is performed to detect closed pores within the structure that account for the porosity but are not exposed to the gaseous reactants. Table 2 shows exemplary results of total and open porosity for strut samples with various pore former concentrations.   Figure 4 shows all results of the porosities collected from selected 3D tomography reconstructions. εopen is presented with error bars because 6 evaluations are obtained per sample (one from each side). For ϕ ≤ 20 vol%, there is no pore connectivity observed as seen graphically in Table 2. For ϕ ≥ 30 vol%, the majority of the pores are connected and the pore network passes through the entire cube sample. For 50 vol%, practically every pore is connected to the pore network since the open porosity and total porosity are almost the same εopen ≈ εstrut. For high RPC porosity, e.g., εRPC = 0.825, the dual-scale porosity changes by less than 0.1 for ϕ = 50 vol% because only 17.5% of the volume is filled with µm -sized pores within the thin struts (see Figure 4). For lower RPC porosities, e.g., εRPC = 0.459, there is a stronger trend for increasing ϕ since there is more solid to be filled with µm -sized pores. The porosities of the model fit for ϕ = 0 vol% represent εRPC as listed in Table 1.
To investigate the transition from closed pores to interconnected pores in a systematic way, spheres with a small Gaussian size distribution of 2-3 voxels are randomly placed within the 3D volume. Pores of such artificially generated structures get connected between εstrut = 0.2 and 0.3. Such connectivity behavior is consistent with empirical correlations of the effective gas diffusivity within porous carbon [35]. Figure 4 also shows the comparison between the numerically determined εopen and the experimentally measured values by mercury intrusion porosimetry (MIP) [36,37], as reported by Furler et al. [15]. The agreement is reasonable well. For ϕ = 30 vol%, the experimentally measured εopen is lower than the numerically determined one due to poor pore connectivity across coating layers, as seen in Figure 5. This anisotropic region results from the 2-step coating applied during the fabrication process, and it is not considered in the determination of morphological properties and effective thermal conductivity within isotropic regions.

Pore Size Distribution
Pore size distribution is determined by applying a morphology-altering algorithm to the digitally segmented 3D structure consisting of an inversion of the solid and void space, followed by erosion and dilation with successively increasing spherical elements of diameter d with In a last step, the algorithm inverts back the solid and void space. The cumulative pore size distribution 1 − F(d) is defined as the ratio of the opening-closing porosity, εoc(d), and the original porosity [27]: The pore size distribution is then calculated as: . Figure 6 shows the cumulative pore size distribution (left y-axis) and pore size distribution (right y-axis) as a function of d.· f(d) agrees qualitatively well with the values obtained by mercury intrusion porosimetry measurements. The mean pore diameter is then defined in such a way to split the area under the pore size distribution curve into two equal areas [27].
The mean pore diameter for porous strut samples with various concentrations of pore former in the range 10-50 vol% is listed in Table 2. As expected, dmean remains independent of ϕ at around 10 µm because the same pore forming material is used to manufacture all samples. Figure 6. Cumulative pore size distribution (left y-axis) and pore size distribution (right y-axis) obtained by morphology operations with spherical structuring elements of diameter d.

Specific Surface Area
The specific surface area (SSA) is determined in three different ways: (1) using statistical two-point correlation function computed on the 3D segmented structures with an in-house Fortran code; (2) resampling of the phase interface area with a surface mesh-based algorithm using the open source software ImageJ (version 1.47v, Java 1.6.0_20 (64-bit)) [38] extended with the free BoneJ plugin (version 1.3.11) [39]; and (3) using the actual phase interface area of tetrahedral 3D meshes generated with an in-house Fortran code [40]. The in-house mesh generator covers the void and solid domain with tetrahedral elements for unstructured body-fitted grids and subsequently refines the elements at the phase boundary. These 3D meshes are later used to solve the steady-state energy conservation equation to determine the effective thermal conductivity. The two point correlation s2(r) is a statistical function that indicates the probability of two arbitrary points Ψ(r) and Ψ(r + rŝ) separated by the distance r to be in the void phase [41]: where Ω is the solid angle and V the cube volume. Porosity and specific surface area are then calculated using the following expressions [27,30]: A code loops once through the entire structure (x-, y-, z-direction), counts each void voxel, and computes s2(r = 0) in a digitally exact manner. Additionally, the code counts for each void voxel the number of direct neighbour void voxels in all 6 directions (-x, x, -y, y, -z, z) and computes 6· s2(r = 1). The volumetric specific surface area is then calculated as: The specific surface area is presented in two different units: per total volume (fluid+solid phases), . The density of ceria is ρCeO 2 = 7.22 g/cm 3 [21,42].
Dual-scale A0 is calculated by multiplying A0 determined for the struts by the solid volume fraction of the RPC: . A0 of RPC with non-porous struts is converted to ssa by: Table 3 lists ssa of RPC with non-porous struts obtained from the original tomography scans and digitally dilated struts. A0 and ssa are plotted as a function of ϕ in Figure 7a,b, respectively. Of special interest is the open ssa, which is directly related to the surface area reachable by the reacting gases for conversion of CO2 and H2O to CO and H2. The open ssa is calculated for the struts without closed pores. Open A0 and ssa are presented with error bars (standard deviation) for the data evaluated. Surface areas calculated from two-point correlation and phase interface area of the 3D meshes (black and white symbols) lie within close proximity to one another, whereas those calculated using ImageJ (grey symbols) are higher. This is because the two-point correlation is based on a statistical model leading to smoothing effects and the 3D meshes actually contain a smoothing algorithm for the phase interface, whereas the resampled phase interface of ImageJ incorporates fine mesh surface irregularities.  Table 2). Trend lines are plotted defining the mean value of all different calculation methods for total (solid) and open (dashed) A0 and ssa.

Heat Conduction Modelling
The governing steady-state heat conduction equations within the solid phase and the stagnant fluid phase are given by: Fluid phase: where ks and kf are the thermal conductivity of the solid and fluid, respectively. The cubic domain is schematically shown in Figure 8. The boundary conditions are given in Equations (13)- (17). An inlet and outlet temperature is set (Thot > Tcold) to provide a steady heat flux through the two phases with length L. Lateral walls of the sample cube are adiabatic. Local thermal equilibrium is assumed at the phase interface. Heat flux across the interface is driven by the temperature gradient and the thermal conductivities in each phase at the interface.
Inlet temperature: Outlet temperature: Adiabatic lateral walls: 0 q  n (15) Local thermal equilibrium at phase interface: sf TT  Heat flux across phase interface: The governing volume-averaged steady-state equation for effective heat conduction within the isotropic porous structure reduces to one equation [22,23]: The effective thermal conductivity is calculated using the 1D Fourier's law and the heat flux determined by DPLS: where keff is the effective thermal conductivity of the cubic porous structure, q the effective heat flux at the inlet or outlet, and Aflux = L 2 is the inlet or outlet area constraint with Thot or Tcold, respectively. The methodology for determination of keff for dual-scale porous structures is schematically shown in Figure 9. In a first step, keff of the strut with µm-sized pores (keff,strut) is determined according Equation (19). In a second step, this keff,strut serves as an input for the solid domain of a further simulation performed with the mm-sized pores of the RPC. Numerical DPLS are performed for RPC with non-porous/porous struts with/without digital strut dilation for different fluid-solid thermal conductivity ratios ranging from 10 −5 up to 1. The cases covered include 4 RPC with non-porous struts (original scan and digitally dilated struts with 2, 5, and 10 voxels) and 16 RPC with porous struts (i.e., 4 RPC, each with 4 different strut porosities ϕ = 10, 20, 30, 50 vol%). Simulations are performed using a commercial computational fluid dynamics (CFD) software (ANSYS ® Academic Research, release 14.0). Initially, grid resolution study is performed, indicating convergence for structures containing element sizes between 0.57 µm at fluid-solid interface to 2.28 µm within bulk for the µm-size pores within the struts, and between 62.0 µm at fluid-solid interface to 247.9 µm within bulk for the mm-size pores of the RPC. Typical number of elements is 20 million, yielding an error of less than 1% compared to the finest mesh tested of 40 million. The correctness of the DPLS was verified by solving simple geometrical cases with exact analytical solutions, while its accuracy was fine-tuned by grid refinement. Figure 10 shows the ratio of the effective thermal conductivity to the solid thermal conductivity vs. the ratio of the fluid-to-solid thermal conductivity for a single porous strut (ϕ = 50 vol%), a RPC with non-porous struts, and a RPC with porous struts (ϕ = 50 vol% and 0 mm strut dilation). The analytical curves for serial and parallel slabs are indicating the maximum and minimum possible heat flux [43,44]. These exemplary simulation results correspond to εstrut = 0.410 and εRPC = 0.825, leading to εdual = 0.897. keff decreases with increasing porosity and decreasing kf/ks. For kf/ks < 10 −3 , keff does not significantly change anymore, indicating heat conduction dominated by the solid domain (e.g., for vacuum applications). In that range, the ratio of keff for RPC with non-porous struts to keff for RPC with porous struts is 2.4. As expected, this ratio approaches 1 for increasing kf/ks as the thermal conductivities of the fluid and solid phases approach each other. The serial and parallel heat conduction mode of the lumped fluid and solid material bracket the minimum and maximum possible heat flux (also called Wiener lower and upper bound), respectively [43,44]. As expected, the simulation results (black symbols) are between the minimum (white symbols with dashed line) and maximum (white symbols with solid line) possible heat flux for each porosity. Table 4 lists various analytical models for keff [31,45]. For simplicity, analytical equations are given in terms of fs η kk  , and eff eff s kk  . Several models allow fitting with geometrical shaping parameters. They were least-squares fitted to three different sets of simulation data: (1) keff,strut for a single porous strut (black squares shown in Figure 10); (2) keff,RPC for a RPC with non-porous struts (black circles shown in Figure 10); and (3) keff,dual for a RPC with porous struts (black triangles shown in Figure 10). To identify the model which agrees best with all data sets, an overall least-squares approximation, keff,all, was fitted for all simulation data. Figure 10. Ratio of the effective thermal conductivity to the solid thermal conductivity vs. ratio of the fluid-to-solid thermal conductivity for a single porous strut (ϕ = 50 vol%), a RPC with non-porous struts, and a RPC with porous struts (ϕ = 50 vol%). The analytical curves for serial and parallel slabs are indicating the maximum and minimum possible heat flux.

Boomsma and
Poulikakos [57]   Calmidi and Mahajan [60]   Dul'nev and Zarichnyak [22,30,61,62]  The root-mean-square error (RMS) is defined to compare keff calculated by the analytical models with that determined by our simulation:  (20) where n is the number of data points per data set over the entire range of keff/ks and fs kk indicated in Figure 10. Table 5 lists the RMS for the three different simulation data sets and for keff,all. Only those models giving an RMS < 5% are considered appropriate. Models 1 to 9, which are not using any geometrical shaping parameter, give RMS > 10%. However, some models perform comparatively well for the prediction of keff,strut with a RMS < 5%: Hashin and Shtrikman upper bound [47] (3.0%), Russell [49] (4.8%), Loeb [50] (2.6%) and Maxwell [51] (3.0%) which is consistent with the findings by Petrasch et al. [31] for SiC foams. Serial slab and Hashin and Shtrikman lower bound give very inaccurate predictions of keff (RMS > 200%). This is because the serial bound model assumes no direct connection of solid paths between heat inlet and outlet area, which is obviously not the case for connected, but tortuous strut paths. The Schuetz-Glicksmann model [54,55] [59] for all fitting parameters. The empirical model of Calmidi and Mahajan [60], shown in Figure 11b, is capable of predicting all three data sets and an overall data sets with a RMS < 5%. The model of Dul'nev and Zarichnyak [62] gives only keff,strut with a RMS < 5%. Dul'nev and Zarichnyak [22,30,61,62] propose a model using a linear combination of the Wiener lower and upper bounds with empirical fitting parameter, f, for weighting linear combination which is also called three-resistor model. However, if keff is fitted individually for each structure (porosity), an inverse trend of f is observed with porosity. Therefore, the three-resistor model is then extended by describing f as a 2nd-order polynomial function with porosity. Such extended three-resistor model, shown in Figure 11c, predicts keff,strut with RMS = 0.3% instead of 3.9%, keff,RPC with RMS = 1.6% instead of 9.2%, keff,dual with RMS = 2.2% instead of 9.8%, and keff,all with RMS = 2.3% instead of 12.6%. The three fitting parameters describing f with a 2nd-order polynomial function (c0, c1, c2) are replaced to allow the serial and parallel resistance, as well as their combination, to linearly scale with porosity, as shown schematically in Figure 12. Least-squares fitting of this modified three-resistor model, shown in Figure 11d, delivers the most accurate predictions: keff,strut with RMS = 0.1%, keff,RPC with RMS = 1.1%, keff,dual with RMS = 1.4%, and keff,all with RMS = 1.3%. The modified three-resistor model shows the best performance in prediction of keff with overall RMS < 1.5%. Fitting parameter a and b allow the lumped fluid and solid parts to deviate from actual ε within the parallel and serial slabs, respectively. Fitting parameter c allows linear combination of the serial and parallel slab to deviate from ε. This gives some degree of freedom for capturing different tortuous regions for a high porosity range (0.09 < ε < 0.9) and predicts the effective thermal conductivity more accurately compared to linear (or non-linear) combination of parallel/serial bounds and to Miller's bound model. Table 5. Root-mean-square (RMS) error of analytical models compared to three simulation data sets and to all simulation data.  Figure 11. Ratio of the effective thermal conductivity to the solid thermal conductivity vs. ratio of the fluid-to-solid thermal conductivity, obtained by our simulation and by least-squares fitted models for structures with a high range of porosities (0.09 < ε < 0.9). Shown are: (a) Miller bound model [59]; (b) Calmidi and Mahajan model [60]; (c) extended three-resistor model; and (d) modified three-resistor model. The modified three-resistor model predicts keff with the lowest RMS.

Summary and Conclusions
High and low resolution computer tomographic scans were performed on complex reticulated porous ceramics (RPC) structures to capture the 3D digital representations of their dual-scale porosity in the mm and µm range. The CT scans were processed with a Gaussian blurring filter for a clustering-based image thresholding of the void and solid phases using Otsu's method. The struts containing µm -size pores were digitally dilated with spherical structuring elements generating structures with different thickness and porosity. The morphological properties analyzed include porosity, pore size distribution, specific surface area, and pore connectivity within representative sample volumes of the isotropic strut regions and of the RPC. The total strut porosity was linearly dependent on the concentration of pore forming agent, and no pore connectivity was observed for concentration less than 20 vol%, consistent with mercury intrusion porosimetry measurements. A well-connected pore network results in high specific surface area and penetration of reactant gas for high fuel production. The effective thermal conductivities of a single porous strut, a RPC with non-porous struts, and a RPC with porous struts (dual scale) were determined by direct pore level simulations of the heat conduction equation with a CFD code. Values were compared to predictions by analytical models over a wide range of porosities. Models without shaping parameters were generally inaccurate (overall RMS > 10%). Miller's model with two shaping parameters predicted keff with RMS error below 2.1% and the modified three-resistor model with three empirical fitting parameters predicted keff with a RMS error below 1.5%. These analytical correlations are applicable to RPC with porosities in both the strut's µm-scale and bulk's mm-scale ranging from 0.09 to 0.9.
The morphological properties and effective thermal conductivity determined in this work serve as an input to volume-averaged models for the design and optimization of solar chemical reactors.

Acknowledgments
We gratefully acknowledge the financial support by the Swiss Competence Center Energy & Mobility, the Helmholtz-Gemeinschaft Deutscher Forschungszentren (Virtuelles Institut SolarSyngas), and the European Research Council under the European Union's ERC Advanced Grant (SUNFUELS-No. 320541).