Flows in porous media
A porous material is a material consisted of solid and voids (cavities, channels or interstices). Clearly, porous media are abundant in nature in forms of rocks and soils. On the other hand, their specific properties (large surface-to-volume area) make them useful for a wide variety of engineering applications, such as fuel cells, chemical reaction surfaces or filters. Hence transport processes need to be studied in them, whether thermal or mechanical.
In this work, we want to demonstrate the CFD simulation in such medium in MATLAB. We start with the medium geometry in form of a set of black-and-white bitmaps, much like those ones after reconstruction of a grayscale Computer-Tomograph (CT) image. In this case, such simulations were performed using the QuickerSim CFD Toolbox for MATLAB while the meshing was done using the Partial Differential Equations (PDE) Toolbox. We additionally used the well-known marching-cubes algorithm in an intermediate step.
Flow simulations – Darcy law vs. DNS
When one thinks of the fluid flow in a porous medium, it is usually the Darcy equation that comes into mind. It is a linear relation between the average flow velocity and the imposed pressure gradient,
Similarly, when the flow velocity gets sufficiently high, the non-linearity in the pressure gradient- velocity relation is observed. Then one speaks of the Darcy-Forchheimer relation, with the so-called viscous and inertial permeabilities,
Still, these two equations are just the macroscopically averaged forms of the momentum equation (Whitaker 1986, 1996). The truth is the fluid flow that is governed by either Navier-Stokes or, if the flow is sufficiently slow and the inertial effects can be neglected, the Stokes equation, at least as long as we deal with the continuous medium. Simply speaking, here we are interested in the direct numerical simulation, DNS, or, this term is sometimes also used, direct pore-level simulation, (Regulski 2015, Parthasarathy 2016).
Geometry and mesh
To perform a simulation in the QuickerSim CFD Toolbox for MATLAB we need to generate the mesh stored in the [p,e,t] format. There are functions in the PDE Toolbox that can create the mesh from the STL-format. However, first we need to create the STL file itself from the binary staircase image of the porous medium. The input geometry is a voxel map of size 40x40x40 pixels with the resolution of 0.1 mm what translates to 4 mm x 4 mm x 4 mm of the entire specimen.
We created a separate tutorial that actually performs that, see Tutorial 28**. What it basically does is to read a set of binary images, stack the binary data in a three-dimensional array and create the surface composed of triangles in the location where the solid turns into the void. Example of such geometry is shown in Figure 1 – we show both the solid skeleton of the structure as well as the void zone where the flow actually occurs.
Note that these both can be generated by exchanging zeros and ones in the array storing the pixel data. Hence if you need to perform e.g. a heat transfer or strength analysis for the solid part of the medium, you can get the mesh easily from here!
Figure 1. The solid skeleton of the structure and the void zone where the flow actually occurs.
**the mesh generation will work only with the PDE Toolbox. However, you can still generate and visualize the STL file if you do not have the PDE Toolbox. You can also get the meshes with the latest release of our Toolbox from < here >.
Flow simulations
The objective of our simulation was to check a few issues. First was just the proper mesh size for our endeavor, then the focus switched to the proper inlet boundary condition. Finally, we wanted to check what the relative difference between the linear and inertial flow actually is.
Hence three scenarios were investigated depending on the numerical model (Stokes/Navier-Stokes) and inlet velocity profile:
Scenario I – Navier-Stokes equations were being solved and constant velocity profile was assumed. Here we focused on the mesh quality.
Scenario II – Navier-Stokes equations were being solved and nonhomogeneous velocity profile was assumed.
Scenario III – Stokes equations were being solved and nonhomogeneous velocity profile was assumed.
You can check out the single flow simulation across such a porous foam in our Tutorial 29.
The fluid flowing through the medium was assumed to be water (with kinematic viscosity 
). Three types of boundary conditions were imposed: inlet with a prescribed velocity profile, outflow (constant pressure) and wall (no-slip condition) on the material walls and on the sides.
Scenario I
Three grids of different resolutions were generated. Basic parameters of each mesh are given in Table 1 below while the grids are shown in Figure 2.
| Grid | Number of points | Number of elements |
| Coarse | 19218 | 82439 |
| Medium | 31590 | 142732 |
| Fine | 43753 | 202202 |
Table 1. Coarse, medium and fine meshes. They are illustrated in Figure 2.
Figure 2. Computational grids with different level of refinement
For each grid, five simulations with different inlet velocity magnitude were performed. The chosen velocity values were 0.025, 0.05, 0.075, 0.1 0.125 m/s. That translates to the Reynolds number of 100, 200, 300, 400 and 500 based on the specimen diameter. Comparison of pressure drop for all cases is shown in Figure 3.
Figure 3. Pressure drops obtained at varying mesh size. Points denote the velocities at which simulations were run.
The discrepancy in pressure drop at the largest velocity between the coarsest mesh and the finest mesh is at the level of 5% while the medium and finest differ around 1%. Also note that all points follow a paraboloid pattern. That indicates that we are in the inertial flow regime and the Darcy-Forchheimer formula should apply. For the finest grid the parabola was fitted using the least squares fit method. It is presented in Figure 4 below.
Figure 4.
Additionally, pressure and x-velocity fields are presented for the entire fluid domain in Figure 5 and in the cross-sectional plane (perpendicular to the y-axis) in Figure 6.
Figure 5. Velocity and pressure distribution.
Figure 6. Velocity and pressure fields in the cross-sectional plane.
Scenario II
Note that for all previous cases an assumption of a uniform distribution of velocity across the inlet plane was taken. Since the no-slip condition on the walls was applied, discontinuities in the velocity field arose, which is non-physical. Hence the nonhomogeneous inlet velocity profile needs to be imposed. For some cases, analytical expressions for the shape of such profiles are developed (see Tutorial 23). For the current simulation velocity magnitude for each inlet, the node was calculated based on quadratic distance from the nearest wall (for more details see Tutorial 29). The profile obtained from this approach is shown in Figure 7.
Figure 7. Inlet velocity profile.
Comparison of pressure drops for both cases for the same values of velocity as in Scenario I is presented in Figure 8. Velocity in case of nonhomogeneous profile means average value velocity across the inlet plane.
Figure 8. Pressure drop obtained in the case of homogeneous and variable inlet velocity profile
Clearly, the discrepancy between pressure drops is negligible. However, the differences between pressure and velocity distributions are significant (see Figures 9,10 and Figures 11,12).
Figure 9. Velocity and pressure distribution in case of homogeneous velocity profile at the inlet.
Figure 10. Velocity and pressure distribution in case of non-homogeneous velocity profile at the inlet.
Figure 11. Velocity and pressure distribution in the cross-sectional plane in case of homogeneous velocity profile at the inlet.
Figure 12. Velocity and pressure distribution in the cross-sectional plane in case of non-homogeneous velocity profile at the inlet.
Although the qualitative effect is similar, imposing nonhomogeneous inlet velocity profile improves simulation accuracy.
Scenario III
So far in order to compute pressure and velocity fields, Navier-Stokes equations were being solved. It’s worth pointing out that for all cases velocity magnitude didn’t exceed the value of 0.35m/s, so these are basically flows where inertial forces due to an advective motion of fluid are comparable with viscous forces. For sufficiently low-velocity magnitudes Stokes equations may be used as a mathematical model. To examine this approach a simulation with an inlet velocity of 0.025 m/s was performed. Results compared with previous cases are shown in Figure 13 where the Stokes pressure drop line is obtained by extrapolation of the single pressure drop data.
Figure 13. The discrepancy of pressure drops with Stokes and Navier-Stokes models.
Note that the difference between the linear and inertial term is gigantic. Hence this velocity range is clearly in the inertial regime. In order to have a closer look for smaller values of velocity, a series of simulations was carried out. You can see the results in Figure 14 below. Note that the discrepancy between the linear and non-linear approach diminishes from 50% at U=0.02m/s to 10% at U = 0.004 m/s what we wanted to demonstrate.
Figure 14. The discrepancy of pressure drops with Stokes and Navier-Stokes models for low Reynolds numbers.
Summary
This case shows that QuickerSim CFD Toolbox for MATLAB is capable of simulating of flows through materials with complex, porous structure. In this way, the whole process of importing files containing the porous geometry, generating the computational mesh and running CFD simulation can be performed in MATLAB software.
Still, mind that the flows in porous media can be computed using the averaged approaches – you can check out the samples here: Tutorial 3 and Tutorial 11.
References
S. Whitaker. Flow in porous media I: A theoretical derivation of Darcy’s law. Transport in Porous Media, 1(1):3–25, 1986. ISSN 01693913
S.Whitaker. The Forchheimer equation: A theoretical development. Transport in Porous Media, 25(1):27–61, 1996. ISSN 0169-3913.
W.Regulski, J. Szumbarski, Ł. Łaniewski-Wołłk, K. Gumowski, J. Skibiński, M. Wichrowski. Pressure drop in flow across ceramic foams—A numerical and experimental study, Chemical Engineering Science 137: 320-337, 2015
P.R. Parthasarathy, P. Habisreuther, and N. Zarzalis. A study of pressure drop in reticulated ceramic sponges using direct pore level simulation. Chemical Engineering Science, 147:91–99, 2016. ISSN 00092509.