Introduction

The article is devoted to the Darcy’s law research for an ideal gas in three-dimensional porous medium. The article is a continuation of [1], where the two-dimensional case was considered. It is assumed that it is possible to establish relation between the flow speed and the gas filtration rate. The flow speed and other macroscopic values are calculated using statistical estimates. More about this is described in [2, 3, 4].

The kinetic model of an ideal collisionless gas is a dynamic system consisting of a statistically large number of particles noninteracting with each other but interacting with boundaries. We have created software that allows calculating a huge number of particle trajectories by Monte Carlo method. In the modeling process such a problem, there is a huge amount of numerical data. Data visualization in required value graphs and gas flow images allows presenting in a visual form and qualitatively evaluating the results.

1. Ideal gas filtration model description

The mathematical model of ideal gas filtration in porous medium is described by the following equation system [5-9]:

(1)

where  is the permeability of a medium;  is the dynamic viscosity of the fluid;  is the particle mean speed;  is the mean free path.

Next we consider the one-dimensional case where the gas flow is in a steady state

(2)

After integration (1) the continuity equation is as follows:

(3)

By substituting the first and the last equations from the system (1) to (3) with subsequent integration the formula of ideal gas filtration in porous medium is defined by

(4)

a) configuration 1 ()

b) configuration 2 ()

c) configuration 3 ()

Figure 1. Computational scheme: the porous medium is red; the elementary domain boundaries are green; the particles are yellow

where è  are the elementary domain numbers,  is the distance between the elementary domains (the porous medium part length).

The elementary domain is a selected part of some simulated domain. Each subdomain may have a certain particle number at the time moment. The statistical estimate macroscopic values are calculated by using the local particle parameters belonging to the elementary domain.

2. Computational problem statement

Let the simulated domain is a rectangular parallelepiped unlimited in the positive direction of the x-axis. The domain is divided into 3 subdomain.

At the initial time, the first empty subdomain is uniformly filled with gas particles. The particle speed direction is also determined randomly with a uniform distribution. The ideal gas flow is directed from left to right. The second subdomain is porous medium consisting of solid spheres, the centers of which are located in the cubic grid nodes. The sphere radius is . Figure 1 shows a computational scheme with three types of the porous medium configurations (see table 1). The last subdomain is empty and unlimited.

Table 1. Porous medium configuration parameters

On the other hand, the simulated domain contains 6 elementary domains (see figure 1). The elementary domain boundaries are highlighted in green. Elementary domains are numbered from left to right along the x-Axis starting from zero.

The sphere number for the th configuration of the porous medium is determined by

(5)

3. Software description

The trajectory of particle, which can interact with the boundaries, is described by

(6)

(7)

where  is particle coordinates;  is particle velocity after the th interaction with the boundary;  is time step;  is the collision number of the particle with the boundaries;  is a function of the particle interaction with the boundary;  is a surface normal of the th boundary at the collision point at time  of the th step.

If the boundary is a plane, then the time moment  of crossing the particle trajectory and the boundary is defined by

(8)

where ; .

In other case when the boundary is a sphere, the crossing time moment  is as follows [10]:

(9)

where ; ; .

Let the intersection point of the particle trajectory and the boundary was found by using (8), (9). As a computational result this point can be on the other side of the half-space bounded by the surface, due to the computational error. This can lead to reflection and incorrect particle trajectory (see figure 2). The solution is to add a correction coefficient  that moves the particle back along the trajectory.

Figure 2. Illustration explaining the use of the correction coefficient : the thick line is the boundary; dotted line is the particle trajectory; the line is the particle trajectory including the correction coefficient

In the case when the particle is reflected from the boundary according to the specular reflection [11], the particle velocity after the collision with the moving boundary is

(10)

Figure 3 shows the scheme calculating the trajectories of the particles, which can interact with the boundaries.

Figure 3. The scheme calculating the trajectories of the particles, which can interact with the boundaries

This scheme allows accelerating calculations by using vectorization operation [12]. Although the data sorting cannot be vectorized, the data swapping can reduce the operation number by moving the particles noninteracting with the boundaries at the group end during time step .

It is known that division and root extraction operations are slower than addition (subtraction) and multiplication operations. If it's possible, the "slow" operations before additional conditions in cycles are transferred inside the condition, changing the condition of addition and multiplication additional operations. It can lead to an increase in the calculation speed. The solution presented in [13] is applied to quickly extract the inverse square root, which is used in the normalization of vectors.

Let the porous medium consists of spheres located in the cubic grid nodes. Figure 4 shows a cell belonging to such a grid. If the particle can move less than half of the cell length per time step, it can be argued that the particle can interact only with one of 8 spheres, the centers of which are located at the cell nodes. This also leads to a significant increase of the particle trajectory calculation speed.

Figure 4. The cubic grid cell of porous medium

Figure 5 shows the general scheme of software [14].

Figure 5. The general scheme of software

The model initial state is defined as the parameter initialization of particles by using a pseudo-random number generator, boundaries and elementary domains. For each elementary domain is accumulated the sums of the particle parameters (the sums of particle mass, the velocity components and the velocity component squares) belonging to the elementary domain at the time moment . These sums are applied to calculate macroscopic parameters statistical estimates by using the formulas presented in [5]. Intermediate data saving includes phase space and particle parameter sums saving, which is applied at the last stage for data visualization and calculation of macroscopic parameter statistical estimates. The last stage (data processing) may also include plotting graphs and the absolute and relative error calculation.

Another idea of computational acceleration is splitting the particle array into large blocks. Each program accumulating particle parameter sums works with a block. Macroscopic parameter statistical estimates are calculated by other program after computing the particle parameter sums for all blocks.

It is known that Linux OS is applied to manage of modern supercomputing systems. The developed software is oriented for such OS. Shell scripts are applied to coordinate all programs, each of which has a special purpose. Algorithms that calculate the trajectories of the particles interacting with boundaries and compute macroscopic parameter statistical estimates were implemented in the C programming language, allowing relatively efficiently (compared with most other languages, the exceptions are Fortran and C++) to use the computational resources of modern central processors. The C language allows to abstract from the processor architecture by using the features such as preprocessor directives and special options at the compilation stage. The Mersenne Twister was applied as a pseudorandom numbers generator, which exists in the MKL library of Intel. Gas flow is visualized of using the OpenGL API. Visualization of the field of macroscopic parameter statistical estimates is made using the TecPlot application package, which requires files with input data in a special format. Scripts written in the Python language are applied to plot graphs using Matplotlib 2D plotting library. This language supports dynamic typing and contains a huge number of extra libraries (for example, Matplotlib).

4. Modelling and result analysis

Computational experiment series with different particle number for three porous media configurations were conducted. Figure 6 presents a gas flow visualization in porous medium where .

à) ;

b) ;

c) ;

Figure 6. Gas flow for the different porous medium configurations at various time moments

As a part of the computational experiment the inverse problem of finding a coefficient  linking the flow speed (the hydrodynamic velocity statistical estimation) and the filtration rate (4) depending on the configuration of the porous medium was solved.

Since the coefficient  contains some values that can be calculated (11), it is sufficient to find the following relation for each configuration .

(11)

The coefficient  calculation is divided into the following stages:

1. For each pair of velocities  and  with the same numbers  and  the coefficient value  is given by the equation:

(12)

The equation (12) is based on the least square method with a time interval , when the gas flow has a quasi-steady state, where time step is . The values  are presented in table 2, depending on  and  for different particle number.

Table 2

The coefficient values  depending on the particle number , porous medium configuration and elementary domain numbers ,

2. The obtained values  are averaged for each porous medium configuration. Figures 7-10 show some results, where the coefficient values  for each porous medium configuration are taken from table 2.

a) ;                                              b) ;

c) ;                                               d) ;

e) ;                                              f) ;

Figure 7. Graphs of flow speed  and filtration rate  depending on time  for different porous medium configurations with different particle number

a) ;                                               b) ;

c) ;                                               d) ;

e) ;                                               f) ;

Figure 8. Graphs of flow speed  and filtration rate , pressure and temperature statistical estimation depending on time  for

a) ;                                              b) ;

c) ;                                              d) ;

e) ;                                              f) ;

Figure 9. Graphs of flow speed  and filtration rate , pressure and temperature statistical estimation depending on time  for

a) ;                                              b) ;

c) ;                                              d) ;

e) ;                                              f) ;

Figure 10. Graphs of flow speed  and filtration rate , pressure and temperature statistical estimation depending on time  for

The macroscopic parameter statistical estimate field is constructed by dividing the simulated domain into subdomains. Each subdomain is an elementary domain. Figure 11 illustrates the temperature and pressure distribution for the three porous medium configurations. The grid size (the elementary domain number in each dimension) is equal to 128õ32õ32.

The software performance evaluation is executed on Intel Xeon E5-2690 V2 processor with a peak performance of 240 GFlops DP for 10 cores.

Table 3 presents the software performance per core for calculating the particle trajectory and including data saving where the particle number is , the time step number is 300, the elementary domain number is 6.

Table 3

The performance evaluation on Intel Xeon E5-2690 V2 processor per a core depending on the porous medium configuration and the time step

Conclusion

We created software for numerical research of an ideal collisionless gas kinetic problem in the three-dimensional space. The software allows to visualize gas flow.

It is possible to establish the following facts by using data visualization, which substantially facilitated the research:

1. The filtration rate and gas flow speed graph comparison allowed to confirm qualitatively that there is a dependence on these values for porous medium configurations.

2. The gas flow visualization helped to explain the initial divergence of the flow speed and the gas filtration rate. The reason was the lack of a sufficient gas particle number in porous medium, i.e., the gas flow has a quasi-steady state.

The relation between the flow speed (a statistical assessment of the hydrodynamic particle velocity) and the filtration rate of an ideal collisionless gas has been identified. Darcy’s law depends on the porous medium structure for an ideal gas without viscosity.

This work was supported by RFBR Grants 16-29-15105, 18-01-00343.

References

1. Galkin V. A., Bykovskikh D. A., Gavrilenko T. V., Stulov P. A. Filtratsionnaya model dvizheniya idealnogo gaza v poristoy srede // Vestnik kibernetiki. Elektr. Zhurn. 2016. 4 (24). pp. 50–57.

2. Bykovskikh D. A., Galkin V. A. O vychislitelnom teste dlya odnoy modeli besstolknovitelnogo idealnogo gaza // Vestnik kibernetiki. Elektr. Zhurn. 2017. 3 (27). pp. 100-120.

3. Cercignani C. Mathematical methods in kinetic theory. Milano: Politecnico di Milano, 1969.

4. Galkin V. A. Analiz matematicheskikh modeley: sistemy zakonov sokhraneniya, uravneniya Boltsmana i Smolukhovskogo . M.: BINOM. Laboratoriya znaniy, 2011.

5. Leontev N. E. Osnovy teorii filtratsii. M.: Izd-vo TSPI pri mekh.-mat. fak-tu MGU, 2009.

6. Masket M. Techenie odnorodnykh zhidkostey v poristoy srede. M.: NITS «Regulyarnaya i khaoticheskaya dinamika», 2004.

7. Wu Y.-S., Pruess K., Persoff P. Steady and Transient Analytical Solutions for Gas Flow in Porous Media with Klinkenberg Effects. California: Lawrence Berkeley National Laboratory, 1996.

8. Trapeznikova M. A., Belotserkovska M. S., Chetverushkin B. N. Analog of kinetically consistent schemes for simulation of a filtration problem // Matem. Mod. 2002. Vol. 14, No 10, pp.  69-76.

9. Schmidt B. E. Compressible flow through porous media with application to injection: internal report for Caltech hypersonics group FM. Pasadena: California Institute of Technology, 2014.

10. Tsissarzh V. V., Marusik R. I. Matematicheskie metody kompyuternoy grafiki. K.: FAKT, 2004.

11. Galperin G. A., Chernov N. I. Billiardy i khaos. M.: Znanie, 1991.

12. Bik A. J. The Software Vectorization Handbook: Applying Intel Multimedia Extensions for Maximum Performance. Intel Press, 2004.

13. Robertson M. A Brief History of InvSqrt. Department of Computer Science & Applied Statistics, 2012.

14. Bird G. A. Molecular gas dynamics and the direct simulation of gas flows. N.Y.: Oxford University Press, 1994.