## Abstract

We describe a three-dimensional imaging and analysis study of eight industrial cellular foam morphologies. The foam morphologies were generated by differing industrial processing methods. Tomograms are acquired on an X-ray micro-computed tomography facility at scales of approximately equal to at resolutions down to 7 μm. The image quality is sufficient in all cases to measure local structure and connectivity of the foamed material, and the field of view large enough to calculate a range of material properties. Phase separation into solid and porous components is straightforward.

Three-dimensional structural characteristics are measured directly on the porous and solid phases of the images. A number of morphological parameters are obtained, including pore volume-to-surface-area ratio, connectivity, the pore and solid phase size distributions defined by maximal sphere openings and chord length measurements. We further calculate the pore size distribution associated with capillary pressure via simulating of mercury drainage on the digital images.

The binarized microstructures are used as a basis for calculations of transport properties (fluid permeability, diffusivity and thermal conductivity) and elastic moduli. From the data, we generate property–porosity relationships for the range of foam morphologies imaged and quantitatively analyse the effects of porosity and microstructure on the resultant properties of the foams.

We compare our numerical data to commonly used theoretical and empirical property–porosity relationships. For thermal conductivity, we find that the numerical results agree extremely well with an empirical expression based on experimental data of various foams. The upper Hashin–Shtrikman bound also provides an excellent prediction of the data across all densities. From simulation of the diffusivity, we can define the tortuosity of the pore space within the cellular solid. We find that different processing methods lead to strong variations in the tortuosity of the pore space of the foams. For elastic properties, our results show that for the Young modulus, *E*, both the differential effective medium theory and the classical correlation give a good correlation. Assuming a constant Poisson's ratio leads to reasonable agreement. The best correlation for is given by assuming a slight variation in as a linear function of porosity. The permeability of the foams varies over three orders of magnitude. Correlations for permeability based on the classical Kozeny–Carman equation lead to reasonable agreement, except at the lowest porosities. Permeability estimations based on mercury porosimetry give excellent agreement for all foams.

## 1. Introduction

Manufactured cellular materials (e.g. polymer, ceramic or metallic foams) are an extremely attractive option as materials engineered for a range of applications ranging from light weight structures to packaging, insulation and crash protection. The useful properties of cellular solids are a direct consequence of their microstructure. It is important therefore to link the physical properties of cellular solids to their density and complex microstructure in order to understand how their structure can be optimized for a given application.

Relevant aspects of the structure of cellular materials include porosity or apparent density, pore and solid shape and size, and the type and frequency of interconnections between pore and solid regions. These features, many of which unfortunately lack precise definition, comprise the morphology of the cellular solid. Accurately predicting properties from microstructural information requires an accurate quantitative description of the material and the ability to solve for elastic and transport properties on large three-dimensional grids. In the absence of this complete structural characterization, attempts to correlate properties of cellular solids have been limited to bounding methods (Hashin *et al*. 1962), effective medium theories (Schwartz 1994; Berryman 1995) and empirical relationships (Scheidegger 1974; Gibson *et al*. 1988; Bauer 1993). None of these are entirely satisfactory.

Effective medium theories (Hashin 1983) are based on formulae derived by extending exact results for small fractions of spherical or ellipsoidal inclusions to higher porosities. A drawback of this approach is that the microstructure corresponding to a specific formula is not precisely known and may not be realistic; agreement or disagreement with data can neither confirm nor reject a particular model (Berge *et al*. 1993). An advantage of bounding methods (Milton 1981) is that they incorporate microstructural information and can be applied to arbitrary complex structures. Bounds are most useful when the constituent materials have similar properties. Unfortunately, in the case of porous cellular solids, bounds are quite broad due to the large contrast in elastic and transport properties between the pore and the matrix phase. This severely limits their predictive power.

From a design viewpoint, empirical formulae are most popular; they provide a simple and convenient, but deceptive, form of summarizing extensive experimental data. In most empirical formulae for cellular solids, the primary microstructural quantity controlling the properties of a cellular solid is the apparent density, , where is the density of the cellular material and is the density of the corresponding bulk phase. For example, experimental results indicate that the Young modulus *E* of a cellular solid with solid modulus is given by the relation (Gibson *et al*. 1988)(1.1)where *C* is a constant. A similar form (Bauer 1993) is used to estimate the reduced thermal conductance of foamed materials, where and are the thermal conductances of the foamed porous material and solid matrix, respectively. Diffusivity in a cellular solid is also estimated by a power law relationship (Tomadakis & Sotirchos 1993). Permeability is commonly estimated via the Kozeny–Carman formula (Scheidegger 1974):(1.2)where the porosity , *S* is the specific surface area of the material, *τ* is the tortuosity of the medium and *C* is a constant. Most experimental and theoretical work has been devoted to establishing the phenomenological correlation coefficients *n*, *τ* and *C* in equations (1.1) and (1.2) for different classes of material. In these equations, all the explicit microstructural information describing the cellular solid is subsumed in the apparent density and all other geometrical and topological details are ignored. Lacking a rigorous connection with microstructure, the formulae do not offer predictive or interpretive power, and can at times fail to offer physical insight. For this reason, while the engineering potential of various materials is considerable, our ability to design and optimize structures is still very much *ad hoc*, since local structure and mechanical/transport properties have not been measurable.

An alternative approach is to computationally solve for the elastic and transport properties directly on digitized models exhibiting complex microstructures (Adler *et al*. 1990, 1992; Roberts & Garboczi 2000, 2000; Knackstedt *et al*. 2003). Computational techniques have progressed to the point where material properties such as conductivity (Arns *et al*. 2001*b*), diffusivity (Schwartz *et al*. 1994), permeability (Martys *et al*. 1999; Arns *et al*. 2004*b*, 2005) and linear elasticity (Roberts & Garboczi 2000; Arns *et al*. 2002) can be calculated on large three-dimensional digitized grids (over 10^{9} voxels). As input to these methods, statistical models of complex porous materials have been proposed. Most model systems are based on matching lower-order morphological information (e.g. volume fraction and two-point correlation function). Results have shown that although the two-point correlation functions of a reference and reconstructed system are in good agreement, this does not ensure that the microstructures and resultant properties of the two systems will match (Adler *et al*. 1990, 1992; Roberts & Garboczi 1996; Biswal *et al*. 1999).

Direct three-dimensional imaging techniques now enable researchers to accurately describe the complex structure of materials via X-ray micro-computed tomography (μCT; Flannery *et al*. 1987; Dunsmuir *et al*. 1991; Spanne *et al*. 1994; Elmoutaouakkil *et al*. 2002; Sakellariou *et al*. 2004*a*,*b*). These techniques provide the opportunity to directly measure the complex morphology of cellular solids in three dimensions at resolutions down to a few micrometres. In parallel with computational modelling on large grids, it is now possible to replace synthetic images derived from statistical models with real images of complex cellular solids and base calculations directly on the three-dimensional microstructure. This has been done previously for the study of the geometric, transport and mechanical properties of sandstones (Arns *et al*. 2001*a*,*b*). In the present paper, we extend this methodology and describe three-dimensional imaging and numerical analysis study of cellular foam samples. Tomograms generated on the μCT have sufficient resolution to obtain good quality images and a sufficient field of view to calculate a range of structural and material properties. Phase separation into solid and porous components is relatively straightforward. Computational results are derived directly from the digitized tomographic images (segmented images); these include measurements of the pore/solid phase morphology, transport properties (fluid permeability, thermal conductivity and diffusivity) and mechanical modulus of the samples. We compare our numerical data to commonly used theoretical and empirical property–porosity relationships. For conductivity, we find that the numerical results agree extremely well with experimental data for various foams. The upper Hashin–Shtrikman (HS) bound provides an excellent prediction of the data across all densities. Diffusive properties lead to a measure of the pore space tortuosity. Variations over an order of magnitude in tortuosity are noted for foams at the same porosity. Accordingly, no single empirical equation can be used to predict diffusive properties. For elasticity, we find that for the Young modulus, *E*, both the differential effective medium (DEM) theory and the classical correlation (Gibson & Ashby 1997) gives a good correlation. Assuming a constant Poisson's ratio leads to reasonable agreement. The best correlation for is given by assuming a slight variation in as a linear function of porosity. For permeability, we find that estimates based on a length-scale derived from mercury drainage experiments are considerably more reliable than estimates based on surface area.

In §2, we describe the foam samples studied and detail the methodology of the acquisition of μCT images and the image analysis on the resultant tomograms. We then describe the numerical techniques used in the analysis of the foamed material morphology and the calculation of the transport and mechanical properties. Section 4 presents results of these studies and compares to current theoretical/empirical methods. We conclude and discuss possible future work in §5.

## 2. Experimental

### (a) Foamed solid samples

A set of eight microcellular polyurethane foams with a polyester basis were studied: the systems contain urethane groups and also urea. The eight foams were generated by four different industrial processing methods labelled types 1–4 in table 1. Photographs of two of the samples (A, type 1 and J, type 3) imaged in this study are shown in figure 1. The nominal (industrial design) density for each of the eight samples is also listed in the table. The samples exhibit a broad range of densities, .

A composite of two-photon laser confocal images of the eight samples studied is shown in figure 2. All images are at the scale of 1 mm^{3}. From the figure, we observe that the structure of the eight foams is reminiscent of a porous liquid foam; all mimic a closed-cell structure with near spherical bubble-like pores. The mean and distribution of the bubble sizes clearly vary from sample to sample.

### (b) Image acquisition

Tomography is a technique that generates a dataset, a tomogram, which gives a three-dimensional representation of the structure and variation of composition within a specimen. Each three-dimensional volume element in the tomogram is a voxel. X-rays from a micro-focus X-ray source are used to probe the specimen and an X-ray camera is used to record the X-ray transmission radiograph. To generate the tomogram, a series of radiographs are collected at different viewing angles by rotating the specimen. This set of radiographs, projection data, is processed with a reconstruction algorithm to generate the tomogram.

Both the X-ray source and detector are optimized for high resolution and for maximal field of view. The experimental apparatus was developed and built in-house and has a cone beam geometry (Sakellariou *et al*. 2004*a*,*b*). By moving the position of the rotation stage and X-ray camera, magnifications can be set between ×1.1 to over ×100. The limiting resolution is determined by the X-ray source which is around 2–5 μm, depending on operating voltage (30–120 kV). The X-ray camera is capable of acquiring radiographs with 2048^{2} pixels at a depth of 16 bits per pixel. The specimen rotation stage has an angular accuracy of 0.001°. The large dynamic range of the X-ray camera is ideal for the discrimination of fine features in complex materials. Since the source generates polychromatic X-rays, filters are used to minimize beam-hardening artefacts in the tomogram. The images of the foam were taken at 40 kV.

Typically, the camera distance is positioned at the maximum possible distance (2500 mm), which minimizes the X-ray cone angle, thus reducing reconstruction artefacts (Feldkamp *et al*. 1984; Sakellariou *et al*. 2004*a*) associated with cone beam tomography. One disadvantage is that the X-ray flux on the camera is reduced and image acquisition takes more time. For the work reported in this paper, the specimen distance was 250–350 cm and the source-to-camera distance was 2500 mm. An acquisition time of approximately 7 h was typical for the 1024^{3} voxel tomograms collected. The voxel resolution for each tomogram is detailed in table 1.

The collected data are optimally sampled, i.e. the number of radiographs in a tomographic series is approximately equal to times the number of pixels across the radiograph (Sakellariou *et al*. 2004*a*). For each 1024^{3} voxel tomogram, 1500 projections are required. This equates to the collection of 3 Gbytes over 7 h. The tomographic series are reconstructed using algorithms based on the Feldkamp technique (Feldkamp *et al*. 1984). The data is processed using supercomputer (see http://nf.apac.edu.au/facilities/sc/hardware.php) and takes 16 CPUs approximately 2 hours to generate the 1024^{3} voxel tomogram (Sakellariou *et al*. 2004*a*,*b*). An example of the X-ray intensity histogram and a greyscale X-ray density map for a slice within a tomogram of sample A is given in figure 3*a*,*b*. Density maps for all eight samples are shown in figures 7 and 8. Movies illustrating the full sequence through selected samples can be downloaded at http://xct.anu.edu.au.

### (c) X-ray density map and phase identification

The tomographic image consists of a cubic array of reconstructed linear X-ray attenuation coefficient values, each corresponding to a finite volume cube (voxel) of the sample. An immediate goal is to differentiate the attenuation map into distinct pore and solid phases for each of the samples imaged. We illustrate this process for foam A in figure 3. This foam is characterized by reasonably homogeneous spherical pores embedded within a solid matrix. From the original 1000^{3} image, we extract a subset (the region of the image which contains only foam) for analysis. The intensity histogram (figure 3*a*) shows two distinct peaks associated with the two phases. The peak centred around 27 500 is associated with the polymeric material. The lower peak around 17 000 is associated with the pore phase. For an intensity histogram with two distinct phase peaks, it is sufficient to do a simple threshold (22 500) followed by an isolated solid cluster removal to remove noise artefacts (small bits of solid suspended in pores). We also remove isolated pores within the solid phase that are made up of up to 20 voxels. Comparison of the grey scale and binarized image of a slice of foam A is shown in figure 3*b*,*c*. The resultant solid phase fraction at the measured attenuation cutoff is , which is in excellent agreement with the nominal sample solid fraction which was . Most samples exhibited two distinct peaks and the phase separation method was straightforward. In three cases, samples D, I and J, where the thickness of the solid wall was in some cases only a small number of voxels, a new three-step technique for image enhancement and segmentation was used (Sheppard *et al*. 2004). In all cases, the resultant image density is in good agreement with nominal and experimental values (see table 1).

## 3. Numerical simulation and error analysis

A microstructure defined by a digital image is already naturally discretized and lends itself immediately to numerical computation of any number of properties. In this subsection, we describe the numerical methods used to calculate the morphology, conductivity, diffusivity, permeability and elasticity of the digital images of foamed materials.

### (a) Morphology

#### (i) Mean sizes and connectivity

From a single three-dimensional digital image, one can directly measure a range of morphological parameters, including porosity *ϕ*, pore/solid volume-to-surface area . Evaluations are made by simply counting the number of voxels and faces on the image (Bieri & Nef 1984; Arns *et al*. 2001*a*). The counting does overpredict the surface area; e.g. for a digital sphere, the correction factor is approximately 3/2 (Martys & Bentz 1992). This correction is used to rescale the surface area. Assuming spherical pores/struts, a mean pore/strut diameter can be estimated from the ratio:(3.1)Analysis of the pore space connectivity via a burning algorithm (Sheppard *et al*. 2004) is undertaken to measure the amount of the pore space that is continuous versus the amount of isolated porosity.

#### (ii) Pore and solid size distribution and capillary pressure

A model-free technique for defining the pore and solid size of a three-dimensional image is based on defining locally, for every point within the structure, the diameter of the largest sphere, which fully lies within the pore or solid phase and covers that point. More complete and generic descriptions of the basic concepts and techniques are given elsewhere (Hilpert & Miller 2001; Thovert *et al*. 2001). An example of the covering sphere map for the pore and solid phases of a binarized foam image in figure 4*a* is shown in figure 4*b*,*c*.

Mercury capillary pressure simulations can be performed directly on voxelated images using a similar technique (Coles *et al*. 1998; Hilpert & Miller 2001). At a fixed capillary pressure (pore entry radius), we consider all the spheres which have radius greater than or equal to the equivalent pore entry radius. Starting with the largest sphere and incrementing the sphere radius downwards (incrementing capillary pressure upwards), the mercury saturation *S* is measured as the subset of all spheres that have invaded the pore space from the outer surface of the pore volume. This results in a curve for the imaged volume. When simulating mercury porosimetry, we use a surface tension () and assume all outside faces are connected to the non-wetting phase. From the simulation of mercury capillary pressure, one can also measure a critical length-scale associated with the diameter of the largest sphere which can percolate across the pore system (Katz & Thompson 1987). This length-scale, , has been used in empirical correlations for the permeability of porous materials.

#### (iii) Pore and solid chord lengths

A strong discriminator of structure of porous materials is based on the chord length distribution, which gives the probability that a randomly chosen chord in phase *i* has length *r*. A chord is defined as any line segment which lies entirely in phase *i* with endpoints at the phase interface. We measure the chord lengths distributions, and the mean chord length, , along the three orthogonal axes defined by the digital tomograms for both the pore and solid phases.

### (b) Material properties

Material properties are calculated directly on the binarized images of the eight cellular foam specimens.

#### (i) Thermal conductivity simulation

The thermal conductivity calculation is based on a solution of the Laplace equation with charge conservation boundary conditions. Radiant heat transfer between the pore walls is assumed negligible compared to conductive heat transfer in the gas and matrix phase. The three-dimensional voxel microstructure is first converted into a network of resistors by connecting each pair of adjacent voxels by a resistor. The conductance of each resistor represents the effective thermal conductivity between the voxels. We assign to resistors in the matrix phase of the polymer foam the thermal conductivity and to resistors in the pore phase the thermal conductivity of air, . The effective thermal conductivity at the boundary between phases is based on the harmonic mean of the conductivities of the two voxels. A temperature gradient is applied in each coordinate direction and is held constant during the entire calculation. The temperature field is established by calculating the temperature for each node based on the heat flux equation. The system is relaxed to a steady state using a conjugate gradient technique to evaluate the field. The overall thermal conductivity of the sample is reported as a function of porosity along the three orthogonal axes.

#### (ii) Diffusivity simulation

Diffusion within the porous medium is also defined by the Laplace equation, where a concentration difference is applied across the system and the total mass flux can be related to an effective diffusion coefficient. In the simulation, we set the free diffusion coefficient in the pore space to and the diffusion within the solid matrix phase to zero. The effective diffusivity is reported, where is the diffusion coefficient within the porous cellular material.

From the diffusion coefficient, the tortuosity *τ* of the conductive flow channels or of the diffusion process can be defined in the pore space (Schwartz *et al*. 1993);(3.2)

#### (iii) Permeability simulation

The air/fluid phase permeability calculation is based on the lattice-Boltzmann method (LB; Qian *et al*. 1986; Martys & Chen 1996). The LB approach is a mesoscopic approach for computational fluid dynamics, where the macroscopic dynamics of the solution of a discretized Boltzmann equation can be shown to match the Navier–Stokes equation. Owing to its simplicity in form and adaptability to complex flow geometries, like the presence of solid–fluid and fluid–fluid boundaries, one of the most successful applications of the LB method has been to flow in porous media (Ferreol & Rothman 1995; Martys & Chen 1996). In this method, the fluid is modelled as a collection of particles moving on a regular lattice, with each lattice point connected to its nearest and next-nearest neighbours. In our implementation, we use the D3Q19 (three-dimensional lattice with 19 possible momenta components; Qian *et al*. 1986). Mirror image boundary conditions (Martys *et al*. 1999) are applied in the plane perpendicular to the flow direction, and all simulations were performed on an system; permeability is measured in the central (*L*^{3}) subset. The physical boundary condition at solid–fluid interfaces is the no-flow condition, which in the LB methods is most simply realized by the bounce-back rule (Martys *et al*. 1999). The pressure gradient acting on the fluid is simulated by a body force (Ferreol & Rothman 1995).

#### (iv) Elasticity simulation

We use a finite-element method (FEM) to estimate the elastic properties of model systems. FEM uses a variational formulation of the linear elastic equations and finds the solution by minimizing the elastic energy using a fast conjugate gradient method. Each voxel is taken to be a tri-linear finite element. A strain is applied, with the average stress or the average elastic energy giving the effective elastic modulus.

In previous work (Garboczi & Berryman 2000; Arns *et al*. 2001*a*), materials with large solid fractions were considered and the digital image was assumed to have periodic boundary conditions. For the cellular systems exhibiting low solid fractions , one cannot impose these boundary conditions. To compute the elastic properties of these low-density materials, we simulate the system using effective medium boundary conditions (EMBC). We first assign a third (boundary) phase to the tomogram of considerable thickness. For systems studied in this paper (approx. equal to subsets), the boundary phase is set to a thickness of 15 voxels, leading to the simulation of a 230^{3} volume of material. This boundary phase is assigned a less stiff modulus than the solid voxels within the image volume. As an initial guess we use the value given by the HS upper bound for the modulus of the boundary phase. The FEM is then solved for the three-phase (pore/solid/boundary) system and the modulus of the composite computed. This new composite modulus is then assigned to the boundary phase, and the process repeated until the difference between the assigned boundary modulus and the resultant composite modulus is less than 1%. Convergence is achieved within less than 10 iterations (see figure 5*a*). The final result is also independent of the initial choice of modulus. This method was tested for periodic arrays of solid spheres across a broad range of density, . The result for the modulus given by the EMBC and the value given from the original FEM method on the periodic system gives excellent agreement (figure 5*b*).

### (c) Numerical error and choice of representative sample size

In order to obtain accurate numerical results, it is necessary to estimate and minimize three sources of error: finite size effects, statistical fluctuations and discretization errors. This is generally not done due to limitations in computational speed. However, if one wishes to compare computations directly to experiment and rigorous theoretical bounds, such error analyses must be carried out (Roberts & Garboczi 2000; Arns *et al*. 2001*b*). We discuss the three sources of error separately.

The size of the tomogram compared to some statistical length-scale (e.g. correlation length, average cell size) controls how many equivalent cells lie within a single tomogram. The full image has thousands of cells within it. Errors will occur if we use too small a field of view to simulate the real foamed properties. Optimally, we would choose a system size which has an acceptable finite size error, but is still small enough to allow computational prediction of properties. Previous work (Roberts & Garboczi 2000; Arns *et al*. 2001*a*,*b*) indicates that samples of greater than or equal to 6 times the length-scale of major features are needed to obtain reasonable results. From inspection of the foam images, we note that primary features are of the order of –—samples at a scale of >1 mm should therefore give representative data. In the results section, we show that all properties of the foam samples give good averaged behaviour at scales of or voxels. Simulations at this scale are feasible on modern workstations. Excitingly, choosing smaller sample sizes gives one the ability to generate an ensemble of *independent* samples. Measurement at a scale of approximately equal to enables us to generate a number of independent data points from each foam image. We have previously shown (Arns *et al*. 2001*a*,*b*, 2004*a*,*b*) that choosing a number of independent samples from a single tomogram enables one to measure physical properties of sedimentary rock across a range of phase fractions. A single sample can, in principle, give one a representative property versus phase fraction relationship across a range of solid/pore fractions and gives an indication of the variation in properties with or rather than a single data point. A scatter plot of the data obtained for the elastic properties of foam A from all subsamples is given in figure 6*a*.

Second, we discuss statistical uncertainties. To estimate errors, we bin the measurements made on the computational cells at approximately equal to as a function of porosity with bin sizes , dependent on the heterogeneity of the sample. The error bars reflect twice the standard error (s.e.=, where is the standard deviation). The binned data from figure 6*a* is given in figure 6*b*. There is a 95% chance that the ‘true’ result lies between the indicated standard error bars.

We finally consider the discretization error; the error due to the use of discrete voxels to represent continuum objects. Errors due to discretization can include, for example, inaccurate description of curved grain boundaries and closing of narrow pores. It has been previously shown that one can accurately estimate the continuum value of the conductivity (Roberts & Knackstedt 1996; Arns *et al*. 2001*b*) and elasticity (Arns *et al*. 2001*a*) of tomograms by varying the discretization used in the computational scheme and assuming(3.3)where , is the measurement of the thermal conductivity and elasticity at resolution and are the ‘continuum’ results. To measure this effect, we generate realizations of the original tomographic datasets at integer multiples of the resolution of the original image. To generate the images at poorer resolution, we bin voxel clusters of sizes , using a simple majority rule (for 50% occupation, we randomly assign the upscaled binned voxels). In figure 6*c*, we illustrate the effect of finite resolution on the prediction of the bulk elastic modulus *K* of foam A. The scaling of the discretization error with along with the eventual continuum values derived using equation (3.3) are also shown. It is important to note that at the original resolution of the tomograms, the predictions differ from the ‘continuum’ result. At 10 μm resolution, the error can be more than 10% and the error increases with decreasing resolution.

In the remainder of this paper, we will report the ‘continuum’ prediction , . While finite size and discretization effects are noted for the conductivity and elasticity, the results of the permeability simulations exhibit only a small discretizaion error (Arns *et al*. 2004*a*), and the data derived on the original tomograms are reported as the continuum prediction.

## 4. Numerical calculation of properties

### (a) Morphology

The pore and solid size distributions defined by pore/solid volume-to-surface area (3.1), by the maximal sphere algorithm and by the chord length measurement, are summarized in tables 2 and 3 and figures 7 and 8.

Samples A–D (types 1 and 2) exhibit reasonably homogeneous pore size distributions as defined by the embedded sphere algorithm. In contrast, samples I–L (types 3 and 4) exhibit a small number of pores of extremely large pore radii, which skew the distributions to larger values and increase the variances of the distributions. In figure 9, we compare the mean pore size as a function of porosity for the four different types of foams. Type 1 and type 4 foams exhibit very similar mean pore sizes. We note that foams of type 3 have the smallest pores. In all cases, the pore size tends to increase with increasing porosity. The pore chord length distribution exhibits a slight anisotropy and mimics a log-normal distribution. All systems with porosities exhibit a nearly completely interconnected pore space. The lowest porosity sample still exhibits significant (85%) interconnected porosity (see figures 7 and 8 and table 2).

The solid size distribution for all foams is reasonably homogeneous, with the broader distributions shown by the foams of higher solids content. In figure 9, we plot the solid size distribution as a function of porosity for the different foam types. The solid chord distribution exhibits a log-normal shape with a broad distribution of values (see figures 7 and 8). In many cases, the chord length within the solid phase in the *z*-direction shows a slightly larger mean value and distribution than in the other orthogonal planes (see table 3).

### (b) Physical properties

#### (i) Thermal conductivity of foam samples

The thermal conductivity assigned to the solid polymer phase was 0.18 W m^{−1} K^{−1} and the air phase 0.025 W m^{−1} K^{−1}. Figure 10 shows results for the *continuum* thermal conductivity of the foam images measured along three orthogonal directions on subsets of the original tomographic image at sizes (in voxel units) of . The results show a very consistent increasing trend in with increasing solid phase fraction despite the very small variation in for each independent subset. A small amount of anisotropy is noted in some samples with larger values of the thermal conductivity exhibited in the *z*-direction. This is consistent with the increased mean chord length in the solid phase along this axis. Samples A, C and I, which exhibit little anisotropy in the mean solid chord length, exhibit the smallest degree of anisotropy in . In contrast, samples B and J, which exhibit the strongest anisotropy in mean chord length, exhibit a notable shift in .

In figure 11, we plot all the data from the eight samples averaged along the -direction. Notably, there is no clear difference between the thermal conductance/reduced density trend for the eight samples despite the differences measured in foam type and cell size/morphology.

Empirical formulae have been developed to describe the density dependence of the thermal conductivity of foams. A commonly used empirical correlation was derived directly from experimental measurements of the conductivity of a range of porous liquid foam morphologies (Bauer 1993), with . The correlation is given by(4.1)where is a shape factor (Bauer 1993) based on fitting the experimental data to the range of foams. In figure 11*a*, we compare this empirical fit to the numerical data derived from the three-dimensional foam images. The results are in good agreement, with the numerical data approximately lower than given by equation (4.1).

The most common theoretical estimates of conductive properties are given by rigorous bounds. Hashin & Shtrikman used a variational principle to obtain bounds on the magnetic permeability of macroscopically homogeneous and isotropic mixtures of multiple phases. These HS bounds hold for the thermal conductivity of a two-phase material:(4.2)where and and define the two-phase fraction. These rigorous bounds are widely utilized; they are the most restrictive bounds dependent solely on the phase fraction and otherwise independent of detailed microstructure. A second theoretical method for treating the transport properties of inhomogeneous materials is based on an effective medium approximation (EMA). The virtue of the method is that it provides a simple method for dealing with complex geometries and is not limited to low concentrations of inhomogeneities or to weakly varying conductivities. For disordered systems exhibiting bi-connectivity and assuming the cellular material is made up of spherical inclusions, one should use the symmetric EMA (Bruggeman 1935) to estimate the resultant :(4.3)

In figure 11*b*, we compare data from all eight samples with the prediction of bounds and EMA. We note that all the foam data lie within the bounds and that the upper bound gives a *very* good agreement with the numerical data. The EMA prediction tends to underestimate the simulation data. Given the excellent fit of both the upper HS bound and equation (4.1), we compare the two predictions directly in figure 11*a*. This shows that for high densities, equation (4.1) slightly overestimates rigorous bounds. Based on this we conclude that for the cellular solids of the type studied here, the upper HS bound may give the best estimate of thermal conductance across the broad range of .

#### (ii) Diffusivity of foam samples

In figure 12, we plot the effective diffusivity of the pore space as a function of porosity . We note a more pronounced anisotropy in samples of types 2 and 3 (samples B, I, J and L). This is consistent with the increased mean chord length in the pore phase along the *z*-direction of the samples. From equation (3.2), we can derive the tortuosity of the pore space within the samples. This is shown in figure 12*b* for diffusivity along one of the axes. We note that the tortuosity of the samples of types 2 and 3 is low for all values of , while foam types 1 and 4 (samples A–C and K, L) show a significantly larger value for ; for smaller , . The difference in the pore connectivity of two foams with similar porosity, but widely varying tortuosity, over a similar sample size is visualized in figure 13. Foam J () exhibits dense interconnectivity of pores, while foam L () exhibits poor interconnectivity and long path lengths. One can also observe a visible difference from comparing the two-dimensional images in figures 7 and 8. Foam J shows strong interconnectivity between pores in the two-dimensional slice, while the pores are distinctly disconnected in foam L. This is the strongest indicator yet of the difference in the pore structure of the different foam types. Clearly, the different pore structure of the foams will have a strong effect on the resultant fluid transport properties within the samples.

Given the distinct tortuosity of the foam types leading to strongly varying diffusive behaviour, it is clear that one cannot fit the diffusive properties of the foam data via a single empirical fit. The difference in across foam types indicates that this measurement may be very important as an indicator of pore morphology.

#### (iii) Elastic properties of foam

The elastic properties assigned to the solid phase comprising the polymer foam structure were Young's modulus and Poisson's ratio . Figure 14 shows results for the reduced bulk modulus of the foam images simulated on subsets of the original tomographic image. The data show a trend of decreasing modulus with porosity.

From a specification of the volume fraction and constituent moduli, one can obtain rigorous upper and lower bounds on the elastic moduli of any composite material. The so-called HS bounds are given by(4.4)(4.5)where , are the bulk and shear moduli of the phase, respectively. Upper and lower bounds are computed by interchanging the moduli of the high and low modulus components. In the current case, where the pore phase has zero elastic moduli, the lower bound becomes zero, and so only the upper bound is meaningful. Young's modulus and Poisson's ratio are obtained via the relations and .

One common effective medium theory, the differential (DEM) method, is constructed by incrementally adding inclusions of one phase into the second phase with known constituent properties. DEM does not treat each constituent symmetrically, but defines a preferred host material. From the composite host medium, at some porosity value is known. One then treats as the composite host medium and as the effective constant, after a small proportion of the composite host has been replaced by inclusions of the second phase. For a solid matrix host, the coupled system of ordinary differential equations for the moduli is given by (Berryman 1992)(4.6)(4.7)with initial conditions and , and where and are shape-dependent geometric factors (see, for example, table 4.9.1 of Mavko *et al*. (1998)). For the present study, as the pores in the model morphology are near spherical, we use the geometric factors for spherical porous inclusions.

In the self-consistent model (SCM) of Hill (1965) and Budiansky (1965), the host medium is assumed to be the composite itself. The equations of elasticity are solved for an inclusion embedded in a medium of unknown effective moduli. The effective moduli are then found by treating as tunable parameters. The result is given in general form (Berryman 1980) as(4.8)(4.9)where the indices to *P* and *Q* indicate the inclusions of fluid ‘’ and solid ‘’ into a background medium of effective moduli and . The solution for the effective bulk moduli is found iteratively. In the SCM study, we considered the geometric factors for spherical pores and a number of granular inclusion shapes. The SCM produces a single formula in which both components are treated uniformly, with neither material being distinguished as the host to the others. Such a symmetric formula has been thought to be more appropriate in complex aggregates like consolidated granular rocks and was shown (Berge *et al*. 1993) to accurately predict the mechanical behaviour of a sintered glass bead pack. We compare the above three theories to the numerical predictions for *E* and in figure 15.

The SCM predictions strongly underestimate the numerical data for *E*, as the theory predicts a vanishing modulus for . The DEM shows a remarkable agreement with the data. It has previously been noted (Berge *et al*. 1993) that the DEM gave a good fit to data for the elastic behaviour of a porous glass foam with large open pores dispersed in a continuous glass background (Walsh *et al*. 1965). However, for reduced densities of 60% or less, the match to the experimental data became poor with DEM significantly overestimating the predictions. This difference may have been due to the lower-density glass foams exhibiting more needle-like pores. In a second numerical study of a model foam morphology based on overlapping spherical (disc-like) pores (Roberts & Garboczi 2000) with porosities ranging from , it was found that DEM overestimated the data. The results shown in figure 15 show that the prediction of DEM for the Young modulus is in excellent agreement with the numerical data. The DEM prediction for Poisson's ratio slightly underestimates the numerical data. DEM predicts as . In contrast, the Poisson ratio for the foam decreases slightly from 0.3 with decreasing .

In figure 16, we compare our data to empirical structure–property formulae recently derived for model cellular foamed solids. A general empirical fit to the Young modulus has been described by the form (Roberts & Garboczi 2000)(4.10)where and *m* are empirical fitting parameters. Values derived for the fitting parameters for Boolean models of overlapping spherical pores were and for overlapping oblate spheroids (Roberts & Garboczi 2000). In figure 16, we compare these predictions to the simulation data for the cellular foams. In both cases, the fits given in Roberts & Garboczi (2000) strongly underestimate the data. Roberts & Garboczi (2000) do note that the parameters derived for their model systems would be best for . We find that the best fit of the data to equation (4.10) is given by . We plot for interest the classic prediction of Gibson & Ashby (1997) for *open*-cell foams; . While the morphology exhibited by the cellular foams is not reminiscent of a classic open-cell foam, the prediction is surprisingly good (correlation coefficient=0.992). For comparison, we also show the prediction of DEM on the plot, which also gives a relatively strong match to the data.

The empirical form of the fit to Poisson's ratio was given by Roberts & Garboczi (2000):(4.11)In previous work (Arns *et al*. 2001*a*, Roberts & Garboczi 2002), it has been noted that the flow diagram for converged towards for several morphologies and in agreement with effective medium theories and the rigorous behaviour in two dimensions (Day *et al*. 1992). The best fit for the model morphology made up of spherical pores gave (Roberts & Garboczi 2000). This behaviour is not noted in the samples studied. In figure 16, we compare the Poisson ratio derived from simulation on the tomograms to the empirical predictions and give an error estimate. The prediction for a model system made up of spherical pores (Roberts & Garboczi 2000) gives the poorest fit (), while the assumption that the Poisson ratio remains constant is satisfactory (). The best fit for the foam data to equation (4.11) is ().

Based on these results we conclude that for these cellular solids, one would use either DEM theory or the classical correlation to estimate Young's modulus. The best correlation for the Poisson ratio is based on equation (4.11) with , but assuming a constant Poisson ratio may lead to reasonable agreement.

#### (iv) Permeability of foam

In figure 17, we plot the computed continuum permeability *k* of eight foam samples as a function of porosity. Permeability is reported in units of Darcies (). We observe an enormous range of permeability from the eight samples with permeability varying across four orders of magnitude. This range of permeability is commonly observed in geological systems across a smaller range of *ϕ* (Arns *et al*. 2004*b*). In the high porosity samples, the permeabilities of the three foams A, D and K are comparable. Type 3 samples, which exhibited a much smaller tortuosity than types 1 and 4, do not exhibit greatly enhanced permeability. Permeability depends, however, on both the tortuosity of the pore space and on the size of the permeable pathways. Permeability has the dimensions of area and may be thought of as representing the cross-section of an effective channel for fluid flow through the pore space. This smaller difference in the permeability of foam J versus L is therefore primarily due to the smaller pore sizes in type 3 foams (recall figure 8).

Estimates of permeability must involve a measurement of a length-scale relevant to fluid flow. Here we compare two formulae that are widely used to link permeability to pore size in the materials literature. The most common formula used to estimate permeability is the Kozeny–Carman formula, which is given in equation (1.2). Assuming cylindrical pores one can rewrite equation (1.2) as(4.12)We compare the prediction of this equation directly to the numerical data for flow through the high permeability foams in figure 18 for each subset. The prediction of the empirical equation is reasonably good for high and intermediate porosity samples, but strongly overestimates the permeability at lower porosities.

The Kozeny–Carman relationship defines the effective channel size defined by , which has been argued to be poorly representative of the fluid flow within rocks. Katz & Thompson (1987) argued that permeability can be correlated to , a critical pore diameter corresponding to the diameter of the smallest pore of the set of largest pores that percolate through the material;(4.13)where *c* is a constant that depends on the distribution of pore sizes. This method for estimating the permeability may have significant merit, as can be directly estimated from mercury intrusion experiments (Katz & Thompson 1987). The best-fit prefactor in equation (4.13) for the foam data is given by *c*=0.04. This is in good agreement with the prefactor derived for many geophysical samples, where we found *c*=0.042 (Arns *et al*. 2005). Note that this value is significantly different to the value originally proposed by Katz & Thompson (1986) (, see equation (4.13)). We plot the prediction of equation (4.13) in figure 18 and see that the agreement with the experimental data is of much higher quality than the conventional Kozeny–Carman relationship. This is particularly noticed at lower to intermediate porosity samples.

## 5. Conclusion and discussion

It is important to link the physical properties of foams to their density and complex microstructure, and to understand how such properties can be optimized for a given application. These results demonstrate the feasibility of combining three-dimensional experimental imaging of complex foamed materials and computational analysis of the resultant images to generate property–porosity correlations for a range of transport and elastic properties. Combined with the analysis of pore/solid morphology, this may assist the engineer in designing foamed materials for specific applications and gain a better understanding of the role of microstructure and processing variables in determining the material behaviour.

To this point we have shown how material properties such as thermal conductivity, diffusivity, elasticity and permeability can be calculated on large three-dimensional digitized images. We have also illustrated that simulations on microtomographic images agree with experimental property predictions on similar porous foam morphologies. Numerical ‘experiments’ can also give important local information that cannot be measured experimentally. In figure 19*a*, we map the local stress distributions resulting from the numerical simulation of the elastic properties of the foam A sample (Saadatfar *et al*. 2004). Identifying these regions can help identify weak points within the material, giving an immediate indication of the yield strength of the material.

One can also solve for local flow properties within the complex foam images. We do this by solving Navier–Stokes equations using the LB approach. To observe the complex flow paths, we generate three-dimensional visualizations of the fluid pathways on high-end graphics workstations. In figure 19*b*, we show a visualization of the flow through the pore space, where the density of the streamlines represents the flow velocity. Movies illustrating the full rotation sequence are available at http://xct.anu.edu.au.

In this paper, we have shown that microtomographic imaging hardware and computational techniques have progressed to the point where properties such as diffusivity, elasticity, permeability and conductivity can be calculated on large three-dimensional digitized images of real cellular solids. Microtomographic imaging at higher resolutions with faster acquisition times (allowing dynamical effects to be imaged) in parallel with faster computational processing speeds will allow further analysis of the behaviour of these materials. The parallel development of these experimental and computational methods should lead to the ability to image and analyse materials in three-dimensional in minutes. This will provide materials scientists and engineers with a real-time virtual property testing laboratory for complex materials.

## Footnotes

- Received November 2, 2005.
- Accepted January 6, 2006.

- © 2006 The Royal Society