The article is devoted to the Darcy’s law
research for an ideal gas in threedimensional porous medium. The article is a
continuation of [1], where the twodimensional 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.
The mathematical model of ideal gas
filtration in porous medium is described by the following equation system
[59]:
(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 onedimensional 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.
Let the simulated domain is a rectangular
parallelepiped unlimited in the positive direction of the xaxis. 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

Sphere radius,

Sphere number,

1

0.5

160

2

0.25

720

3

0.125

4000

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 xAxis starting from zero.
The sphere number for the th configuration of the
porous medium is determined by
(5)
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 halfspace 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 pseudorandom 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).
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 quasisteady 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 ,






1

2

3

4

1

2

3

4

1

2

3

4

10^{5}

2

0,0616




0,0596




0,0462




3

0,1498

0,0462



0,0648

0,0578



0,0598

0,0506



4

0,2132

0,1540

0,0796


0,0686

0,0624

0,0575


0,0646

0,0592

0,0527


5

0,2514

0,2288

0,1741

0,0902

0,0315

0,0263

0,0203

0,0126

0,0665

0,0618

0,0579

0,0528


0,1449

0,0461

0,0572

10^{6}

2

0,1476




0,1079




0,0596




3

0,2290

0,1912



0,1331

0,1169



0,0648

0,0578



4

0,2648

0,2431

0,1944


0,1402

0,1264

0,1100


0,0686

0,0624

0,0575


5

0,2886

0,2692

0,2437

0,2117

0,1443

0,1321

0,1197

0,1110

0,0692

0,0638

0,0598

0,0569


0,2283

0,1242

0,0620

10^{7}

2

0,1946




0,1149




0,0604




3

0,2417

0,2288



0,1299

0,1181



0,0664

0,0603



4

0,2684

0,2541

0,2291


0,1390

0,1275

0,1169


0,0691

0,0632

0,0580


5

0,2847

0,2683

0,2463

0,2254

0,1446

0,1337

0,1240

0,1161

0,0699

0,0645

0,0600

0,0569


0,2441

0,1265

0,0629

10^{8}

2

0,2018




0,1152




0,0607




3

0,2450

0,2359



0,1308

0,1201



0,0663

0,0600



4

0,2702

0,2561

0,2340


0,1396

0,1284

0,1175


0,0692

0,0633

0,0584


5

0,2863

0,2699

0,2485

0,2295

0,1451

0,1343

0,1241

0,1160

0,0702

0,0648

0,0605

0,0573


0,2477

0,1271

0,0631

2. The obtained values are averaged for each
porous medium configuration. Figures 710 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 E52690 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 E52690
V2 processor per a core depending on the porous medium configuration and the
time step
Configuration


,
ñ.

,
GFlops DP

, %


0.1

3.548

6.829

28.5

0.01

2.408

7.703

32.1

0.001

2.091

8.323

34.7


0.1

3.928

6.820

28.4

0.01

2.517

7.875

32.8

0.001

2.208

8.731

36.4


0.1

4.143

6.791

28.3

0.01

2.468

7.910

33.0

0.001

2.125

8.902

37.1

We created software for numerical research
of an ideal collisionless gas kinetic problem in the threedimensional 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 quasisteady 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
162915105, 180100343.
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. 100120.
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.: Izdvo TSPI pri mekh.mat. faktu 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. 6976.
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.