Numerical Study of the Flow Through Porous Structures Built from Gray–Scott Patterns

Numerical solutions of the reaction–diffusion system of equations in three dimensions provide a systematic procedure to generate a variety of porous media structures of a chosen porosity. A way to characterize such geometries is by their Minkowski functional densities and their average pore diameter. We generated multiple porous geometries of the granular and foam-like types. Hydraulic flow through these series of geometries is simulated to study the dependence of the permeability k and the Forchheimer parameter β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} on the chosen geometrical identifiers. Both k and β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} exhibit two distinct behaviors aligned with the two distinct generated types of geometries when scaled with the Euler characteristic. A good correlation for dimensionless k as a function of the Minkowski functional densities and the average pore diameter was found for both types of geometries with an accuracy of 15% and 18%, respectively, while no good correlation was found for β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}, implying that more geometrical information is needed beyond the proposed one. From the correlations, it was found that, together with the porosity, the mean curvature is also a good parameter to characterize the permeability. Porous media generation based on solutions to reaction-diffusion equations Model porous structures characterized through Minkowski functionals and average pore diameter Scaling based on the Euler characteristic reveals two distinct flow behaviors. Porous media generation based on solutions to reaction-diffusion equations Model porous structures characterized through Minkowski functionals and average pore diameter Scaling based on the Euler characteristic reveals two distinct flow behaviors.


Introduction
In the phenomena of flow-through porous media, the geometrical characterization of the pore space on a representative elementary volume (REV) relevant to determine flow behavior remains a fundamental open question of great practical implications.Transport phenomena in porous media are involved in a wide range of applications, from biophysics (Khaled and Vafai 2003) to petroleum engineering (Parker 1989), and include porous structures of extreme diversity.The search for a set of geometrical parameters that fully characterize a porous structure in order to obtain constitutive relations for the transport properties for a wide range of applications seems to be a formidable task.
Originally porous structure was thought to be sufficiently parameterized by its porosity , pore size d and tortuosity , and most commonly used approximate relations for physical properties such as permeability, formation factor, electric tortuosity or capillary pressure involve these parameters (Bear 1988).However, dependence on , d, and alone has proved to be insufficient given the diversity of porous structures, and the predictive power of such relations is very limited.During this century, the inclusion of Minkowski functionals (porosity is one of them) in porous geometry characterization has been explored with varying degrees of success (Armstrong et al. 2019;Mecke 2008;Slotte et al. 2020;Vogel et al. 2010).
To determine flow properties on porous media, specially for the permeability dependence on geometrical characteristics of porous structures, researchers have employed computational fluid dynamics simulations for the past decades.For the simplicity involved in including complex geometries and parallelization, the lattice-Boltzmann method is a common choice for flow simulation in porous media (see most of the cited references in the rest of this paragraph).To perform flow simulations, it is necessary to determine the geometries of the flow domain, the pore space, and to this end one can either reconstruct real porous material or artificially create such intricate domains.Information from real media samples is commonly obtained through X-ray tomography, from which a three-dimensional direct measurement of the pores space can be obtained and used to study the flow properties (Vogel et al. 2010;Mahmoodlu et al. 2016;Santos et al. 2022).Artificially created porous geometries for flow studies include the use of networks of tubes (Blunt and King 1990), fractal geometries (Lenormand 1990), and volume tessellation techniques (Xiao and Yin 2016).
In this work, we take advantage of the known feature of the solutions of some reaction--diffusion equations to mimic patterns observed in nature, to generate geometries that resemble those of porous media (Halatek and Frey 2018).We choose to work with the Gray-Scott system of equations to produce the basis from which families of model porous geometries were generated.Geometrical parameters, including the Minkowski functional and the average pore diameter, were computed for the set of model porous structures.Lattice-Boltzmann simulations were performed to compute the flow through the created geometries and to study the relation of the volumetric flux to the pressure gradient, and determine any possible relation with the geometric parameters.
The manuscript is organized as follows: first in Sect.2, we do a brief introduction of the porous geometry parameters that are commonly used to characterize it.In Sect.3, a detailed description of the porous geometry generation is presented.Section 4 presents the flow simulation details.Results for permeability as a function of the Minkowski functionals and the average pore diameter are found in Sect.5, and similar results for the Forchheimer coefficient can be found in Sect.6.Finally, in Sect.7 a summary and conclusion is presented.

Porous Geometry Characterization
Let us start with a brief review of Minkowski functionals.A porous geometry can be thought as a d-dimensional compact object X embedded in a d-dimensional Euclidean space Ω .Hadwiger's characterization theorem (Klain 1995) demonstrates that any additive, motion invariant, and conditionally continuous functional on a finite union of compact convex sets in Euclidean space can be written as a linear combination of d + 1 geometric measures known as intrinsic volumes or Minkowski functionals (MF) (Mecke 2008).In d = 3 , the four Minkowski functionals are defined as: with dv the volume element on Ω , ds the surface element on X with principal radii of curvature given by r 1 and r 2 , X the boundary of X, ( X) the Euler characteristic of the bounding surface, and (X) the Euler characteristic of the solid object (Ohser and Mück- lich 2000).The MF comprises a set of measurements in various dimensions: M 0 corre- sponds to the volume, M 1 to the surface area, M 2 to the mean curvature, and is a topo- logical quantity (dimensionless) related to the number of isolated objects (N), redundant loops (O), and cavities (C) through (X) = N − O + C .Permeabilities are not an additive quantity, and consequently, they are not covered by Hadwiger's theorem.However, if they happen to be functions of a set of additive measures, then they will also be a function of the MF.Assuming the main dependence of the permeabilities come from their MF dependence, we will make use of them as a way to parametrize any given porous geometry.We expect the functional dependence of the permeability to be independent of the scale of the material sample so we will work with the volume-normalized MF defined as with V the volume of Ω .Note that m 0 coincides with the standard definition of porosity .The computation of the MF densities of the computer-generated geometries to be used will be performed through the algorithm outlined in (Michielsen and Raedt 2001).The MF densities are properties associated with the integral geometry of the sample so the information regarding the detailed pore structure is not contained by them.Additionally, one of the basic elements of pore structure information is the average pore diameter size d defined as: where f(r) is the normalized pore diameter frequency distribution, i.e., ∫ ∞ 0 f (r)dr = 1 , such distribution can be discretely estimated through the morphological operations of erosion and dilation as outlined in (Vogel et al. 2010).It will be useful for later convenience to write dimensionless versions of the MF densities and the average pore diameter.To that end, we will use the following topological length H as characteristic length for this purpose1 Given the topological nature of H, its computation carries no associated numerical error and H itself is insensitive to the roughness of the solid-fluid inter-phase.Using H, we can define the independent dimensionless variables These variables will be relevant when analyzing the dependence of the permeability with the MF densities.Note that despise using H we will still use d as the relevant characteristic length when computing the Reynolds number.

Porous Geometry Generation
The systematic generation of porous geometries is based on a special treatment, outlined in the rest of this section, of a class of solutions to the Gray-Scott reaction-diffusion system of equations given by where D u and D v are diffusion coefficients associated with the species concentration u and v, respectively, f is known as the feed parameter, k is known as the kill parameter.Periodic boundary conditions were used in all boundaries, and spherical initial distributions were given by: where H e (x) is the Heaviside step function, U [0,1] is a random variable sourced from a uni- form distribution defined between [0, 1], and (r, , ) are spherical coordinates centered at the cube center, and r c is a constant.Table 1 shows the values of the coefficients used in ( 9)-( 10) that produced a series, or families, of distinct structures that we named from A to H.
The system (9)-( 10) is solved numerically via a finite difference scheme on a cube Ω L GS of size L GS with periodic boundary conditions.As time evolves, the patterns generated by the solutions reveal a class of solutions whose iso-surfaces are reminiscent of the internal boundaries encountered on porous media, see Fig. 1.
The initial configuration evolves into complex patterns that extend all over the computational domain, like those shown in figure 1.Once a steady state has been reached the   1.The solutions were obtained on a cube of length L GS = 192 computed distributions u, or v, can be transformed into a porous geometry with porosity through the following procedure.2 1.The support of the numerical solution {u, v} is reduced to a discrete set of points located on a cube ( Ω L � ) of size L ′ < L GS centered on the original L GS sized cube Ω L GS .Therefore, we will consider only a subset of the found solution to generate the porous geometry.This removes some undesired boundary behavior that causes the geometry to not have the properties of a representative elementary volume (REV).3 2. We find an iso-surface u c , or v c equivalently and treat it as a basic fluid-solid boundary of the porous geometry.This means that the domain where u l < u c , or u l > u c alternatively, is taken as the pore space, while its complement its taken as the solid space, see Fig. 2. The choice u l < u c or u l > u c , to decide which domain is the pore space, will produce two distinct families of geometries, that are in a sense complementary and that show very different flow properties.This is the reason why, in Table 1, there are two families per choice of parameters in the Gray-Scott system of equations ( 9)-( 10), and so family A is complementary to family B in the sense stated above.This results in a preliminary porous geometry defined on a cubic lattice of size L ′ .The information for this prelimi- nary geometry is saved on a function P(x) with x ∈ Ω L � and P(x) ∈ {0, 1} where 0 and 1 denote the belonging to the solid and pore space, respectively.3. The preliminary porous geometry might have pore diameters as small as a single computational node causing potential problems on those sections when a flow is run through Fig. 2 Transverse cut of the two solutions shown in figure 1 was only the L = 96 central cube has been considered as support for the solution.On the left part of each image, a particular iso-surface u c is shown, and on the right-hand side, we see in black the union of iso-surfaces u l < u c and in gray we have the union of iso-surfaces u l > u c .In both cases, if the black section is chosen as the pore space a porosity of = 0.55 is found otherwise a porosity of = 0.45 is found.Note that the two geometries that can be generated from a single solution have different structure, see figure 6 for more examples.It might seem that if the black section is taken as the pore space it will not be connected; however, we have to remember that the geometry is three dimensional and the connection will occurred through other planes it. 4To remove this potential issue, we increase the resolution of the sample by expanding the sample onto a cube Ω L s of length L s = sL � with s > 1 .This is done by turning each point x on the original cube Ω L � into a cube of length s.This procedure increases the minimum pore radius to s 3 computational nodes.We found {s = 3, L � = 96} to be a good compromise between resolution and computational speed.4. The potential porous geometry is further enhanced by removing any sections of the geometry that are disconnected from the main geometry body.The first step involves choosing one pair of parallel faces of Ω L s and identifying one as the inlet and the other as the outlet.We then take all pore points x belonging to the inlet and check which pore points y are connected to them.Any pore point z not belonging to this group of connected points is converted into a solid point.This procedures ensures that the whole pore space is connected to the inlet.5. We use the algorithms of (Michielsen and Raedt 2001) to compute the MF densities of the potential porous geometry.To verify that the sample corresponds to a representa-  4 We noted that a couple of these pathological points existed in all generated geometries.tive elementary volume (REV), we repeat this computation for all possible subsamples of size 0 < l < L s .If the variations of the MF densities and the average pore diameter around the size of the sample are small, we can consider our sample to be larger than a single REV, see figure 3. 6.Once we have verified our sample behaves and is larger than a single REV, we can use the algorithm presented in (Vogel et al. 2010) to compute the average pore radius of the geometry, see Fig. 4, to visualize the pore distribution f(x) as well as the cumulative pore size d(x) = ∫ x 0 rf (r)dr .The porous geometry is now ready to be subject to hydraulic flow.
Following the previously outlined procedure, we generated 154 different geometries classified into 8 different geometry families.Each family corresponds to a different solution to the reaction-diffusion system (9)-(10).All solutions were obtained on a cube of side L through a finite difference scheme taking Δx = 0.01 and Δt = 0.9Δx 3 4D u .In the table below the parameters and initial conditions in lattice units are shown for each of the families,5 see Fig. 6 to visualize the structure of some of these families.The range of the dimensionless densities 1−3,d and the number of data sets analyzed for each family set are specified in Table 2.
In summary, the pore space of each family is generated from level surfaces of a given Gray-Scott solution.The level surfaces are chosen such that the porosity/volume associated with their union gives out a specific porosity 0 .Due to this, we expect that, for a given family, the rest of the Minkowski functionals and the pore diameter to be functions of the porosity This dependence can be seen in Fig. 5. From this figure, we can also see that there is a clear division between geometries in two groups or types.We will group together geometries A, C, E, and G into one group and the remaining B, D, F, and H into another one.Each type of geometries corresponds to one of the two ways of constructing a porous structure from the solution to the Gray-Scott system, as described in step 2 of the outlined procedure at the beginning of this section.Some clear differences between the two types of geometries are: 1.The sign of d 1 ∕d 0 is different for each of the types of geometries, see Fig. 5a. 2. For {A, C, E, G} families 2 takes positive and negative values, while for {B, D, F, H} families it only takes positive values.From figure 5b), we can notice that both types do not overlap on the ( 2 , 0 ) space.We should note that they also do not overlap on the ( 2 , d ) space, see Fig. 5d. 3.For {A, C, E, G} families d d ∕d 0 changes sign in the range of porosities that were studied, while for {B, D, F, H} families geometries d d ∕d 0 remain positive.See Fig. 5c. 4. From a visual point of view, both geometries are also quite distinct.{A, C, E, G} families geometries resemble those generated from lumping together grains, see Figs. 6a andc.{B, D, F, H} families resemble those of foams, see Figs. 6b and d.
Although properties 1 to 3 are probably exclusive to geometries generated through the method outlined in this section, it is possible to classify other porous geometries into one of the two types by using the visual cues discussed in point 4 of the previous list.For example, a porous geometry made out of packed grains will share similarities with {A, C, E, G} families, while a foam-like porous geometry will share similarities with {B, D, F, H} families.We will then refer to the first group of geometries as granular geometries, and to the second group of geometries as foam-like geometries.

Hydrodynamic Setup
The geometries generated using the previously outlined algorithm define a cubic lattice of length L = 288 computational nodes, where each node can be part of either the fluid or the solid domain.Note that by construction the solid domain is a connected space.We now want to study the flow through these geometries in the case where one of the faces of the cube acts as an inlet at fixed pressure p i , while its opposite face acts as an outlet at fixed  6 To simulate this setup, we used a D3Q19 lattice-Boltzmann method (LBM) where a distribution function f k is defined and computed over the lattice identified as the porous geometry.The distribution function is then used to compute the fluid velocity u at the lattice nodes.Lattice spacing as well as time steps can be conveniently set to unity.At every node r in the lattice, the distribution functions evolve in time according to The coefficient represents a relaxation time and is related to the fluid kinematic viscosity = − 1 2 3 .We choose = 0.55 , a common value for the method to be stable at moderate Reynolds number (Zhao 2013).The local equilibrium distribution function f eq k is given by: The equilibrium distributions depend on the macroscopic fields u , and , the mass density and must be computed every time step through The constants k in the equilibrium function definition (15) take the values 0 = 1∕3 , k = 1∕18 for k=1,...,6 and k = 1∕36 for k = 7, ..., 18 .The set of microscopic velocities {e k ∶ k = 0, ..., 18} is given by: Equation ( 14), with the chosen microscopic velocities, provides an algorithm for updating all the distribution functions f k at a given node in the lattice, as long as its 18 nearest neighbors in the lattice are inside the fluid domain.For nodes close enough to a solid wall, the distribution functions coming from neighboring nodes outside the fluid domain must be provided as a boundary condition for the method.When the neighboring nodes correspond to the solid domain of the porous geometry, including the four closed off faces of the cubic domain, the half-way bounce-back conditions will be used (Zhang et al. 2012).When the neighboring nodes go outside the computational domain7 at the inlet and outlet, we will use the boundary conditions proposed in (Zou and He 1996;Hecht and Harting 2008) that fix the pressure to a constant value at both inlet and outlet.On the aforementioned setup, and for a slow enough flow, Darcy law (Darcy 1856;Bear 1988;Whitaker 1986) is expected to hold, where Δp is the pressure gradient between inlet and outlet, L is the length of the sample cube, is the dynamic viscosity of the fluid, Q is the discharge rate at the outlet, and K is the permeability.As discussed in the introduction, the permeability is assumed to depend exclusively on the geometric properties of the sample.To be more precise regarding the validity of ( 19), it is useful to define the pore level Reynolds number Re as where the average pore diameter d , see ( 6), is used as the pore characteristic length, Q L 2 corresponds to the seepage velocity and is used as the characteristic speed at the pore level.
We then expect Darcy to hold for Re  < 1 , while for Re  >> 1 a quadratic term known as the Forchheimer term is added and ( 19) is expected to take the quadratic form where we will refer to as the Forchheimer parameter.It will be convenient to rewrite the general form (21) as a function of the pore level Reynolds number as with the dimensionless permeability k = K H 2 , dimensionless Forchheimer coefficient b = H , and Δ * given by To compute the permeability k and Forchheimer coefficient b, we just need to analyze the relation between pressure gradient and discharge rate once the flow has reached a steady state.Some examples of steady stream lines, in different geometry families, obtained from the LBM method are shown in Fig. 6.We can notice that the shape of the stream lines differs for granular-and foam-like geometries, see Fig. 6.For granular geometries, the stream lines are what you expect from a fluid avoiding obstacles on its path; however, for foamlike geometries we can see that the stream lines take advantage of the three-dimensional space and become more tortuous than their granular counterparts.We expect tortuosity to be an alternative classifying factor for the two types of observed geometries; however, we did not explore this possibility on this work. Δp

Permeability
The permeability associated with any given data set can be obtained by computing the discharge rate at the outlet for multiple stationary numerical solutions on a range of small Reynolds number Re  ≲ 1 .Fitting the data to Darcy law ( 19) allows for k = K H 2 to be computed, see Fig. 7 for some examples.
We will now assume that8 k = k 0 , 1 , 2 , d , and our aim will be to find a good enough fit for the functional form of k 0 , 1 , 2 , d and try to get a picture of how all Minkowski functionals affect the permeability.As mentioned in Sect.3, due to the way the geometries were generated the Minkowski functionals are not independent.This implies that for any given family the dimensionless permeability k can be written as a single parameter function,9 for example we can choose the porosity as the independent parameter and write k = k( 0 ) .We exemplify this in Fig. 8, where we show k vs 0 and k vs. 2 , together with a fit for both cases for each of the families.In Fig. 8, all the generated structures are included, grouped in geometry families as defined in Sect.3. In plotting k vs. 0 and k vs The fit parameters can be found in Table 3.The fit k( 0 ) was inspired by the Kozeny-Car- man equation.The fit satisfies the expected11 relation k( 0 → 0) = 0 .For the fit k( 2 ) we ( 24) found that, if plotted in semilog scale, a polynomial behavior on 2 is needed to account for the multiple inflection points on the curve k( 2 ) , while the exponentiation is necessary to guarantee that 12 k( 2 → ∞) = 0 .Note that a power law in 2 might not work in general as 2 can take both positive and negative values.
For foam-like geometries, we will instead use the following fits: The fit parameters can be found in Table 4, note that due to the small number of data points we will use a three-parameter fit in contrast to the four parameter fit used in ( 24)-( 25).
Given the rapid increase in k with 0 that occurs on foam-like geometries, we used a power law combined with an exponential term.For k( 2 ) , we kept the exponential form but removed the cubic term, the reasoning for this is that there are not enough data points to correctly fit such parameter and there are no inflection points in the range of explored 2 .
From the Δ( i ) columns of Tables 3 and 4, we can note that for most families k( 2) is a slightly better fit than k( 0 ) .This seems to indicate that the curvature 2 is probably a natural parameter to describe the permeability in conjunction with the porosity.We can also notice that the fit parameters for all granular geometries are roughly of the same order.The same holds true for foam-like geometries.We will then assume that we can incorporate all family geometries into a single fit k( 0 , 1 , 2 , d ) .Based on the fits ( 24)-( 27) as well as the fitting proposed by Slotte et al. (2020) for two-dimensional porous we were able to fit the numerical data for granular geometries into the trend function 13   with the best fit values ( 26) In Fig. 9, a plot of the relative error for the dimensionless permeability can be found, and from it, we can see that the trend line approximates the expected value for each of the types of geometries with an error of less than 15% for granular geometries and less than 18% for foam-like geometries.
We can finally note that for both fits ( 28) and ( 30) the limit k(0, 1 , ∞, d ) = 0 is satisfied, and there is an exponential dependence on the dimensionless curvature 2 and a power law on the porosity 0 .

Forchheimer Parameter
The Forchheimer parameter associated with any given data set can be obtained by computing the discharge rate at the outlet for multiple stationary numerical solutions on a larger range of Reynolds number Re  > 1 .Fitting the data to the Forchheimer correction to Darcy law (21) allows for b = H to be computed, see Fig. 7 for some examples.Ideally we would like to find if there is a good enough fit for the functional form of b( 0 , 1 , 2 , d ) for both types of geometries and identify the main dependence on the different Minkowski functionals.Just as for the permeability, it should be possible to first do single parameter   24) and k( 2 ) using ( 25) together with the maximum relative error Δ( i ) ≡ 100max| | with k the permeability as given by the numerical simulations  26) and k( 2) using ( 27) together with the maximum relative error with k the permeability as given by the numerical simulations Family a fits; however, we were only able to find simple fits for b( 2 ) .In Fig. 10b, b( 2 ) is shown for all geometries together with its fit.For completeness, we show in Fig. 10a the behavior of b( 0) , from it we can notice that b( 0 ) rapidly increases as 0 decreases and, as expected, the behavior of each family seems to resemble the behavior of other families belonging to the same type.For granular and foam-like geometries, we used a fit of the forms with k 3 set to zero for foam-like geometries given the lower quantity of numerical data in comparison with granular geometries.The fit parameters for both types of geometries, together with the maximum relative error Δ( 2 ) , are presented in Table 5.
We can already noticed that the fit is not as good as the one for the permeability for all types of families.However, it is quite notable that the exponential form (32) can capture most of the exhibited behavior of b( 2 ) .Unlike k, we were unable to find a good fit func- tion for each of the types of geometries solely using the Minkowski functionals and the average pore diameter as in (28) or (30).Given that we expect each family to share geometric features, we are inclined to conclude that additional geometric parameters that specify the more detailed structure of the geometry become relevant.

Conclusion
Making use of the reaction-diffusion equation, we were able to successfully implement a procedure that systematically generates porous geometries with distinct features for any given porosity.Such geometries were then characterized by their Minkowski functional densities and their average pore diameter.We generated 154 distinct geometries, arranged in 8 families sharing similar visual structure features and divided into two types of geometries with distinct geometrical features.The first type resembles granular porous geometries, while the second one resembles foam-like geometries.We subjected them numerically to a pressure gradient-driven hydraulic flow.Based on the mean pore diameter, the Reynolds number at pore level explored ranged from 0.2 up to 17.
From the observed relation between the discharge rate and the pressure gradient, we were able to compute the permeability and Forchheimer parameter for all generated geometries.To analyze a possible relation between these coefficients, and the Minkowski functional densities and average pore diameter, we used a dimensionalization involving the length derived from the Euler characteristic density.Under such dimensionalization and for each of the families, we were able to write down single parameter fits for the permeability as function of the porosity 0 and the dimensionless mean curvature 2 , and for the Forchheimer parameter in terms of the dimensionless mean curvature 2 .Additionally the dimensionless permeability arranged itself with the two types of families allowing us to find a fit as function of the dimensionless Minkowski functionals and the dimensionless average pore diameter for each of both types of geometries with an expected error of less than 15% for the granular geometries and less than 18% for foam-like geometries.
The two types of geometries are (generation and visual wise) complementary to each other and occupy different regions of the ( 2 , 0 ) space, as seen in Fig. 5b), and different regions of the ( d , 2 ) space, as seen in Fig. 5d).It is unclear whether this is the only geo- metric distinction between them, although looking at the stream lines, see Fig. 6, it seems to be an indication that tortuosity might be a better discriminatory quantity between the two observe types of geometries, and we leave the exploration of this possibility to future work.
The Forchheimer parameter also grouped itself into two types of families; however, we were unable to find a good fit for each type of geometry in terms of the Minkowski functionals and the average pore diameter.This seems to indicate the Minkowski functionals and the average pore density is not enough to determine the Forchheimer parameter, and further geometry descriptors are probably needed.However, for each of the families a single-parameter fit, analogous to the one used for the permeability, in terms of the dimensionless mean curvature 2 was possible.The nature of 2 and its relevance on all found fits seem to indicate that it is a key Minkowski parameter when predicting either the dimensionless permeability k or the dimensionless Forchheimer parameter b, recalling that a dimensionalization using the Euler characteristic was used.It will be of interest to study the relevance of 2 on generic porous geometries.
Our results show that reaction-diffusion systems have the potential to produce model porous structures of arbitrary geometrical characteristics.Such structures are suitable for the study of the flow properties dependence on the geometric characterization of porous media.Numerical simulations suggest that scaling lengths with the Euler characteristic reveal a different behavior of the dimensionless permeability, and Forchheimer parameter, of the two types of porous structures studied.Our work can be extended to include inhomogeneities and anisotropies and can be used to numerically study flow in a wide range of porous structures.

Fig. 1
Fig. 1 Iso-surfaces for solutions to the reaction-diffusion system for two set of parameters {D u , D v , f , k} , the parameters on the figure on the left(right) correspond to family A/B(C/D) in Table1.The solutions were obtained on a cube of length L GS = 192

Fig. 3 Fig. 4 a
Fig. 3 Porosity as a function of the sample length size for four of the used data sets.The full sample at L s=3 = 288 is located beyond the point where the sample can be considered a REV.The same fluctuation behavior is observed for the rest of the data sets and for the remaining MF densities and average pore diameter

Fig. 5 a
Fig. 5 a Dimensionless-specific surface area 1 vs porosity 0 for all geometry families.b Dimensionless mean curvature density 2 vs porosity 0 for all geometry families.c Dimensionless average pore diameter d vs porosity 0 for all geometry families.d Dimensionless average pore diameter d vs mean curvature density 2 for all geometry families

Fig. 7 a
Fig. 7 a Data fit to Darcy law for eight data sets representing each of the distinct geometries.b Data fit to Darcy law with the Forchheimer correction for eight data sets representing each of the distinct geometries

Fig. 8 a
Fig.8a Dimensionless permeability k as a function of porosity 0 together with a fit of the form (24) for granular geometries and (26) for foam-like geometries.b Dimensionless permeability k as a function of dimensionless curvature density 2 together with a fit of the form (25) for granular geometries and (27) for foam-like geometries.All dimensional rescaling was made by using the length H = (−m 3 ) − 1 d = 1.58, k 0 = 3.52, k 1 = 5.12, k 2 =3.15, and k 3 = 0.88.k 1 =2.34, and k 2 = 1.85.

Fig. 9
Fig.9Relative error of the trend functions (28) and (30) for all used geometries

Table 2
Range of dimensionless MF densities (using the length H) and dimensionless mean pore diameter for all geometry families

Table 3
Granular geometries fit parameters for k( 0 ) using (

Table 5
Granular and foam-like geometries fit parameters for b( 2 ) using (32) together with the maximum relative errorΔ( i ) ≡ 100max|