CC BY
7
Вычислительные технологии
Том 22, № 3, 2017
Numerical simulation of the percolation threshold of the electric resistivity
E.P. Shurina1,2'*, M.I. Epov1,2, A.Yu. Kutischeva1,2 Novosibirsk State Technical University, Novosibirsk, 630073, Russia
2Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, 630090, Russia *E-mail: shurina@online.sinor.ru
Theoretical estimation, laboratory experiments and numerical simulations are the typical approaches to homogenization. In this paper, we propose a numerical method for obtaining effective resistivity. This approach is based on the simulation of the electrostatic field in heterogeneous media. We employ heterogeneous multi-scale finite element method (FE-HMM) to solve this problem in complex media (i. e. we consider materials with high contrast micro inclusions). FE-HMM is a two stage method. On the first stage a macro-scale solver is selected and missing macroscopic data are estimated. This makes it possible to study percolation (the flow of the electric current through the complex media) without any simplifications of the model.
In this paper we consider a 3D problem that requires to determine effective resistivity of an object with cuprum plate inclusions necessary for the study of the percolation processes. Developed FE-HMM algorithm allows us to obtain effective scalar characteristics (e.g., thermal conductivity, electrical resistivity etc.) of the materials containing high-contrast inclusions. Developed algorithms are implemented in C++ for parallel architectures and are adaptable for the use on modern supercomputers. In order to obtain effective resistivity of the samples with inclusions with different shape and physical characteristics, we carried out a number of computational experiments. Our study indicates that the surface square of each inclusion is inversely related to the percolation threshold of the material with the inclusions. Our studies also indicate that if inclusions are localized in the subdomain of the media, a simplification of the model to the uniform distribution will be unacceptable.
Keywords: effective resistivity, heterogeneous multiscale finite element method (FE-HMM), heterogeneous media, percolation threshold.
Introduction
Heterogeneous media are the macroscopic systems consisting of homogeneous components, delimited by boundary surfaces (e. g., rocks, composite materials). Such media can be described by the system "matrix/inclusion" [1]. In the heterogeneous media studies it is necessary to determine media that are multi-scale (i.e. the object size significantly exceeds the size of the inclusions) and multi-physics (i. e. inclusions and matrix have contrasting characteristics).
© ICT SB RAS, 2017
Studies of the properties of the composite materials show that if the concentration of the inclusions in an object exceeds critical value, effective properties of this medium may abruptly change [2, 3]. This process is called percolation [4]. The percolation theory was proposed in [4]. According to this theory, the system "matrix/inclusion" is considered as a finite or infinite system where material properties are concentrated inside the nodes [5]. Neighbouring points with similar properties are grouped into clusters that can be analyzed as single structures [6]. The lowest concentration of the conductive inclusions, at which a continuous path from one border to the other is formed, defines the threshold of percolation.
At present the percolation theory is most frequently used in the studies of the filtration processes [7- 10] and of the conductivity of the composite materials (e. g. polymer composites containing nanotubes [11, 12] or nano-plates [13, 14]). In order to study the influence of various factors on the percolation threshold, various approaches based on special theoretical model [15, 16], laboratory experiments [17] and numerical simulations of heterogeneous media are used. However, such approaches do not permit to determine the effective characteristics of the media near the percolation threshold. Numerical methods, based on the mathematical models of the processes in heterogeneous media make it possible to establish threshold concentrations of inclusions, as well as to calculate effective characteristics of the media considered.
There exist various numerical methods for the study of the scalar effective properties that can be described by the stationary model (e. g. permeability, electrical conductivity, thermal conductivity, etc.). In most cases, in order to solve such problems, Galerkin formulation of the finite element method (FEM) in [18] is used. Due to the geometrical complexity of the problems considered, the finite element method in its classical formulation cannot be used, as it will lead to the significant increase of the dimensionality of the problem. There exist a number of modifications of the FEM for the solution of the multi-scale and multi-physics problems [19]. One of these specialized methods designed to investigate properties of the materials that have multi-scale physical and geometrical properties is the heterogeneous multi-scale finite elements method (FE-HMM [20-22]). It is based on the multi-scale finite element method (MsFEM [23]) and generalized finite element method (GFEM [24-27]). FE-HMM belongs to a growing class of multi-scale methods based on the decomposition of the entire space of solutions into the sum of two subspaces: the coarse space responsible for the effective properties of the media, and the fine space where all small scale features can be considered (e.g., pores, fibrous structures, inclusions, cracks, etc.) [2, 28]. In the FE-HMM polyhedrons can be used as finite elements in the spacial discretization, unlike the MsFEM, where parallelepipeds are the only elements that can be used on the coarse grid in the three-dimensional case. We can apply the FE-HMM to solve such problems in domains with complex geometries.
1. The mathematical model
In this paper we consider a three-dimensional problem of computing effective electrical resistivity of the heterogeneous media under the direct current. The mathematical model of the process considered is a homogeneous elliptic equation with Dirichlet and Neumann boundary conditions.
-V ■ (p£ (x)-1 Vu£ (x)) =0 in Q,
(1)
where p£ (x) is electrical resistivity (Ohm ■ m), Q = {x G R31 (x G Qi) V (x G Q0)} is the solution domain which consists of contrast micro inclusions Q1 and a matrix Q0.
We suppose that Q is a cylinder with a boundary dQ = r1 U r2 U r3 U dQjn, where r1, r2, r3 are external boundaries, dQin is the boundary between inclusions Q1 and matrix Q0. We set up inhomogeneous Dirichlet boundary conditions on the boundaries r1 and r3 (the top and bottom faces) that correspond to the potential difference. On the boundary of r2 (side faces) homogeneous Neumann conditions are set up
p£(x) du£ (x)
-1 du£ (x) d n
= 0, [u£ (x)U =0,
dQi
d n
(2)
r2
= 0, u£ (x)|ri = 0, u£ (x)|r3 = 1,
where n is the external normal to the boundary, [.] is the jump operator.
Therefore it is necessary to solve homogeneous elliptic equation (1) with the boundary conditions (2) on the domain Q with contrast micro inclusions with different geometries.
2. The heterogeneous multiscale finite element method
The FE-HMM is a generalization of the classical grid multi-scale method and the mesh-free generalized finite element method. By analogy with the MsFEM in the FE-HMM the solution is defined as a linear combination of the non-polynomial shape functions defined on the finite supports. However, these functions are not explicitly constructed on the entire macro element (the coarse element). Instead, on each macro element a quadrature rule and integration points are selected and the neighbourhoods are formed by joining fine grid elements around each integration point, with the following restrictions: each neighbourhood is located strictly inside the macro element and the integration points cannot be on the boundary of the appropriate neighbourhood. Non-polynomial shape functions are constructed as a decomposition on the basis of the fine grid in the neighbourhood of each integration point, but not on the entire macro element. Therefore, the coupling of the non-polynomial shape functions is based on the quadrature rule. Therefore, the FE-HMM is a "quasi mesh-free" method. This approach complicates assembly procedure of the global solution and makes it impossible to use standard techniques typically used in the FEM and the MsFEM.
Our approach permits to work on the coarse grid level without additional builds of the free form finite supports in order to solve problems with complex internal structures. This approach can be seen as an approximate virtual element method [29]. Non-polynomial shape functions are constructed and processed independently that allows for efficient parallel implementation of the algorithm.
2.1. The multiscale shape functions
Let S (Q) be a conformal decomposition of Q C R3 into polyhedral elements Kp such that
NK
IJ Kf = Q, where Nk is the total number of elements in S (Q). Usually, S (Q) is called ¿=1
the coarse grid, and Kp is called the macro element. We will call [Kp, F, E, A} the finite element [30], where Kp C Q is a simply-connected polyhedron; F C H1 (Q) is the space of non-polynomial multi-scale shape functions , i = 1,p, (Fig. 1) defined in Kp and
Fig. 1. An example of the shape function: a — the macro element; b — the shape function in plane section
reflecting an inner microstructure of the media, dimF = p; E = F* is the subspace of degrees of freedom, dim E = p; A is the quadrature rule. Let the quadrature rule A for the polyhedron Kp be
nK
f{x)dx = (xi) , (3)
K p 1
where uk is the number of the integration points, u and xi G R3 are weights and nodes for the numerical integration, respectively. In the FE-HMM it is necessary that the quadrature rule A satisfies the following conditions:
1) Vxi, 1 = 1^ : xi /QKP;
2) the order of the quadrature rule A is constant for all elements Kp G S (Q);
3) let polyhedrons KX,K2 G S (Q) have a common face or edge (i.e. the elements K and K2 are neighbours), then the integration points xi1 G Kx, = 1,u^i and xi2 G K2, l2 = 1,uk2 are symmetrical about the common face or edge, respectively.
In order to construct the non-polynomial multi-scale shape functions , i = 1,p the finite element grid should be generated according to the following algorithm:
1) an adaptive nonuniform conformal simplicial decomposition (Kp) = {k} of all macro elements Kp is generated;
2) around each integration point xi ,1 = 1, uk we construct a sphere Sph (xi) C Kp with the center at xi; the sphere radius is chosen based on the physical characteristics of the problem being solved and the geometrical properties of internal structures;
3) the polyhedron I (x{) = G ^ (Kp) : fc C Sph (xi)| is created taking into account
the following condition xi G dl (xi ).
Therefore a polyhedron I (xi), I = 1, nx consisting of the elements of the simplex decomposition (Kp) is constructed around each integration point xi, I = 1, uk .
We define linear interpolation shape functions ^^ (x) on ^ (Kp) and discrete subspace
Hi (Kp) = span i-à*
(Kp) = span {-f :i=1,PK} C H (Kp) , (4)
where PK is the number of degrees of freedom of the simplicial decomposition (K).
The non-polynomial multi-scale shape functions p, i = 1,p are introduced as an expansion on the basis of the space (4)
p K
<pF (x) = ^ ql€ (x) V = T^,
r=1
where p is the number of degrees of freedom of the macro element Kp.
By using solution in each neighbourhood I (xl), I = 1,n? we will find the expansion weight
-V ■ {p£ (x)-1 VV?P ) = 0 in I (Xl), ^ (x)\dI{xi) = (x) ,
where rqf (x) is the solution of the reduced 2D problem on each face Es C R2 of the polyhedron I (xl), and , i = 1,p is the non-polynomial multi-scale shape functions.
V2 D ■ {p£ (x)-1 V2DV? (x)) = 0 in Es, (x)\dEs = if (x) ,
where i = 1,p is the number of degrees of freedom of the macro element Kp,\Js Es = dl (xl) and Vs,r : Es H Er = 0, £E (x) is a function, that is defined in Kp and is equal to 1 on the z-th node and 0 on the other nodes of Kp. The function (x) is linear on edges of Kp.
2.2. The variational formulation
Consider Hilbert spaces:
rl^n\ _ L. . V7„. (T \\3
H1
(Q) = {u G L2(Q) : Vu G (i2(n))3} (u(x),v(x)) = u(x)v(x) dQ,
n
F01(n) = G H 1(Q): = 0} .
For the elliptic boundary value problem (1), (2), the variational formulation takes the following form:
find u£ G + uo (dty such that
J p£ (x)-1 Vu£(x) • Vv(x) dQ = 0 Vv G Hl(Q). n
The solution of the global multi-scale problem is a linear combination of the multi-scale non-polynomial shape functions
u£'h (x)= £ it bK <pK (?),
Kp&(Q) ¿=1
where S (Q) is the decomposition of the domain Q, p is the number of degrees of freedom of the macro element Kp, and <fKP is the non-polynomial multi-scale shape function defined in Kp.
The variational formulation of the FE-HMM for the elliptic boundary value problem (1), (2) has the following form:
find us'h E Vh(tt)+u0(dfi) such that
Z^»-' V^-W-M^o W EVh(n), * E n
n
The corresponding system of linear algebraic equations (SLAE) in matrix form can be expressed as follows:
^global u _ global
where Aglobal is the global matrix of the SLAE. This matrix is obtained by assembling local matrices AK P>local for VKP gS (fi):
A— = / (✓(*»-' '* EK(5)
K p
Each element (5) is determined according to the following quadrature rule (3):
nK
AKP,local = (p£(*i) )-' VvK(*) • VVf (*).
1=1
3. Numerical experiments
In order to determine the percolation threshold (the critical volume concentration of inclusions), a series of computational experiments for parallelepiped samples (30x30x 100 mm) with inclusions of different form, concentration and physical properties (see table) were carried out. Electrical resistivity of the matrix in all experiments was set to pmatrix = 1.4 Ohm-m that corresponds to the quartz sand saturated NaCl.
Classification of inclusions
Type Electrical resisti- Volume con-
No. of inclusions vity of inclusions centration of Size of inclusions, mm
pincl, Ohm-m inclusions, %
1.1 Spheres 10-8 10-44 d = 3.85
2.1.1 Plates 10-1 0.7-10.5 a = 0.3; b = 5 L = 10
2.1.2 10-2
2.1.3 10-4
2.1.4 10-6
2.1.5 10-8
2.2 10-8 1.5-13.7 a = 0.6; b = 5 L = 10
2.3 10-8 1.4-6.6 a = 0.3; b = 5 L = 20
3.1 Deformed plates (curve centered at 90°) 10-8 1.8-7.1 a = 0.3; b = 5 L = 10
3.2 10-8 1.3-11.9 a = 0.6; b = 5 L = 10
3.3 10-8 1.6-4.9 a = 0.3;6 = 5; L = 20
Fig. 2. Examples of objects that contain plates 0.6 x 5 x 10 mm (type 2.2) with 11 % concentration: a — uniform distribution; b — exponential distribution for A = 0.4
In the samples, inclusions are randomly distributed according to the three-dimensional uniform (Fig. 2, a) and exponential Exp (A) distributions (Fig. 2, b).
The calculation of the effective resistivity according to the formula [2] was verified in physical and computational experiments [31, 28].
Peff = (||J||/||grad(«)||)
-1
(
In
| J|
1
\
|grad (u) |2
where J is a current density; u is a scalar potential.
3.1. The influence of the electrical properties of inclusions on the percolation threshold
Fig. 3 shows the calculation results of the effective resistivity for the media, that contain uniformly distributed plates 0.3x5x 10 mm (the types 2.1.1-2.1.5) (Fig. 2, a). The resistivity of the inclusions varies from 10-1 to 10-8 Ohm • m.
The percolation threshold occurs when the concentration of the inclusions is ~ 5 % (conductive plates 2.1.3-2.1.5). In samples with small conductive inclusions (the plates 2.1.1, 2.1.2) percolation does not occur. Therefore, all further experiments will be performed under the condition of the high contrast of resistivity of the matrix and inclusions. This result is consistent with the experiments [32].
nJ G
'B tj 0)
o
W
Pi Ji
O
§
>
1.00E-01 1.00E-01 1.00E-03 1.00E-05 1.00E-07
■ ____ 7"!- —•-•- -•
\ v \ \ \
% •
2 4 6
Volume concentration of inclusions, %
•Plates 2.1.1 Plates 2.1.2 Plates 2.1.3 Plates 2.1.4 Plates 2.1.5
Fig. 3. The variation in the effective resistivity for inclusions with different resistivity
3.2. The influence of the inclusions shape on the percolation threshold
The investigation of the heterogeneous materials shows dependence of the effective resistivity of the media on the shape of the particles, even under the low concentrations of the inclusions [33]. We now consider a number of samples when inclusions have equal volumes but square surfaces of the inclusions are different. We consider spheres (type 1.1), plates 0.6 x 5 x 10 mm (type 2.2) and 0.3 x 5 x 20 mm (type 2.3), deformed plates 0.6 x 5 x 10 mm (type 3.2) and 0.3 x 5 x 20 mm (type 3.3). Inclusions are uniformly distributed (Fig. 2, a). Electrical resistivity of inclusions is set to pinci = 10-8 Ohm ■ m.
The results of the computational experiments (Fig. 4) show that for the spherical inclusions the percolation threshold is achieved at the conducting inclusion concentration of ~ 40 %, while for the plate inclusions the percolation occurs already at 3.5-4.5 %. It should be noted that all inclusion types have equal volumes but the ratio of the surface square of the inclusions is in the range from 2 to 4.
3.3. The influence of an inclusion distribution on the percolation threshold
A significant part of the natural and artificial materials have random nonuniform distribution of inclusions in volume (e.g. viscous media are characterized by a settlement of particles). However, when considering heterogeneous media with inclusions it is assumed that the particles in an object are uniformly distributed over the entire volume [34, 35]. To assess the influence of the distribution model on the percolation threshold, we consider samples with plate inclusions 0.6 x 5 x 10 mm (type 2.2), exponentially distributed as follows Exp (A) so that density distribution increases from the top side to the bottom (Fig. 2, b). The resistivity of the inclusions is pinci = 10-8 Ohm ■ m.
In case considered exponential distribution can be seen as a model where particles settle out in viscous media. If the value of A increases, the degree of localization will increase. In the lower part of object there are more inclusions than there are particles in the upper part of the object. Value A = 0.01 corresponds approximately to the uniform distribution of inclusions in the sample. From Fig. 5 it follows that with increasing A the threshold concentration increases. For A = 0.4 the threshold concentration is higher in 2 than the
itf c
u G
'S à
Tj o
> £
'S -S .ti
ta vi
m S
l.e+01 1.e+00 l.e-01 l.e-02 l.e-03 l.e-04 l.e-05 l.e-06 l.e-07
0 20 40 60
Volume concentration of inclusions, %
¿Su___
s
1 * __
1__ I
1 1
1 1
•i
m
l.E+01
1.E+00
l.E-01
l.E-02
l.E-03
l.E-04
l.E-05
l.E-06
l.E-07 O.E-OO
mïïm
- \ \
: 1 V \ : \ \ \
\ \ \1 Xv \
• » ; \ \ • V
— Spheres 1.1
- Plates 2.2 Plates 2.3
Deformed plates 3.2 Deformed plates 3.3
2.E+04 4.E+04 6.E-04 Total square surface of inclusions, mm2
b
a
Fig. 4. The variation in the effective resistivity of the media with increasing concentration of inclusions. Inclusions have equal volumes but different square surface: a — dependence of the effective resistivity of the volume concentration of inclusions, b — dependence of the effective resistivity of the total square surfaces of inclusions
Fig. 5. The threshold concentration of inclusions for different values of A
threshold concentration of the uniform distribution. Therefore for the simulation of materials where localization of inclusions is high enough it is necessary to consider the distribution type, since the use of the uniform distribution model leads to the significant errors in the results.
3.4. The mixtures of different inclusions
In addition to the different types of the distributions of inclusion, there exist mixtures of particles with different geometries (Fig. 6, a) and with a variety of contrasting electrical and physical properties (Fig. 6, b).
3.4-1- The mixtures of inclusions with various electrical resistivities
Consider the objects (Fig. 6, a) containing uniformly distributed plate inclusions 0.3 x 5 x 10 mm (type 2.1). Electrical resistivity of the half of the inclusions is p1 (Ohm • m), and for the rest of the inclusions electrical resistivity equals p2 (Ohm • m). We define the following mixtures of the particles:
1) all inclusions have high conductivity p1 = p2 = 10-8 Ohm• m (type 2.1.5);
2) all inclusions have low conductivity p1 = p2 = 10-2 Ohm• m (type 2.1.2);
3) a mixture of the low and high conductivity particles p1 = 10-2 Ohm • m (type 2.1.2), P2 = 10-8 Ohm• m (type 2.1.5);
Fig. 6. Examples of samples with various inclusions with concentration of 5%: a — particles with different electrical resistivity; b — particles with different geometry
ca g
o .
'B g
O J3
^ O 0»
<D
o
EH w
1.00E+01 1.00E+00 L00E-01 1.00E-02 L00E-03 1.00E-04 1.00E-05 1.00E-06 L00E-07
"Ni V
\ \ \ \ "4
\ * *
t i \ :
\ X
Plates 2.1.5 Plates 2.1.2 Plates 2.1.2 and 2.1.5 Plates 2.1.4 and 2.1.5
0 2 4 6 8 10 12 Volume concentration of inclusions, %
Fig. 7. The change of the effective electrical resistivity of the media with mixtures of the inclusions with various resistivity
® £
o .
'S a
O JH
JS O
ID
> P
O £ W
1.00e+01
1.00E+00
1.00E-01
1.00E-02
1.00E-03
1.00E-04
1.00E-05
1.00E-06
1.00E-07
0 2 4 6 8 10 12 14 16
Volume concentration of inclusions, %
L
* i.
"J K X
"t
* i i i >
i v1,
• f »1 • 1
j
—■a— Deformed plates 3.2 —Mixture of 3.2 and 2.2 - Plates 2.2
Fig. 8. The change of the effective resistivity for media with mixtures of inclusions with different geometries
4) a mixture of the high conductivity particles p\ = 10 6 Ohm■ m (type 2.1.4), p2 = 10-8 Ohm■ m (type 2.1.5).
Figure (Fig. 7) illustrates that in case of the mixture with particles of low and high conductivity the percolation paths are formed when high conductivity inclusions reach critical concentration (triangular markers in Fig. 7). If media contains a mixture of high conductivity inclusions with different resistivities, two abrupt changes of the effective resistivity (the rhomb markers in Fig. 7) will be observed.
3.4-2. A mixture of inclusions with different geometry
In order to examine mixtures of inclusions with different geometries we consider samples that contain deformed and undeformed plates 0.6 x 5 x 10 mm (there are types 3.2 and 2.2) in equal parts. Simulation results (Fig. 8) showed that the use of this mixture gives a small change of the percolation threshold compared to the samples with only one type of inclusions.
Conclusion
In the study of the heterogeneous materials consisting of a matrix and inclusions, one of the most important problems is to simultaneously determine the threshold concentration and effective properties of the media near the percolation threshold. Developed FE-HMM algorithms allow us to obtain effective scalar characteristics (e.g., thermal conductivity, electrical resistivity, etc.) of the materials containing high contrast inclusions. The use of the free-form finite (polyhedron) elements also makes it possible to work with the samples that have a complex external boundary. Developed algorithms were implemented in C++ for parallel architectures and the code has been adapted to run on modern supercomputers.
In order to obtain effective resistivity of the samples with inclusions with different shape and physical characteristics, a number of computational experiments have been conducted. Our results can be summarized as follows:
1) if the media are saturated by low contrast inclusions (the ratio of the matrix resistivity and the inclusions resistivity is less than 102) the percolation threshold will not be achieved;
2) the surface square of each inclusion is inversely related to the percolation threshold of the material with the inclusions (this result is in agreement with the geometric assumptions that objects with larger surface square form clusters at the lower volume concentrations);
3) if inclusions in the media are localized in the subdomain of the media, a simplification of the model to the uniform distribution will be unacceptable;
4) for the mixture of inclusions with different resistivities several step changes of the effective resistivity of the media may occur (there are several percolation thresholds);
5) if the samples contain deformed plates (type 3.1-3.3.) effective resistivity of the material will be less than the effective resistivity of the materials with the undeformed plates, however the percolation threshold will vary slightly.
Acknowledgements. This work is supported by the Grant FPMM No. 43 (2016). References
[1] Dulneev, G.N., Zarichnyak, Yu.P. The thermal conductivity of mixtures and composite materials. Leningrad: Energiya, 1974. 264 p. (In Russ.)
Дульнев Г.Н., Заричняк Ю.П. Теплопроводность смесей и композитных материалов. Л.: Энергия, 1974. 264 с.
[2] Epov, M.I., Shurina, E.P., Artem'ev, M.K. Numerical homogenization of electrical characteristics for the media with contrasting small-scale inclusions // Doklady RAN. 2012. Vol. 442, No. 1. P. 118-120. (In Russ.)
Эпов М.И., Шурина Э.П., Артемьев М.К Численная гомогенизация электрических характеристик сред с контрастными мелкомасштабными включениями // Докл. РАН. 2012. Т. 442, № 1. С. 118-120.
[3] Durmaz, S. A numerical study on the effective thermal conductivity of composite materials. IZMIR, 2004. 240 p.
[4] Broadbent, S., Hammersley, J. Percolation processes I. Crystals and mazes // Proc. of the Cambridge Philosophical Society. 1957. Vol. 53. P. 629-641.
[5] Stauffer, D., Aharony A. Introduction to percolation theory. London: Taylor & Francis, 2003. 180 p.
[6] Bunde, A., Kantelhardt, J.W. Theory and applications / A. Bunde, J.W. Kantelhardt (Eds). Diffusion and Conduction in Percolation Systems. Berlin; Heidelberg: Springer, 2005. P. 895-914.
[7] Su, B., Fujimitsu, Y., Li, C. Application of pseudo-percolation threshold theory in evaluation of low-porosity and low-permeability reservoir // J. of Petroleum and Gas Exploration Res. 2012. Vol. 2, No. 7. P. 125-134.
[8] Mukhopadhyay, S., Sahimi, M. Calculation of the effective permeabilities of field-scale porous media // Chemical Eng. Sci. 2000. Vol. 55. P. 4495-4513.
[9] Friedman, S.P., Jones, S.B. Measurement and approximate critical path analysis of the pore-scale-induced anisotropy factor of an unsaturated porous medium // Water Resources Res. 2001. Vol. 37, No. 12. P. 2929-2942.
[10] Rodriguez-Valverde, M., Ramon-Torregrosa, P.J., Cabrerizo-Vilchez, M. Estimation of percolation threshold of acid-etched titanium surfaces using Minkowski functionals // Microscopy: Science, Technology, Applications and Education. 2010. Vol. 3. P. 1978-1983.
[11] Bauhofer, W., Kovacs, J.Z. A review and analysis of electrical percolation in carbon nanotube polymer composites // Composites Science and Technology. 2009. Vol. 69. P. 1486-1498.
[12] Dang, Z.-M., Shehzad, K., Zha, J.-W., Mujahid, A., Hussain, T., Nie, J. Complementary percolation characteristics of carbon fillers based electrically percolative thermoplastic elastomer composites // Composites Science and Technology. 2011. Vol. 72. P. 28-35.
[13] Syurik, J., Alyabyeva, N., Alekseev, A., Ageev, O.A. AFM-based model of percolation in graphene-based polymer nanocomposites // Composites Science and Technology. 2014. Vol. 95. P. 38-43.
[14] Ma, P.C., Li, J., Sze, C.W., Kai, T.C., Tang, B.Z., Kim, J.-K. Percolation threshold of polymer nanocomposites containing graphite nanoplatelets and carbon nanotubes // Proc. of the ICCM Intern. Conf. on Composite Materials, Kyoto, Japan, 2007. P. 1-8.
[15] El-Tantawy, F., Deghaidy, F. Effect of iron oxide on vulcanization kinetics and electrical conductance of butyl rubber composites // Polymer Intern. 2000. Vol. 49. P. 1371-1376.
[16] Yi, Y., Tawerghi, E. Geometric percolation thresholds of interpenetrating plates in three-dimensional space // Physical Review. 2009. Vol. 79. P. 041134-1-041134-6.
[17] Bilotti, E., Zhang, H., Deng, H., Zhang, R., Fu, Q., Peijs, T. Controlling the dynamic percolation of carbon nanotube based conductive polymer composites by addition of secondary nanofillers: The effect on electrical conductivity and tuneable sensing behaviour // Composites Science and Technology. 2013. Vol. 74. P. 85-90.
[18] Bondeson, A., Rylander, T., Ingelstrom, P. Computational Electromagnetics. Berlin: Springer, 2000. 231 p.
[19] Reed, W.H., Hill, T.R. Triangular mesh methods for the neutron transport equation // Tech. Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973. P. 1-23.
[20] E, W., Engquist, B. The heterogeneous multiscale methods // Communications in Math. Sci. 2003. Vol. 1, No. 1. P. 87-132.
[21] E, W., Ming, P., Zhang, P. Analysis of the heterogeneous multiscale method for ellpiptic homogenization problems // J. of the American Math. Soc. 2003. Vol. 8. P. 121-156.
[22] Abdulle, A., Schwab, C. Heterogeneous multiscale FEM for diffusion problem on rough surfaces // SIAM Multiscale Modeling and Simulation. 2005. Vol. 3, No. 1. P. 18-19.
[23] Hou, T., Wu, X.-H. A multiscale finite element method for elliptic problems in composite materials and porous media // J. of Comput. Phys. 1997. Vol. 134. P. 169-189.
[24] Babuska, I., Banerjee, U., Osborn, J.E. Generalized finite element methods — main ideas, results and perspective // Intern. J. of Comput. Methods. 2004. Vol. 1, No. 1. P. 67-103.
[25] Babuska, I., Caloz, G., Osborn, J.E. Special Finite Element Methods for a Class of Second Order Elliptic Problems with Rough Coefficients // SIAM J. Numer. Anal. 1994. Vol. 31, No. 4. P. 945-981.
[26] Duarte, C.A., Oden, J. An hp adaptive method using clouds // Computer Methods in Applied Mechanics and Engineering 1996. Vol. 139. P. 237-262.
[27] Duarte, C.A., Babuska, I., Oden, J. Generalized finite element methods for three dimensional structural mechanics problems // Comput. Struct. 2000. Vol. 77. P. 215-232.
[28] Kutischeva, A.Yu. The solution of elliptic problems with contrast coefficients the heterogeneous multiscale finite element method // Sbornik Nauchnykh Trudov NGTU. 2013. Vol. 1, No. 71. P. 52-58. (In Russ.)
Кутищева А.Ю. Решение эллиптических краевых задач с контрастными коэффициентами гетерогенным многомасштабным методом конечных элементов // Сб. тр. НГТУ. 2013. Т. 1, № 71. С. 52-58.
[29] Beirao da Veiga, L., Brezzi, F., Cangiani, A., Manzini, G., Marini, L., Russo, A.
Basic principles of virtual element methods // Math. Models and Methods in Appl. Sci. 2013. Vol. 23, No. 1. P. 199-214.
[30] Brenner, S., Scott, R. The mathematical theory of finite element methods. New York: Springer-Verlag, 2008. 400 p.
[31] Shurina, E.P., Epov, M.I., Kutischeva, A.Yu. Determination of the effective coefficient resistance for materials with microinclusions using heterogeneous multiscale finite element method // Fizicheskaya Mezomekhanika. 2016. Vol. 19, No. 3. P. 93-102. (In Russ.) Шурина Е.П., Эпов М.И., Кутищева А.Ю. Определение эффективного коэффициента электрического сопротивления материалов с микровключениями гетерогенным многомасштабным методом конечных элементов // Физ. мезомеханика. 2016. Т. 19, №. 3. P. 93-102.
[32] Heagy, L.J., Oldenburg, D.W., Chen, J. Where does the proppant go? // GeoConvention 2014: FOCUS, 2014. 7 p.
[33] Kalaitzidou, K., Fukushima, H., Drzal, L.T. A Route for polymer nanocomposites with engineered electrical conductivity and percolation threshold // Materials. 2012. Vol. 3. P. 1089-1103.
[34] Yu, Y., Song, G., Sun, L. Determinant role of tunneling resistance in electrical conductivity of polymer composites reinforced by well dispersed carbon nanotubes // J. Appl. Phys. 2010. Vol. 108, No. 8. P. 084319.
[35] Hautot, S., Tarits, P. Effective electrical conductivity of 3-D heterogeneous porous media // Geophys. Res. Letters. 2002. Vol. 29, No. 14. P. 14-1-14-4.
Received 18 April 2016 Received in revised form 22 April 2017
