Abstract: The quality control of pore sizes and pore distributions are critical to its success in providing an alternative for pore-scale simulations. However, care must also be exercised, while checking the mesh quality, as this also controls accuracy in all numerical simulations. There are limitations to the technique and some remedies to overcome these limitations. To produce reliable solutions at the pore-scale, adequate amount of grid points across the pores are required. In low porosity porous media, the void volume occupied by fluid is much less than in the high porosity ones. Hence, with prescribed uniform mesh spacing, this is equivalent to fewer cells and less resolution available for solving the Navier-Stokes equations. To meet this requirement of sufficient resolution, the technique of adaptive meshing has been incorporated in this study to locally refine meshes further in those regions that require it. This technique is robust in accurately predicting transport in porous media, while accounting for pore-scale phenomena.

INTRODUCTION

Although, the modeling of fluid transport in porous media has been extensively studied in the past, the length scales involved and the complexity of pore structure make such studies mathematically and numerically challenging even up to the present day. The detailed geometry of the pore structures, aside from statistical (or manufacturing) variations, can be given through microtomography or microscopy. Such studies are limited to only small representative samples, giving information for only small portions of perhaps a more complicated porous structure. Even if one knows the detailed structure, a general exact solution describing the mass and heat transfer in these complex pore geometries is difficult to obtain. Most historical and current reported research can be categorized as being either theoretical simplified or approximate, experimental, numerical or combinations of the above. This study adopts a numerical approach. In the numerical modeling approach, most models employed so far treat the porous medium as a single phase continuum, using the volume averaging approximations. For these models, a priori knowledge of spatially averaged transport parameters (Vafai and Tien, 1981; Vafai et al., 1985; Calmidi and Mahajan, 1999) is required. Often, those parameters can only be obtained through experiments. As a result of the cost and time required for experiments, these models are not generally applicable for materials over a wide range of pore composition, shapes and concentrations. In addition, there is a real concern as to whether the averaged parameters at the macroscopic level can really be representative of multiple scales of fluid transport in this type of problem. If the length scale of interest is small by comparison with the scale over which the parameter is averaged (such as occurs in the dispersion effect), then averaged parameters are not suitable for describing any phenomena of smaller scale.

Some other researchers have developed numerical models that involve directly solving the momentum equations in the microstructure of the porous medium. However, these models involve a variety of simplifying assumptions. For example, a bank of parallel cylinders or a bed of parallel spheres in lattices are often used to represent fibrous or granular materials, respectively. Assumptions of lattice periodicity and material homogeneity are common in these models. Some limiting cases of this sort can be solved analytically. For example, Larson and Higdon (1989) studied Stokes flow in a periodic array of cylinders. In these parallel lattice models, fluid flows around a regularly spaced array of separated solid objects, instead of flowing through tortuous interstitial gaps in a solid matrix. Thus, the phenomena of mixing and of recirculation of local fluid (which are classified as secondary flow effects) and the resulting thermal dispersion (Calmidi and Mahajan, 2000) are neglected. In addition, the exclusion of wall effects, variable porosity and heterogeneous materials in this kind of approach limit the applicability of analyses to only special cases. A recently developed network model (Thompson and Fogler, 1997; Thompson, 2002) and the emerging Lattice Boltzmann model have been applied to simulating porous medium transport and their merits in this regard have been demonstrated by Martys (2002) and Verger and Ladd (1999). However, in Lattice Boltzmann models, the lack of precise knowledge of the location of the pore/solid interface is a potential source of error. In addition, the handling of interfaces between solid and fluid in Lattice-Boltzmann model can be tedious and computing intensive due to its rigid lattice structure. In general, these aforementioned numerical efforts involve one or more of the following assumptions, which limit the applicability of the approach to special situations:

Mesh-based microstructure representation algorithm: Some earlier porous medium transport simulation models (Martys, 2002; Verger and Ladd, 1999) have accounted for the microstructure by incorporating both the solid phase and the fluid phase as distinct phases in the geometry of the domain, but they are generally accompanied by simplifying approximations in the transport equations. Typically, a simplified set of transport equations is solved in the fluid regions, possibly together with the heat conduction equation in the solid regions. In the conventional approach, a prescribed physical porous domain of interest is mapped to a computational domain, sometimes by periodic extension of an optical scan of a microtomograph of a representative porous sample, followed by a labeling of solid and fluid regions. This is in itself a difficult task, especially when attempting to numerically represent a three-dimensional porous microstructure. Next, the computational domain is geometrically discretized by mesh generation and finally the governing equations are discretized on the mesh and solved for unknowns at each grid point or cell center by using an appropriate numerical algorithm, subject to corresponding boundary conditions. In the second step of mesh generation, major difficulties arise from the need to respect the geometry of the fluid-solid interface in order to accurately capture the pore-scale transport phenomena. If a structured multiblock mesh is generated, the difficulty of gridding the numerous tortuous pore passages are extreme. Even if an unstructured mesh is generated, the need for the mesh to conform closely to a very complex boundary structure places great demands on the mesh generator. Furthermore, the meshes on the fluid and solid sides of an interface must conform to each other at the interface, or else computationally expensive interpolation will need to be carried out at the interface, often with a resultant loss of accuracy. These difficulties in generating a quality mesh that respects a physically prescribed convoluted fluid-solid interface have in the past made a full Navier-Stokes solution of flow in porous media seem infeasible. Therefore, typically either the microstructure geometry that has been vastly altered so that it can be described by a regular periodic lattice, or only a global parameter such as porosity is employed, which masks the dynamics within the porous region.

MATERIALS AND METHODS

In this study, the difficulties outlined above are circumvented by a new approach to the computational domain definition and mesh generation, which reflects the random nature of the spatial distribution of the pores. In the conventional approach, a prescribed physical domain with existing fluid-solid interface is mapped to a computational domain and the latter is geometrically discretized by mesh generation. In the new approach developed here, the reverse sequence of steps is followed. First, the outline of the computational domain is selected to represent the physical domain of interest. Next, the computational domain is discretized with a fine enough mesh, which is generally an easy task since there are no convoluted interfaces to conform to. Then mesh locations are chosen as seed locations of solid material; these locations are chosen randomly, but with prescribed rules. The solid regions are then allowed to grow, again randomly but in accordance with prescribed rules. The rules are established so that the eventual statistical properties of the porous region matches closely the geometric properties (such as the porosity and the pore size) of the physical porous medium that it is sought to simulate. The solid and fluid-filled pore regions are then assigned the thermodynamic and transport properties of the corresponding physical materials. As a natural byproduct of this process, the interface between solid and pore regions always lies at the boundary between mesh cells. Thus, it results in a complete conforming at the solid and fluid interfaces. Because the mesh cells act as a bed or base from which the porous medium is constructed, this approach is termed a mesh-based microstructure representation algorithm. This algorithm has 3 advantages over the conventional approach in that the mesh generation is an easy task, the meshes automatically conform to the solid-fluid interface and to each other at the interface and no laborious transfer of the solid-fluid interface geometry from the physical to the computational domain is required.

This methodology entails three steps, namely generating the mesh, selecting seeds and growing the solid matrix from seeds. A mesh generator is used to generate a bed of unstructured or structured meshes in the computational domain. This mesh bed is used both for generating the solid matrix and for performing the numerical computations. The mesh used herein is triangular and quadrilateral for two-dimensional and tetrahedral and hexahedral for three-dimensional geometries. Seeds, which are initial solid elements are selected within the mesh bed. The seed locations are chosen using a random number generator, which has been well validated for its uniformity of randomness (Matsumoto and Kurita, 1994). These seed cells are grown spatially into a geometry determined by imposed constraints; for example a fiber-like or sphere-like solid matrix as used in the current computations. The processes of seeding and growing are iterated until the desired pore quality and porosity are achieved. In summary, the procedure is:

RESULTS AND DISCUSSION

It is expected that the more the mesh is refined, the more closely the micro-scale can be simulated. In the vicinity of the solid-fluid interface, there exist large velocity and temperature gradients and an adequate number of mesh points is necessary to accurately resolve these variations. Also, the use of fine grids in the flow direction may be critical when the effects of axial conduction and convection are of concern. However, the desire for finer meshes is limited by available computing resources. Hence, it was necessary to conduct a grid-independence test before this mesh-based methodology could be used for further computations. The test was done by varying mesh sizes while holding porosity to a constant. Table 1 shows the results obtained in this grid dependence study.

Table 1:	Five cases are computed with various combinations of mesh sizes in axial and transverse directions for grid dependence study

The studied case is a 2D channel flow with porous sample located in the middle of the channel. The flow is fully developed before, it reaches the porous section. The inflow, outflow and boundary conditions are kept the same throughout this numerical sensitivity test. A structured mesh was used, with fibrous solid matrix, porosity of 0.83 and Reynolds number R_e of 1093. The first column in Table 1 refers to the spacing in the x and y directions. The second and third columns refer to the total number of grid points in the x and y directions, respectively. The fourth column is the dimensionless pressure drop across the porous medium. The case number is used for discussion.

Grid points are clustered along the y direction in the regions near the bounding walls. The Δy used here is the local length scale of the first grid point away from wall. Both Δx and Δy are in meters. The dimensionless pressure drop, B, is used to measure the required mesh resolution in both x and y directions. It can be concluded that Δy of 10^-4 m is fine enough in the y-direction, i.e., the cross flow direction, if case 1 and 3 are compared. For these 2 cases, Δy varying from 2.5x10^-4-1.9x10^-4 has negligible effect on the computed pressure drop. Since initial pore locations (seeds) are chosen by a random number generator, two sets of initial seeding lead to almost same amount of pressure drop. It shows that flow is insensitive to where seeds are planted so long as the porosity and the shape of solid matrices are kept the same. Hence, Δy of 2.5x10^-4 m, which has been shown to be adequate for obtaining a converged solution, will be used for all computations to minimize the required computing resources. Similarly, the Δx of 2.0x10^-4 m is chosen for the same reasoning.

CONCLUSION

This new modeling technique for porous medium transport not only allows for better quantification of how microscopic properties affect macroscopic transport, but also has been demonstrated to be quite robust. The geometrical shape of the solid phase is determined by the type of mesh of which it is composed. The pore density is currently controlled by random number generator; it can be extended to include defined geometries. The geometry of the numerical pores/granules generated in this study would not be a perfect circle or sphere and would possibly be sharp-edged ligaments in open-cell form or exhibit a bit of zigzagged fibrous shape. Nevertheless, such characteristics mimic those of the real manufactured porous media which come with slightly varying pore size and sometimes even with rough edges (Lu et al., 1998). There are limitations to the technique and some remedies to overcome these limitations. To produce reliable solutions at the pore-scale, adequate amount of grid points across the pores are required. In low porosity porous media, the void volume occupied by fluid is much less than in the high porosity ones. Hence, with prescribed uniform mesh spacing, this is equivalent to fewer cells and less resolution available for solving the Navier-Stokes equations. To meet this requirement of sufficient resolution, the technique of adaptive meshing can be incorporated to locally refine meshes further in those regions that require it. The second limitation of the current technique is that the flow is no longer continuum in nature if the local characteristic length (finest mesh scale) is comparable to the mean free path of the fluid Knudsen number >>1. Though the behavior of the fluid flow near the wall is not yet fully understood, replacement of the non-slip condition with a Maxwellian slip condition at the interface or wall (Arkilic et al., 1997) is a common practice in the range of 10^-3 ≤Kn ≤0.1 (Barber and Emerson, 2002). In conclusion, this technique is preferable in accurately predicting transfer in porous media while, accounting for pore-scale phenomena. It is computationally more expensive than conventional simplified approaches, but has a great potential with the continuously evolving cheaper and faster computing resources.