Three-Dimensional Digital Reconstruction of Ti2AlC Ceramic Foams Produced by the Gelcast Method

A digital reconstruction technique is presented that generates three-dimensional (3D) digital representations of ceramic foams created by the foam-gelcasting technique. The reconstruction process uses information that is directly extracted from Scanning Electron Microscopy (SEM) images and offers a 3D representation of the physical sample accounting for the typically large pore cavities and interconnecting windows that are formed during the preparation process. Contrary to typical tessellation-based foam treatments, a spherical representation of the pores and the pore windows of the foams is assumed and a novel hybrid algorithm that combines a variation of Lubachevsky-type and Random Close Packing of Hard Spheres (RCPHS) algorithms has been developed to obtain near-optimum solutions to the packing problem of the spheres that represent the pores. Numerical simulations are performed directly on the 3D reconstructed foams to determine their gas permeability. The model predictions are compared with experimental gas permeability data that were obtained for the physical samples. The pore wall thickness can be treated as the single fitting parameter in the entire reconstruction process, although it is shown that images of sufficient resolution could eliminate the need even for that. The foams that are produced by this method yield quantitatively similar pressure drops with experiments for various superficial velocity values, with a very small deviation in the range of 1.7–2.8%. The proposed methodology could be utilized for the prediction of the permeability and transport properties of complex foamy porous structures, similar to the gelcast-type of foams, from a single SEM image of the foam sample without resorting to serial tomography or other structural information, thus saving considerable time and effort from experimental work.


Introduction
Layered materials are widely available in nature and have been shown to exhibit new and technologically attractive properties when they undergo excessive thinning, delamination or exfoliation, in comparison to their bulk counterparts. As a result, a wide variety of ceramic materials comprised of complex layered structures has appeared, such as MAX phases, characterized by the M n+1 AX n (n = 1-3) chemical composition formula, where M is an early transition metal, A is an A-group element of the periodic table and X is either C or N [1]. Ti 2 AlC is a member of the MAX-phase group of layered ternary carbides and nitrides. The chemical bonding in Ti 2 AlC is anisotropic in nature and can be described as metal-covalent-ionic because of the distribution of the charge density. Moreover, the chains of directional covalent bonding between titanium and carbon (Ti-C-Ti) and layers of the closed that are used in the literature for the generation of dense packings of spheres are classified into two categories-a) the static construction methods and b) the dynamic construction methods. In the static case, the spheres are spatially allocated following a usually sequential deposition process and are considered immobilized after their positioning. In contrast, the dynamic approach uses the initial position of the spheres as an input parameter and assumes short-or long-range interactions between spheres, generating displacements of the spheres and spatial rearrangement of the whole packing, ultimately converging to a packing with the desired properties.
Overall, the static algorithms are very popular due to their direct nature and simplicity. The large value of the foam porosity necessitates the generation of spherical pore packings with a large number density. This can be achieved through several algorithms that are used to generate dense packings of hollow spheres. The majority of those models are based on a simple yet efficient conception, namely that the spheres are placed sequentially in the working domain under the effect of a usually unidirectional force-field and, then, remain fixed in space [19,20]. The static methods, however, can generate relatively homogeneous packings. The packing fractions that are obtained are smaller than the maximum possible fractions for disordered sphere packings [20]. Contrary to the static methods, the dynamic methods are more flexible, albeit more intricate and perplexing. Their complexity is attributed to the fact that the final position of the particles is dependent on the process, which can be either collective or individual. Moreover, the resulting packing fraction strongly depends on the evolution process of the sphere positions.
The dynamic simulation processes are further divided into two major subcategories, namely event-driven simulations and molecular dynamics simulations. In the event-driven simulations, the components of the system under analysis (disks in two dimensions or spheres in three dimensions) evolve independently, unless an event takes place [21]. The sphere evolution is based on the integration of Newton's equation of motion that is used to predict the trajectory of the spheres [22]. Thus, the temporal displacement of moving spheres is directly proportional to their individual velocity. Collision events interrupt this evolution process and the lineal motion of spheres [22]. These events may correspond to an elastic or inelastic collision between two spheres or to the collision of a sphere with the bounding box walls. The collisions introduce an abrupt interchange on the momentum of the spheres. The duration of body contact during collision is usually assumed to be zero. Using this algorithm, the whole system of spheres is updated after each event [22].
The reconstruction process that is developed here yields three-dimensional, digitized representations of different Ti 2 AlC foams that have different pore size distribution and different porosity values. The final goal of this work is the ability to create 3-dimensional representative volumes of the ceramic foams from minimum input data and to derive useful information regarding their internal structure. Effective transport properties can then be extracted accurately using computational treatment only, thus avoiding exhaustive experimentation. The distinctive advantage of this reconstruction algorithm over other algorithms found in the literature [17,18] is that this algorithm can digitally reconstruct foam structures where each pore cell may be interconnected with others via multiple pore windows and other pore cells. Moreover, this algorithm can be used in order to predict the 3D porosity as a function of the wall thickness of the spherical caps (a quantity that can be easily extracted from Scanning Electron Microscopy (SEM) images). This work aims to speed up the characterization procedure for ceramic foams and provide guidance to the tailoring of their structural properties. A modified event-driven algorithm is proposed, which involves an overlap minimization process that produces an optimized hollow sphere packing. The gas permeability of the reconstructed foams is obtained via numerical methods and is found to compare very satisfactorily with experimental values.

Materials and Methods
In this work, a methodology for the 3-dimensional reconstruction of dry foams is developed, followed by the numerical extraction of their effective flow and transport properties. The reconstruction approach is based on a packing of non-overlapping or partially overlapping spheres, whose radii are sampled from the size distribution that is extracted by pore-space analysis on a single SEM image.
Depending on the application, some disk-packing or sphere-packing generation algorithms are found to perform better than others in some aspects of the simulation procedure yet fall behind in some other aspects. Thus, it is desirable to combine algorithms in a way that exploits their strengths and manages to sidestep their weaknesses, formulating a new, hybrid, "composite" algorithm. The hybrid algorithm that is developed here consists of two well-established algorithms for the generation of a hollow sphere packing. The first constituent algorithm is the algorithm proposed by Lubachevsky who worked first on billiard simulations [23] and then on random disk packing problems [24]. The second constituent algorithm is the Random Close Packing of Hard Spheres algorithm (RCPHS) developed by Wu et al. [25], which is a rearrangement algorithm with an optimization subroutine.
The Lubachevsky algorithm uses the concept of following the motion of a number of spheres in a closed domain, similar to a billiard simulation. However, the diameter of the spheres grows gradually in time. Since the volume of the domain is finite, there will be a time when the spheres cannot grow any further, yielding a jammed pack. This algorithm has one parameter representing the growth rate of the spheres compared to the mean velocity of the spheres. In practice, low values of the growth rate parameter are associated with sphere packing that tend to approach the ordered state; however, for larger values of the growth rate a more chaotic dispersion is attained. The Lubachevsky algorithm also has the appealing characteristic of producing dense jammed packaging for monodispersed sphere systems; however, this claim cannot be extended to polydispersed sphere systems [26]. Moreover, due to the fact that this algorithm has been strictly designed for hard spheres, it is not capable of handling jammed packings of soft, partially overlapping spheres without inherent modifications.
The RCPHS algorithm considers an arrangement of spheres with overlap between neighboring spheres. Then, a procedure is applied to relocate each sphere so that overlapping is reduced. When the total sphere overlapping volume cannot be further reduced, all spheres are incrementally shrunk. By consecutive repetition of the relocation and shrinking steps, a non-overlapping or, in any case, a minimized overlapping packing is eventually obtained. However, in the core of this algorithm lies a shrinking step that manipulates the radius of the spheres in order to facilitate the minimization of the overlapping between the spheres. The shrinking step of this algorithm is undesirable in our implementation since the resulting distribution of the sphere radii is not guaranteed to follow the desired one.

Proposed Hybrid Algorithm for Packing Generation
In this work, a hybrid algorithm that combines a variation of the Lubachevsky and the RCPHS algorithms is developed to obtain a near-optimum solution to the sphere packing problem. Because of the fact that the RCPHS algorithm requires an initial arrangement of spheres in order to operate on them and minimize their overlap, a variation of the Lubachevsky-Stilling [23] algorithm is used. At first, the simulation domain is defined paying attention to the conditions that should apply at the boundaries. The spheres can either rebound off the walls of the simulation domain (rigid boundary) or pass through the wall to re-enter through the opposite side (periodic boundary). In the present case, the rigid boundary condition is adopted since it combines robustness and simplicity.
Then, the simulation domain is populated by spheres that are randomly placed in space and have as radii a very small fraction of the finally desired ones. Their velocity components are represented by randomly oriented unit vectors. To simplify the computations, the radii and positions are further normalized by the characteristic domain length. The spheres move with the unit velocity into random directions and are subject to elastic collisions with other spheres or with the domain boundaries. As the sphere sizes gradually increase (exponential growth is assumed) and the domain size is fixed, there will be a moment at which a jammed configuration is obtained. More specifically, starting from a random spatial distribution of very small spheres, they are allowed to transform gradually into non-overlapping or partially overlapping spheres, the radius of which is increased at a constant expansion rate, defined as b rate , while their positions evolve over time according to Newtonian mechanics. The expansion of the spheres is assumed to occur simultaneously on all spheres after every collision event. Thus, as time progresses, each sphere experiences an increasing number of collisions due to the fact that the available free space is reduced gradually. When a jammed stage is reached and no further repositioning according to the previous rules is possible, the spheres are magnified to their original (desired size). The RCPHS algorithm is subsequently employed to minimize the overlapping of spheres by rearranging their positions.

Sphere Collisions
When the spheres are not colliding, they move in straight paths; no forces act on them and there is no acceleration. The collision time between two spherical particles can be determined as follows. Let us assume two spheres i and j with radius R i and R j and velocity → v i and → v j , respectively. If these two spheres are in a collision trajectory, they are going to collide at time ∆t c and the following equation will be satisfied, Equation (1).  (1), is obtained by squaring the vectorial expression for the binary collision of spheres with constant velocity in time. However, during the computation of the roots of Equation (1), numerical errors accumulate during the calculation, attributed to the finite precision of floating point calculations. These errors, if left untreated, inevitably lead to negative time-steps or complex solutions of the quadratic. They also affect the prediction of the positions of the colliding spheres at the time of impact, which results either in spheres overlapping or in positioning of spheres at a small distance from each other or from the wall, at the time of collision. The sphere overlapping case is extremely problematic [22,27], since the particles have numerically entered the infinite energy hard-core [27]. Errors resulting from → r ij ≥ R i + R j do not affect the stability of the algorithm. In order to deal with the floating point precision problem, the algorithm that was proposed by Bannerman et al. [18] is implemented. In this algorithm, an overlap function, F HS , is used at any time step, defined as the left-hand-side of Equations (1) and shown in Equation (2), which evaluates the overlapping between two spheres at each time step. This overlap function obtains negative, positive or zero values for a pair of overlapping, non-overlapping or one-point contacting spheres, respectively. Thus, the problem of sphere collisions is transformed into a problem of finding the roots of F HS = 0. The derivative of the overlap function with respect to the collision time may be utilized in order to predict future overlapping between spheres that are not currently overlapping. By considering the aforementioned properties, overlap in each pair occurs when. This overlap decreases, that is, the spheres move away from each other, when ∂F HS /∂∆t c > 0. Conversely, the spheres approach each other and the overlap increases when ∂F HS /∂∆t c < 0. Thus, the search for the earliest occurring collision pair is led by the smallest non-negative ∆t c that is subject to the overlap function conditions: Thus, by following the definition of the overlap function and employing the following algorithm sequentially, for each pair of spheres i, j the ∆t c,i,j is calculated, which is the time increment for the mutual collision of the spheres.
If none of the above is satisfied, then ∆t c is calculated from At the end of the execution of this algorithm, a list of the mutual collision times is formed. In addition to the fact that the spheres are placed inside an orthogonal computational domain, the collision times with the domain walls are also calculated, ∆t c,w . The next colliding interval is determined as the minimum positive value of both lists: The spheres are allowed to move forward in time by a ∆t c time-interval. The collision time, t c = t + ∆t c , represents the moment of collision between a pair of spheres or a sphere-wall collision. Thus, in the next step the pair collision needs to be resolved, to yield the outcome of the collision and define the starting condition for the next rearrangement stage.

Collision Handling
The sphere-to-sphere collisions, as well as the sphere-to-wall collisions, are very fundamental for this algorithm, as they determine the sequential events of the simulation. The classical elastic collisions type is adopted here. However, this property is not ensured since each sphere is gradually growing and, thus, there is an outward surface velocity relative to the sphere center. A corrective algorithm has been developed in order to treat this peculiarity, described in Section 2.1.5. The following two equations are obtained by the employment of the conservation of momentum and kinetic energy for the two colliding spheres: where ∆ → p m is the momentum change of each colliding sphere, j are the post-collisional velocities of spheres i and j andn ij is the unit vector directed from the center of sphere i to the center of the sphere j. However,n ij ·n ij = 1, thus Equation (7) is reformulated into Equation (8): Assuming that all spheres in the simulation domain are identified with the same unity mass although they may have different volume, the expression for the momentum change reads: In this way, all spheres are allowed to experience the same domain volume irrespective of their size. At the end of each collision step, all spheres have their radii increased by the b rate factor. This will inevitably create overlaps between the spheres, as well as between the spheres and the walls. In this study, b rate assumes a value of 1.00005 that is translated to an increase of 0.005% in all sphere radii after each collision event. This value was found to be optimal for the present simulations. Very small values of b rate would result in significantly elevated running time without any noticeable change in the structure features. On the other hand, simulations with larger b rate values would result in encapsulation of small spheres by larger ones, an unfavorable situation for the algorithm. The reason behind those overlaps lies in the fact that, when the spheres are colliding and their post-collision velocities are calculated, their radius is incrementally increased, creating an artificial overlap, as shown in Figure 1a,b, which is inevitable especially when a jammed state is reached while the spheres resume their actual, desired size. shown in Figure 1a,b, which is inevitable especially when a jammed state is reached while the spheres resume their actual, desired size. In order to resolve those overlaps, a local repositioning algorithm is devised to relocate each sphere to a new position and reduce the overlap. In the present work, a modification of the rearrangement algorithm by Wu et al. [25] is proposed. The spheres are moved sequentially in order to reduce the total overlapping volume gradually. Only one sphere may be moved in each attempt, while the rest of the spheres are assumed immobilized. To this end, an optimization routine is employed that searches for a new position for the sphere under consideration within a distance of twice the sphere radius from the original position of the sphere. The new position should minimize its overlapping volume with the surrounding spheres. The overlapping volume between a pair of spheres is a non-linear function of their relative coordinates, Equation (10): Thus, the total overlapping volume of a sphere with the rest of the N spheres is given from Equation (11): This process aims to the minimization of the total overlapping volume, treated as an objective function, over a set of decision variables under a set of constraints. The minimum-overlap problem can thus be formulated as follows: In order to resolve those overlaps, a local repositioning algorithm is devised to relocate each sphere to a new position and reduce the overlap. In the present work, a modification of the rearrangement algorithm by Wu et al. [25] is proposed. The spheres are moved sequentially in order to reduce the total overlapping volume gradually. Only one sphere may be moved in each attempt, while the rest of the spheres are assumed immobilized. To this end, an optimization routine is employed that searches for a new position for the sphere under consideration within a distance of twice the sphere radius from the original position of the sphere. The new position should minimize its overlapping volume with the surrounding spheres. The overlapping volume between a pair of spheres is a non-linear function of their relative coordinates, Equation (10): Thus, the total overlapping volume of a sphere with the rest of the N spheres is given from Equation (11): This process aims to the minimization of the total overlapping volume, treated as an objective function, over a set of decision variables under a set of constraints. The minimum-overlap problem can thus be formulated as follows: min This formulation transforms the problem of sphere-to-sphere overlapping to a problem of coordinates. Thus, the query to be addressed is now converted to the following-what coordinates of the center of sphere i will minimize its overlap with the rest of the spheres j. The optimizer that is used to solve the minimization problem is a Sequential Least-Squares Programming algorithm (SLSQP) [28] that was developed based on the least-squares solver of Lawson and Hanson [29]. The minimization of the overlapping volume concludes either when a predefined threshold fraction of the overall domain volume is reached or when no further overlap reduction is observed. The threshold fraction is set to 10 −6 , as a compromise between speed and accuracy. Smaller threshold values have also been tested but they were found to increase the computational burden substantially without improving the accuracy in practice. The algorithm optimizes successive second-order (quadratic/least-squares) approximations of the objective function via updates, with first-order approximations of the constraints. Moreover, this method was applied to aerodynamic and robotic trajectory optimization by Kraft [28].

Selection of Sphere to Relocate
The selection of the sphere to be relocated is made with a probability that is proportional to the count of the overlaps with neighboring spheres. Hence, spheres that contribute greatly to the overlapping phenomenon are prioritized in the relocating process. To this end, a weight is assigned to each sphere, indicating the overlap count. Repeating this process for all the spheres in the simulation domain provides a distribution of the overlap count that all spheres exhibit with their neighboring spheres, Figure 2. (SLSQP) [28] that was developed based on the least-squares solver of Lawson and Hanson [29]. The minimization of the overlapping volume concludes either when a predefined threshold fraction of the overall domain volume is reached or when no further overlap reduction is observed. The threshold fraction is set to 10 −6 , as a compromise between speed and accuracy. Smaller threshold values have also been tested but they were found to increase the computational burden substantially without improving the accuracy in practice. The algorithm optimizes successive second-order (quadratic/least-squares) approximations of the objective function via updates, with first-order approximations of the constraints. Moreover, this method was applied to aerodynamic and robotic trajectory optimization by Kraft [28].

Selection of Sphere to Relocate
The selection of the sphere to be relocated is made with a probability that is proportional to the count of the overlaps with neighboring spheres. Hence, spheres that contribute greatly to the overlapping phenomenon are prioritized in the relocating process. To this end, a weight is assigned to each sphere, indicating the overlap count. Repeating this process for all the spheres in the simulation domain provides a distribution of the overlap count that all spheres exhibit with their neighboring spheres, Figure 2.

Generation of the Micro-Voids
Observation of several SEM images reveals that the TI2AlC foams produced with the gelcasting method are comprised, typically, by nearly spherical cells connected by circular windows [30]. In the current approach, the pore windows are simulated as circular cutouts on the spherical cells. Such small circles can be created either by intersecting a plane with the pore cells or with spheres whose centers lie on the surface of the cell. In the present algorithm, the pore windows are generated by creating an arrangement of non-overlapping small spheres on the surface of each pore cell. Alternatively, this can be portrayed as the packing of different caps on the pore cell surface. This approach allows one to translate the problem of pore windows generation into a case of nonoverlapping packing of small spheres on a large sphere surface, as described next.
2.1.5. Non-Overlapping Packing of Spheres on a Sphere Surface The pore windows in each pore cell is sequentially generated in the present algorithm. Initially, the number and size of pore windows on each pore cell surface are extracted from the SEM image

Generation of the Micro-Voids
Observation of several SEM images reveals that the TI 2 AlC foams produced with the gelcasting method are comprised, typically, by nearly spherical cells connected by circular windows [30]. In the current approach, the pore windows are simulated as circular cutouts on the spherical cells. Such small circles can be created either by intersecting a plane with the pore cells or with spheres whose centers lie on the surface of the cell. In the present algorithm, the pore windows are generated by creating an arrangement of non-overlapping small spheres on the surface of each pore cell. Alternatively, this can be portrayed as the packing of different caps on the pore cell surface. This approach allows one to translate the problem of pore windows generation into a case of non-overlapping packing of small spheres on a large sphere surface, as described next.

Non-Overlapping Packing of Spheres on a Sphere Surface
The pore windows in each pore cell is sequentially generated in the present algorithm. Initially, the number and size of pore windows on each pore cell surface are extracted from the SEM image and their frequency distributions are constructed. During the reconstruction process, the number and the size of the pore windows are sampled from this distribution for the specific pore cell size at hand. Initially, the sphere arrangement representing the pore windows on the surface of pore cells takes place using a random process. By iterating and minimizing the overlaps of the spheres on the surface of the macrovoids, an overlap-free arrangement of spheres on the surface of the cell is created. This is performed by initially depositing randomly the sphere centers ensuring a uniform initial distribution on the surface of each pore cell. The minimization of the overlapping function is performed using the differential evolution algorithm [31].
However, there is a significant probability that the resulting macro-voids overlap with each other, creating extra spherical caps that would create extra pore windows connecting the macrovoids, a process that must be included in the reconstruction process. In order to generate pore cells that would contain the exact number and size of the pore windows extracted from the SEM, the radius of the spherical cap that relates to the overlapping of the pair of macrovoids is computed. Then, from each overlapping macrovoid the sampled pore windows sizes to be deposited on their surface are modified. The modification process includes the assignment of the radius of the spherical cap that is related to the pair of the overlapped macrovoids to the pore window with the closest radius to it. Moreover, the position of the sphere-representing the pore window connecting two macrovoids is fixed according to where R mv,1 , R mv,2 are the radii of the overlapping macrovoids and → x mv,1 , → x mv,2 the positions of the centers of the two overlapping macrovoids.

Thickness of the Macrovoids
An arrangement of permeable hollow spheres is obtained after the last step. At the end, the thickness of all hollow spheres is adjusted to achieve the desired foam porosity. Alternatively, the shell thickness can be extracted from the SEM image and used as input to the algorithm. In this case the porosity will be estimated directly from the algorithm bypassing the need for its experimental determination. The capability of the proposed algorithm to offer a reliable estimate of the foam porosity is discussed in Section 3.3. The effect of pore wall thickness on the foam reconstruction is portrayed in Figure 3a,b.
Materials 2019, 12, x FOR PEER REVIEW 9 of 24 distribution on the surface of each pore cell. The minimization of the overlapping function is performed using the differential evolution algorithm [31]. However, there is a significant probability that the resulting macro-voids overlap with each other, creating extra spherical caps that would create extra pore windows connecting the macrovoids, a process that must be included in the reconstruction process. In order to generate pore cells that would contain the exact number and size of the pore windows extracted from the SEM, the radius of the spherical cap that relates to the overlapping of the pair of macrovoids is computed. Then, from each overlapping macrovoid the sampled pore windows sizes to be deposited on their surface are modified. The modification process includes the assignment of the radius of the spherical cap that is related to the pair of the overlapped macrovoids to the pore window with the closest radius to it. Moreover, the position of the sphere-representing the pore window connecting two macrovoids is fixed according to Equation (14) mv,1 mv,2 mv,2 mv,1 where Rmv,1, Rmv,2 are the radii of the overlapping macrovoids and x ⃗ mv,1 , x ⃗ mv,2 the positions of the centers of the two overlapping macrovoids.

Thickness of the Macrovoids
An arrangement of permeable hollow spheres is obtained after the last step. At the end, the thickness of all hollow spheres is adjusted to achieve the desired foam porosity. Alternatively, the shell thickness can be extracted from the SEM image and used as input to the algorithm. In this case the porosity will be estimated directly from the algorithm bypassing the need for its experimental determination. The capability of the proposed algorithm to offer a reliable estimate of the foam porosity is discussed in Section 3.3. The effect of pore wall thickness on the foam reconstruction is portrayed in Figure 3a The flowchart of the reconstruction algorithm proposed in this work is presented in Figure 4. The algorithm begins by extracting statistical information from the SEM image of the porous media and its execution provides the user with a 3D digital reconstruction of the porous media under

Flowchart of the Algorithm
The flowchart of the reconstruction algorithm proposed in this work is presented in Figure 4. The algorithm begins by extracting statistical information from the SEM image of the porous media and its execution provides the user with a 3D digital reconstruction of the porous media under consideration. The algorithm proceeds by performing three main sequential steps, namely (i) the creation of a packing of hollow spheres representing the pores, where the pair overlaps are minimized, (ii) generation of pore windows achieved by the deposition of non-overlapping spheres with traces on the surface of each hollow sphere and (iii) adjustment of the thickness of the pore walls to match the porosity of the real sample.

Pressure Drop
In recent years, a number of studies have shown the importance of determining the various flow regimes before evaluating the Darcy permeability, as well as the inertia correction coefficient [32]. Depending on the flow regime, namely, Darcy, inertial or turbulent, different macroscopic equations must be used and different flow coefficients must be evaluated. A common practice that is followed in the literature to define the flow regime is the graphical method, which is also adopted here and is briefly described in Section 2.2.5.

Flow Simulation
Since the present work focuses mainly on the resolution of flow inside the porous medium, employment of a detailed description of the flow inside the porous medium at pore-level was chosen. The flow characteristics are discussed in Section 3.1.3 using the characteristic numbers for the specific foam structures. Because of the very low values of the Knudsen number, well below the threshold limit of 0.01, the continuum approach was adopted. In the continuum flow regime the thermodynamic equilibrium near the walls and the no-slip boundary conditions are inherently valid [11]. Since the Reynolds number was on the order of unity, the laminar regime was also adopted for the viscous description and the momentum and mass conservations were described by the Navier-Stokes equations. Gas compressibility effects were also non-critical due to the very low Mach number of the systems. In the work by Akbarnejad et al. [33], the effect of fluid flow on the resulting pressure gradient was studied during permeametry experiments on unsealed ceramic foam filters. They empirically derived Darcy and non-Darcy permeability coefficients based on Forchheimer's equation. Their simulations were performed in the turbulence flow regime. Flow in the porous domain was described with a Reynolds-Averaged Brinkman-Forchheimer equation without resorting to pore-level description.
OpenFOAM [34,35] was used here to simulate the gas flow in the reconstructed porous domains. OpenFOAM is an open source Computational Fluid Dynamics (CFD) software package that solves complex fluid flows using the Finite Volumes method and handles chemical reactions, turbulence and heat transfer. It uses various CFD algorithms for laminar or turbulent flow modeling. The rhoSimpleFoam [35] has been utilized, which is a steady-state solver used to simulate compressible Reynolds-Averaged Navier-Stokes (RANS) flows employing a pressure-based (predictor-corrector) approach and the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm [36,37]. This solver offers an evaluation of the mass flow and the temperature variation throughout the working domain. In this context, the conservative, differential form of the 3D compressible Navier-Stokes (N-S) equations is solved. The N-S equations regulate the preservation of the compressible fluid flow mass, momentum and energy and are shown in Equations (15)-(17): where U represents any Cartesian component of the local fluid velocity, p stands for the pressure of a fluid element, ρ is the density and E = e + |U| 2 /2 indicates the total specific energy with e the internal specific energy. The closure of the set of partial differential equations is accomplished with three common assumptions regarding the thermodynamic and kinetic nature of the gas, which correspond to the ideal gas assumption, the calorically perfect gas assumption and the Newtonian fluid assumption. The first two expressions are represented by Equations (18) and (19): where c p and c v are the specific heat of the gas at constant pressure and constant volume, respectively. The heat flux used in Equation (17) is calculated considering the Fourier law, expressed in Equation (20): where µ represents the dynamic viscosity of the gas, α is its thermal conductivity and Pr is the Prandtl number defined for laminar flows. The components of the viscous term tensor are given by where I is the unit tensor.

Transport Properties
In order to solve the energy equation, definition of the fluid properties is necessary, namely, the viscosity, µ, the specific heat, c p and the thermal diffusivity, α. For the sake of comparison with experimental data, the present simulations address the isothermal flow of Argon gas close to the ambient pressure, at T = 279 K in Gelcast93 and Gelcast87 samples and at T = 973 K for the case of the Gelcast84 sample. These properties for the flowing Argon gas are assumed to be constant. Therefore, the following values for the dynamic viscosity, the heat capacity and the Prandtl number are used-µ 279 K = 2.24 × 10 −5 (kg m −1 s −1 ), c p,279 K = 520 (J kg −1 K −1 ), Pr 279 K = 0.66 [38,39], µ 973 K = 5.44 × 10 −5 (kg m −1 s −1 ) [40] c p,973 K = 520 (J kg −1 K −1 ), Pr 973 K = 0.66 [38].

Boundary Conditions
The boundary conditions complementing the transport equations are reported next.

•
Inlet-a flat velocity profile is imposed, corresponding to a fixed volumetric flowrate entering the domain equal to the flowrate reported in the experiments. Argon is assumed to enter the porous domain at a constant static temperature of T 0 = 297 K for Gelcast93 and Gelcast87 and T 0 = 973 K for Gelcast84. • Walls-since the pore walls are assumed to be immobile, every velocity component at the walls is set to zero, as in no-slip conditions. The zero gradient boundary condition is specified for the pressure and the walls are assumed to have the same temperature as their environment, thus a constant static temperature is prescribed on them. • Outlet-since the foam is open to the atmosphere and no recirculation is expected, a Neumann boundary condition is used for the velocity there. A uniform fixed value of 1 atm is assumed for the pressure. The zero gradient boundary condition is imposed on the temperature. • Domain faces parallel to the flow direction-the normal derivatives of the velocity, pressure and temperature are set to zero.

Graphical Estimation
In the graphical method the reduced pressure gradient is plotted against the superficial velocity and the form of the curve reveals the flow regime of interest. In the Darcy regime, the pressure drop is prescribed by the linear Darcy equation, shown in Equation (22): where ∆P is the macroscopic pressure drop and L is the medium thickness. However, in the strong inertia regime the quadratic Forchheimer equation is known to hold, shown in Equation (23): It should be noted that the quantity K F−D is usually referred to as Darcian permeability; however, it is not the same as the Darcy permeability, K D , used in Equation (22).

Representative Elementary Volume (REV) Analysis
A Representative Elementary Volume (REV) analysis has been conducted for each reconstructed foam, namely for the Gelcast84, Gelcast87 and Gelcast93 foams, by calculating the porosity using gradual increments in the REV volume. Every control-volume was cubic in shape and different control volumes were obtained by increasing the edge length. Placing this variable control volume randomly in the reconstructed structure and sampling the porosity of the portion of the structure that is enclosed in the corresponding control volume, offered an efficient correlation between the porosity and the volume of the REV.

Results and Discussion
The successful reproduction of the actual structure is essential to the efficient and realistic description of flow phenomena in the computer-generated samples. A comparison between the geometric properties of the computer-reconstructed structure and those of the experimental data is carried out by visually inspecting and digitally characterizing the reconstructed structures. The suitability of the proposed method, covering a broad range of void fractions between 0.84 and 0.93 and a wide range of pore sizes, was examined for the three foam samples. The ultimate test of the reconstruction efficiency is certainly the comparison of the flow properties with experimental data, as discussed later in this section, given the strong dependence of the flow field and the resulting permeability evaluation on the internal foam structure, especially pore geometry, pore sizes and connectivity.

Structural Features of Reconstructed Foams
Qualitatively, the reconstructed images show features similar to the corresponding ones of their reference images, with good connectivity of the void and the solid phases. In a similar manner to their reference images, the reconstructed gelcast foams contain spherical cells with small circular open pore windows interconnecting neighboring cells, as shown in Figure 5a-c, representing the Gelcast93, the Gelcast87 and the Gelcast84 foam, respectively. The spherical cells were a result of the foaming process, as air bubbles were surrounded by Ti 2 AlC particles. The circular pore windows were attributed either to the contact points between bubbles formed in the liquid during the foaming process or due to the fact that the Ti 2 AlC particles were not fully covering the surface of the bubbles during the foaming process. This resulted into uncovered bubble surfaces that were not present during the sintering process [41]. Thin slices of the reconstructed domains reveal that the Gelcast93 was comprised of large spherical cells, with an increased density of pore windows on each cell existed, amounting to an average of 8 pore windows per pore cell, Figure 6a. The Gelcast87 foam, shown in Figure 6b, was also comprised of spherical cells, yet much smaller than in the Gelcast93 foam. Nonetheless, the pore windows density in each cell remained similar. In contrast, the pore cells were much smaller in the Gelcast84 foam, with also much smaller pore window density in each cell, reaching an average of 3.2 pore windows per cell, Figure 6c. The surface cell density is similar to the cell surface density extracted from the SEM image-the Gelcast93, Gelcast87 and Gelcast84 are characterized by a surface cell density of 1.5 cells/mm 2 , 4.0 cells/mm 2 and 8.5 cells/mm 2 respectively.   Figure 7b,c show that the REV can have any size between the volume of 3.8 mm 3 (minimum REV) and 4 mm 3 , the total sample size for the Gelcast84 and Gelcast87 foams. However, in the case of the Gelcast93 foam, a sample size greater than 30 mm 3 was required, Figure  7a. Thin slices of the reconstructed domains reveal that the Gelcast93 was comprised of large spherical cells, with an increased density of pore windows on each cell existed, amounting to an average of 8 pore windows per pore cell, Figure 6a. The Gelcast87 foam, shown in Figure 6b, was also comprised of spherical cells, yet much smaller than in the Gelcast93 foam. Nonetheless, the pore windows density in each cell remained similar. In contrast, the pore cells were much smaller in the Gelcast84 foam, with also much smaller pore window density in each cell, reaching an average of 3.2 pore windows per cell, Figure 6c. The surface cell density is similar to the cell surface density extracted from the SEM image-the Gelcast93, Gelcast87 and Gelcast84 are characterized by a surface cell density of 1.5 cells/mm 2 , 4.0 cells/mm 2 and 8.5 cells/mm 2 respectively. Thin slices of the reconstructed domains reveal that the Gelcast93 was comprised of large spherical cells, with an increased density of pore windows on each cell existed, amounting to an average of 8 pore windows per pore cell, Figure 6a. The Gelcast87 foam, shown in Figure 6b, was also comprised of spherical cells, yet much smaller than in the Gelcast93 foam. Nonetheless, the pore windows density in each cell remained similar. In contrast, the pore cells were much smaller in the Gelcast84 foam, with also much smaller pore window density in each cell, reaching an average of 3.2 pore windows per cell, Figure 6c. The surface cell density is similar to the cell surface density extracted from the SEM image-the Gelcast93, Gelcast87 and Gelcast84 are characterized by a surface cell density of 1.5 cells/mm 2 , 4.0 cells/mm 2 and 8.5 cells/mm 2 respectively.   Figure 7b,c show that the REV can have any size between the volume of 3.8 mm 3 (minimum REV) and 4 mm 3 , the total sample size for the Gelcast84 and Gelcast87 foams. However, in the case of the Gelcast93 foam, a sample size greater than 30 mm 3 was required, Figure  7a.

Pore Size Distribution of Reconstructed Foams
Pores are defined by utilizing the generalized network extraction algorithm developed by Raeini et al. and by assuming a spherical shape to simplify the analysis [42]. The pore size distributions of the Gelcast84, Gelcast87 and Gelcast93 foams are shown in Figure 8a-c, respectively and a statistical description in view of the minimum, the maximum, the median and the mean values of their pore size distributions is given in Table 1.

Pore Size Distribution of Reconstructed Foams
Pores are defined by utilizing the generalized network extraction algorithm developed by Raeini et al. and by assuming a spherical shape to simplify the analysis [42]. The pore size distributions of the Gelcast84, Gelcast87 and Gelcast93 foams are shown in Figure 8a-c, respectively and a statistical description in view of the minimum, the maximum, the median and the mean values of their pore size distributions is given in Table 1.

Pore Size Distribution of Reconstructed Foams
Pores are defined by utilizing the generalized network extraction algorithm developed by Raeini et al. and by assuming a spherical shape to simplify the analysis [42]. The pore size distributions of the Gelcast84, Gelcast87 and Gelcast93 foams are shown in Figure 8a-c, respectively and a statistical description in view of the minimum, the maximum, the median and the mean values of their pore size distributions is given in Table 1.  It is reported in literature that ceramic foams produced by the same gelcasting method are characterized by a bimodal pore size distribution [41], where the large pores are attributed to the foaming process and smaller pores are distributed over the cell walls creating the pore windows, which can be attributable to the contacts among gas bubbles. However, smaller pores can be observed and can be explained by the fact that the powders used for the gelcast generation are not in strict touch with each other; they are distributed over the bubble surface. These powders present a point of contact which tends to be extended during the sintering process, yet it does not fill the entire space among the spheres [32]. The pore size distributions that are obtained from the reconstruction follow a bimodal distribution as shown in Figure 8a-c. The larger part of the macro porosity is attributed to the cells; however, a micro porosity is also present, attributed to the smaller pore windows.
Using the assumption that the macroporosity and the microporosity are independent of each other, the total pore size distribution is regarded as a mixture of the individual pore-size distributions [43] where rA is a random variable designating the pore radius, λ1 is the probability density for a micropore with size rA and Pmicro and Pmacro are the micro and macro pore contributions, respectively. The probability densities of the micro and the macro pores, respectively, with a pore radius rp, follow the Gaussian distribution, where μ1,p and sd1,p are the mean and standard deviation of the pore radius attributed to the microporosity and μ2,p and sd2,p are the mean and standard deviation of the pore radius attributed to the macroporosity. The values for the three foams are given in Table 2.  It is reported in literature that ceramic foams produced by the same gelcasting method are characterized by a bimodal pore size distribution [41], where the large pores are attributed to the foaming process and smaller pores are distributed over the cell walls creating the pore windows, which can be attributable to the contacts among gas bubbles. However, smaller pores can be observed and can be explained by the fact that the powders used for the gelcast generation are not in strict touch with each other; they are distributed over the bubble surface. These powders present a point of contact which tends to be extended during the sintering process, yet it does not fill the entire space among the spheres [32]. The pore size distributions that are obtained from the reconstruction follow a bimodal distribution as shown in Figure 8a-c. The larger part of the macro porosity is attributed to the cells; however, a micro porosity is also present, attributed to the smaller pore windows.
Using the assumption that the macroporosity and the microporosity are independent of each other, the total pore size distribution is regarded as a mixture of the individual pore-size distributions [43] where r A is a random variable designating the pore radius, λ 1 is the probability density for a micropore with size r A and P micro and P macro are the micro and macro pore contributions, respectively. The probability densities of the micro and the macro pores, respectively, with a pore radius r p , follow the Gaussian distribution, where µ 1,p and s d1,p are the mean and standard deviation of the pore radius attributed to the microporosity and µ 2,p and s d2,p are the mean and standard deviation of the pore radius attributed to the macroporosity. The values for the three foams are given in Table 2. Table 2. Mean and standard deviation of the microporosity and macroporosity of the foams examined and the probability density for the micropores.

Characteristic Flow Numbers
For the determination of the dimensional numbers that define the flow regime a strict definition of the characteristic length is required. In this research the hydraulic radius is replaced with the characteristic length of the cross-sectional area of the flow. For one-dimensional flows, the ratio of the cross-sectional area to the wetted perimeter is usually defined as the characteristic length of the flow. By projecting this concept to three dimensions, the characteristic length is defined as the ratio of the total interstitial space to the total interphase surface between the solid and the void phases [44], which is translated to the inverse of the specific surface area of the porous medium. The specific surface area, the characteristic length, the Knudsen, Reynolds and Mach numbers of each foam sample are shown in Table 3. The characteristic length of each foam is similar to its mean pore size shown in Table 1, showing that the utilization of either one as characteristic length would not affect considerably the estimation of the characteristic numbers of the flow. Since the Knudsen number (Kn) is smaller than the continuum threshold, the validity of Navier-Stokes equation is assured, whereas, since Re is of the order of unity, the laminar flow regime description is valid. Mach numbers (Ma) much lower than 0.2 also reveal that compressibility effects are practically not significant for the current flow description. Table 3. Specific surface area and resulting characteristic numbers of the foam simulations.

Domain Meshing
An initial simple hexahedral mesh was created with the blockMesh utility of OpenFOAM. The initial mesh produced by the blockMesh utility was fed in conjunction with the geometries to the snappyHexMesh utility, in order to generate the final hexahedral dominant mesh. The snappyHexMesh utility is an OpenFOAM subprogram that acts simultaneously as a mesh sculptor and a mesh generator [45]. It takes an existing simple mesh like the one generated by the blockMesh and fits the mesh to complex geometries by simultaneously performing mesh refinement and boundary layer addition. The most critical operation was the construction of a high-quality mesh to ensure grid-independent results. To this end, a range of global refinements were tracked throughout the domain to test quantitatively the resulting mesh suitability with the CFD code. The global refinements were performed on the initial hexahedral mesh (using blockMesh) and a range of mesh densities from 25 up to 200 points per dimension was obtained and examined for the mesh independency on the foam with the smaller pore size, namely on the Gelcast84. Local refinement and boundary layers was allowed (using the snappyHexMesh) on each mesh with the same parameters for each resolution. In order to avoid significant numerical errors and extensive use of multiple iterations of non-orthogonal correctors [46], the maximum cell orthogonality in all meshes was set to 65 degrees and the maximum cell skewness was set to 4. The obtained meshes were then tested against their accuracy with respect to the corresponding pressure drop, Figure 9b. The mesh independence study indicated that a mesh with the domain edge discretized by 150 elements, shown in Figure 9a, was sufficient to capture the essential behavior of the model with suitable accuracy. A comparison of the averaged pressure drop, for a superficial velocity of 0.09 m/s, for the Gelcast84 foam and for different discretization densities is shown in Figure 9b. The results indicate that for 150 elements per edge length of the domain are suitable for these calculations, since larger discretization densities introduce a variation of less than 5%. The most suitable meshes were obtained by the grid independence analysis and their statistics for each foam type are shown in Table  4: Table 4. Summary of the cell type distribution of the three different foams. Gelcast93  4603068  3395049  150086  19190  1038743  Gelcast87  5576328  4106789  212812  5656  1251071  Gelcast84  6232131  4532497  293359 5816 1400459  The mesh independence study indicated that a mesh with the domain edge discretized by 150 elements, shown in Figure 9a, was sufficient to capture the essential behavior of the model with suitable accuracy. A comparison of the averaged pressure drop, for a superficial velocity of 0.09 m/s, for the Gelcast84 foam and for different discretization densities is shown in Figure 9b. The results indicate that for 150 elements per edge length of the domain are suitable for these calculations, since larger discretization densities introduce a variation of less than 5%. The most suitable meshes were obtained by the grid independence analysis and their statistics for each foam type are shown in Table 4:  Figure 10a,b illustrates the velocity magnitude distribution of the domain on a through-plane at the middle of the structure, reflecting the flow of air through the spherical pores and the corresponding pore windows. As expected, the bulk of the flow is channeled through high-void space regions, whereas the flow is severely restricted in other areas.

Permeability Estimation
In total, four realizations of each foam were produced and the pressure gradient profiles as functions of the mean fluid superficial velocity were obtained and averaged. Figure 11 shows the interdependence of the pressure drop per unit length with the superficial velocity for all three foam samples examined here. A comparison of simulated results with experimental data is made in the figure that reveals a very satisfactory agreement [8,9]. To be more specific, the quality of the fit of the simulated to the experimental data was measured in terms of a normalized root mean square deviation (NRME) and revealed that the NRME value is 0.028 for the Gelcast93 sample, 0.026 for the Gelcast87 and 0.017 for the Gelcast84. Moreover, the variance introduced by the randomness of the foam realizations shows that Gelcast93 exhibits an average deviation of 6%, while Gelcast87 and Gelcast84 exhibits an average deviation of 8% and 4%, respectively. These values are practically very low in view of the complexity of the foam structures. Figure 10a,b illustrates the velocity magnitude distribution of the domain on a through-plane at the middle of the structure, reflecting the flow of air through the spherical pores and the corresponding pore windows. As expected, the bulk of the flow is channeled through high-void space regions, whereas the flow is severely restricted in other areas.   In total, four realizations of each foam were produced and the pressure gradient profiles as functions of the mean fluid superficial velocity were obtained and averaged. Figure 11 shows the interdependence of the pressure drop per unit length with the superficial velocity for all three foam samples examined here. A comparison of simulated results with experimental data is made in the figure that reveals a very satisfactory agreement [8,9]. To be more specific, the quality of the fit of the simulated to the experimental data was measured in terms of a normalized root mean square deviation (NRME) and revealed that the NRME value is 0.028 for the Gelcast93 sample, 0.026 for the Gelcast87 and 0.017 for the Gelcast84. Moreover, the variance introduced by the randomness of the foam realizations shows that Gelcast93 exhibits an average deviation of 6%, while Gelcast87 and Gelcast84 exhibits an average deviation of 8% and 4%, respectively. These values are practically very low in view of the complexity of the foam structures.
In this study, the Forchheimer coefficients KF-D and KF-F were obtained by fitting the pressure drop against fluid velocity curves with polynomials of the second order on velocity, with a regression coefficient higher than 0.999 in all cases. The only fitting parameter was the thickness of the pore cells, as discussed above. Values of the simulation and experimental data for the permeability are shown in Table 4 for the sake of direct comparison with each other. The resulting Forchheimer permeability and the estimated inertia coefficient are very close to the experimentally determined values considering that the entire simulation was based on a single SEM image of a small sample of the actual foam.
Moreover, the Darcy permeability and the tortuosity for the corresponding foams were calculated, Table 5. The Stokes flow equation was solved for the determination of the Darcy permeability. The tortuosity was determined by simulating the diffusion of non-sorbing species inside the pore space. More specifically, the diffusion process was simulated by performing discrete lattice walks in the void space [47] and the mean-square displacement of the non-sorbing walkers provided estimation for the diffusive tortuosity. Care was taken in order to ensure that the random walks were conducted under the bulk diffusion regime. This was accomplished by setting the mean Figure 11. Pressure drop as expressed in References [8,9] vs. superficial velocity for the (a) Gelcast84 foam, (b) Gelcast87, (c) Gelcast93 foams. Comparison of simulation results with experimental data.
In this study, the Forchheimer coefficients K F-D and K F-F were obtained by fitting the pressure drop against fluid velocity curves with polynomials of the second order on velocity, with a regression coefficient higher than 0.999 in all cases. The only fitting parameter was the thickness of the pore cells, as discussed above. Values of the simulation and experimental data for the permeability are shown in Table 4 for the sake of direct comparison with each other. The resulting Forchheimer permeability and the estimated inertia coefficient are very close to the experimentally determined values considering that the entire simulation was based on a single SEM image of a small sample of the actual foam.
Moreover, the Darcy permeability and the tortuosity for the corresponding foams were calculated, Table 5. The Stokes flow equation was solved for the determination of the Darcy permeability. The tortuosity was determined by simulating the diffusion of non-sorbing species inside the pore space. More specifically, the diffusion process was simulated by performing discrete lattice walks in the void space [47] and the mean-square displacement of the non-sorbing walkers provided estimation for the diffusive tortuosity. Care was taken in order to ensure that the random walks were conducted under the bulk diffusion regime. This was accomplished by setting the mean free path of the walkers to a sufficiently low value that would correspond to Knudsen number at the continuum limit (0.01 or less).  The Darcy permeability, K D , deviated from the experimentally determined Forchheimer permeability, K F-D,exp . This observation confirms the conceptual difference between the Darcy and the Forchheimer permeability. The tortuosity values show that the effective diffusivity inside the foams decreases as the foam porosity decreases. Although the Darcy permeability was found to differ by almost two orders of magnitude between the Gelcase93 and Gelcast84, the tortuosity difference was only 32%. This is attributed to the fact that flow is considerably more sensitive to pore size than bulk diffusion.

Pore Wall Thickness
The only fitting parameter of the present algorithm is the wall thickness of the pore cells. However, it was found that the value that was chosen to fit the experimental permeability data was almost identical to the value that could be extracted from the SEM images of all three foam samples, as shown in Table 6. Consequently, the reconstruction method that was developed here could eventually be applied without having to resort to any fitting parameter whatsoever, yet reproducing the experimental data for the flow coefficients of the actual foams shown in Table 7.

Overlap Minimization Aspects
From the technical viewpoint, it must be noted that the SLSQP optimizer was chosen for the minimization of the overlap during the sphere rearrangement process. In theory, SLSQP is appropriate for finding a local minimum of an objective function, however it does not guarantee a solution near a global minimum. Nevertheless, this is not a concern due to the fact that the local minimum obtained by the SLSQP algorithm is still a strong improvement towards the reduction of the overall overlapping. The successive relocation of multiple spheres ultimately reduces the total overlapping volume. Thus, a more complex, global optimization routine to relocate the spheres would not be efficient as it would probably significantly increase the calculation time without an essential improvement in the final result.

Limitations and Future Work
The proposed reconstruction algorithm focuses on the reconstruction of isotropic ceramic foams produced with the gelcast method; however, for a plethora of foamy materials, such as multiscale polylactic acid structures (with largely varied secondary porosity and pore radii), pore shape distribution and pore size vary greatly. Building a precise 3D digital reconstruction using a single 2D slice is rather difficult. This is attributed to the inability to ensure the representative adequacy of the 2D images at hand.
Moreover, in the proposed reconstruction algorithm the shape of the pores, as well as the shape of the pore windows, are assumed to be spherical. This assumption seems to be valid in the majority of the foams produced with the gelcast method, since the pore shapes are nearly spherical. If pore cells strongly deviate from the spherical shape (e.g., foamy layers produced with the phase inversion method), the proposed approach can be extended to an arbitrary pore shape, provided that there is an analytical expression of the overlapping volume between the pore objects.

Conclusions
A 3D reconstruction method was developed in this work, aiming at the stochastic digital representation of the internal structure of ceramic foams that are produced by the gelcast method using a single SEM foam image. The only fitting parameter that was employed was the thickness of the pore walls, while all other geometrical and topological factors were directly matched to the experimentally extracted ones from image processing. The pore wall thickness can be treated as the single fitting parameter in the entire reconstruction process, although it is shown that images of sufficient resolution could eliminate the need even for that. It was also shown that the pore wall thickness was very close to the value that could be extracted from the SEM images of all three types of foams examined here.
The novel hybrid algorithm developed here combines a variation of Lubachevsky-type and Random Close Packing of Hard Spheres algorithms to obtain near-optimum solutions to the pore-sphere packing problem. The proposed method ensures that the pore cell size distribution matches the one extracted from the imaging data, a factor missing from previous tessellation-related works. In addition, the proposed method approaches the curvature of inherently spherical foam cells using spheres rather than deformed polyhedra, as implemented in previous reconstruction methods.
The bimodal porosity distribution of the foams, in view of the macrovoids of the cells and the smaller porosity of the interconnecting pore windows, is captured with the proposed algorithm for all three types of foams examined here, in accord with experimental observations related to their internal structure description. Moreover, a key advantage of the proposed reconstruction algorithm is its ability to reconstruct realistic foam structures with pairs of pore cavities interconnected via multiple pore windows, a feature that cannot be reproduced with typical tessellation processes.
Numerical simulations are performed directly on the 3D reconstructed foams to determine their gas permeability. The model predictions are compared with experimental gas permeability data that were obtained for the physical samples, showing a satisfactory agreement, as they produced quantitatively similar pressure drops with experiments under variable superficial velocity, with a residual deviation, quantified by the normalized root mean square deviation (NRME), in the very low range of 1.7-2.8%. Moreover, the successful repeatability of the reconstruction method is shown in the low value of the deviation introduced by the randomness of the foam realizations, which shows that gelcast foams studied here exhibit a low deviation in the range of 4-8%.
The difference between the first-order Forchheimer term and the Darcy permeability was confirmed by direct calculation of their values for the same computational domains. Despite the complexity of the structure of these foams, it was found that the 3D reconstruction of pore cells and pore windows connecting the cells was successful judging from the calculation of the flow properties of the samples.
Overall, the proposed reconstruction process can assist in a rapid characterization of the foams produced by the gelcast method, avoiding most of the time-consuming experimental characterization processes. With information extracted solely from a few SEM images, such as the pore size and the window size distribution and the pore wall thickness, the user can obtain an accurate 3D digital representation of the foam at hand. The 3D reconstructed volumes can be further utilized in order to obtain valuable structural information, such as the specific surface area, as well as their corresponding effective resistance to mass transport and flow.