A surrogate model for the prediction of permeabilities and flow through porous media

In this contribution we present an approach to generate a data driven surrogate model for the prediction of permeabilities and flow through two dimensional random micro‐heterogeneous materials. The laminar flow is well described by Darcy's law. In order to achieve an efficient computational tool for the generation of the database (up to 103 realizations), needed for the training of the neural networks, we apply a stochastic model based on the Brownian motion. The stationary state of the resulting stochastic model solves the Darcy equation and can be iteratively solved by a Monte Carlo approach applied to a particle simulation. Improved numerical efficiency can be yield by usage of the related transition matrix.


Introduction
Since millennia humanity learns by experimental science about the phenomena of nature, since centuries by theoretical science finding underlying laws. Since some decades scientists perform computer based science e.g. in form of numerical simulations to get predictions of complex phenomena. Since a couple of years the data based science plays a growing importance in the gain of knowledge, [1]. In the most recent years, due to growing computational power availability and improved algorithms, especially supervised and non supervised deep neural networks (DNNs) have emerged in science and engineering for the prediction of complex questions. Supervised trained DNNs can also be used as efficient surrogate models as an alternative to e.g. traditional reduced order models. As a data based method DNNs need large training data sets, see e.g. [2], which can be obtained by simulations using artificial input data. For an overview of applications in material science where deep learning has become a game-changing technique for the prediction of material properties see [3]. As a model system we have chosen the laminar flow through a porous medium at low Reynolds numbers assuming a linear relation between the flux and the pressure gradient. This linear relation is governed by the permeability tensor k which in an inhomogeneous material depends on its location. We assume that the laminar fluid flow can be modeled with Darcy's equation [4]. In order to characterize the effective flow behavior we have to determine the homogenized (effective) permeability tensork. One approach is to find heuristic or semi-analytical approximations for the estimation of the effective permeability as a function of geometrical and statistical properties of the micro-heterogeneous material, for an overview we refer to [5]. Therefore we need another computational inexpensive simulation scheme (that is, from our motivation, a non Finite-Element-based method) for the generation of a sufficiently large training set for the DNN in order to obtain a predictive surrogate model. For an overview of general requirements to achieve surrogate models within a data based approach using large training sets in a deep learning framework see [6], and [7][8][9][10] as examples of applications of trained DNNs as surrogate models in mechanics. Furthermore, e.g. in the context of optimization or parameter identification problems, where a large number of evaluations has to be performed, a fast surrogate model should be designed for efficiency reasons.
The solutions of Darcy's flow problem can be considered as the stationary state of a diffusion process. In order to achieve a computational efficient scheme we use the fact that such processes can be formulated with a particle based atomistic description, the so called Brownian motion of particles, see [11]. This process is driven by a thermal induced stochastic motion of the particles and can in the averaged form be described by the Poisson equation [12]. Therefore the stationary state of the diffusion process is described by the Laplace equation. The advantage of this particle based approach to solve Darcy's problems compared to Finite Element based methods is that no meshing is needed and the computation of the flux fields is computationally simple and fast. The resulting simulation performs random movements of particles to neighbored cells on a equidistant grid and can be realized either by a Monte Carlo method are by a pure algebraic method using the related stochastic transition matrices as described e.g. in [13]. This simulation method is used for fast generation of pairs of random geometries and resulting flux fields which can serve as a training set for machine learning algorithms. Let B ⊂ IR 3 be the body of interest parameterized in x, the surface of B is denoted by ∂B. The later is subdivided into a part ∂B q where the flux across the surface is described and a remaining part ∂B p where the pressure is described. The decomposition of the boundary satisfies the relations ∂B = ∂B q ∪ ∂B p and ∂B q ∩ ∂B p = ∅ . For a fully saturated porous medium the mass conservation law with constant density ρ and porosity appears as div(ρ q) = −f * with the mass source density f * and the flow rate q. This so-called Darcy's law is suitable for the approximate solution of laminar flows. We use the reformulation div q + f = 0 with q = −k · ∇p and f = f * /ρ . (1) Here, the spatially constant permeability k := k * /µ is introduced, where k * , measured in N −1 s −1 , denotes the tensor valued intrinsic permeability and µ the fluid viscosity in Ns/m 2 . The boundary conditions have to be specified where n is the outward unit normal on ∂B, and [p] SI = Pa = N/m 2 , [q] SI = m/s. The Darcy equation is valid for a slow flux under isothermal conditions. In order to get a large number of solutions for different geometries in reasonable time, we use a simple simulation tool based on a stochastic Brownian motion. The simulation generates random geometries and computes the resulting flux under given periodic boundary condition. These pairs of geometry with the resulting flux build the training set.

Brownian motion based on a Monte Carlo approach
In the following we consider a two dimensional micro-heterogeneous structure consisting of two phases, an impermeable phase A and a permeable matrix B. Phase A is given by the union of the inclusions which are defined by a given number of circles with random centers inside the considered domain [0, l] × [0, h] and with radii r ∈ [r min , r max ] varying uniform randomly. The set S of admissible points, i.e. the matrix domain, is given by where i and j are the indices of the equidistant L · H grid points. At the initial state all fluid particles are positioned randomly at admissible points at the left boundary, i.e. at the left column of the discretized region given by j = 1 and i = 1...H and (i, j) ∈ S. Each particle can move potentially to all four directions provided the target position is admissible. In each iteration step the following rules of motion are applied to each particle: 1. Exactly one of the four direct neighbors (left, right, lower, upper) can be visited with probability 1 4 .

This movement is rejected if the chosen neighbor is not in S.
3. If the particle wants to cross a horizontal boundary edge i = 1 or i = H the particle is moved to the corresponding position on the upper/lower side, see Fig. 1 green lines. 4. If the particle wants to cross a vertical boundary (j = 1 or j = L) the particle is moved to a randomly chosen admissible point on the left boundary, see Fig.1

red lines.
If the move is rejected, i.e. the target position is not in S, the particle stays for this step at the current position. An iteration step inside this algorithm can be written in pseudo code as: These transition rules, see Fig. 1, assign a cylindrical topology and define the left edge as inflow and the right one as outflow boundary. It is important to observe that these transition rules keep the number of moving particles constant and that the sum of probabilities for the possible moves of each particle is 1. Fig. 2 shows two snapshots of a Monte Carlo simulation with 1000 × 1000 grid points and 800.000 particles. The most time critical part of the algorithm is the generation of the random numbers. The simulation of 100.000 time steps needs on a standard laptop 10 minutes.
be the probability that a particle located at position J moves to position I. M ∈ IR L·H×L·H is a sparse asymmetric matrix of dimension L · H × L · H, see Fig. 3 as an example. Let consider p I ∈ [0, 1] be the probability that a particle can be found at position I, M p with p ∈ IR L·H is the probability distribution after applying the transition rules given above. Due to the definition of M the vector1 ∈ IR L·H (all entries are 1) is a left eigenvector of M with the eigenvalue 1. Because the 1-norm ||M|| 1 is 1, the number 1 is also the largest eigenvalue of M, see e.g. [14]. Having a left eigenvalue with eigenvalue 1 there must also exist a right eigenvector associated to this eigenvalue. This means there is a vector * p with M * p = * p, the so called stationary state, which can be computed by forward iteration starting with an arbitrary probability distribution p 0 and iterating p n+1 = M p n . The initial probability distribution p 0 at position I is in correspondence to the Monte Carlo approach given by: where m is the number of admissible positions at the left boundary. The computational complexity is nearly the same as in the Monte Carlo approach but due to the usage of fast sparse matrix-vector operations faster by a factor of three: to get a nearly converged stationary state on a grid with L = H = 1000, see Fig. 4, both methods get with 100.000 iteration comparable results. The matrix based implementation needed about three minutes on a standard laptop.

Cascadic stochastic and multi-grid solver
An essential improvement of numerical efficiency can be archived by using a cascade of models which is produced by coarsening the model recursively by a factor of two. This cascade is used by taking recursively the stationary state of coarser level www.gamm-proceedings.com as a starting state on the finer one. This reduces the number of iterations needed on the (computational expensive) largest level to some ten. Because the iterations on coarse level need negligible amount of time the computing time for the same resolution of 1000x1000 is reduced to approximately ten seconds. An even faster convergence we get by applying a multi-grid solver defined by its algebraic components: restriction operator from level j to level j+1 the prolongation the restricted transition matrix the smoothing (for a given k > 0) and the coarse grid solver By construction the restricted matrices M j are also stochastic, I − M jmax +1 ⊗1 has the same eigenvalues as M jmax and is invertable. To achieve a satisfying convergence the smoothing operator must be used with a rather large exponent k > 100. The resulting computational time for the resolution of 1000x1000 reduces to approximately four seconds.

Extraction of flux and effective permeability from the stochastic model
We can identify the pressure as the density which is just the spatial probability distribution p. The probability that a single particle moves from position J = (i, j) to I = (i, j + 1) is given by the number M IJ . The horizontal component of the macroscopic flux q I at position I is the effective amount of particles moving from J to I given by 4e-5 2e-5 0e-5 Setting the dynamic viscosity to 1, we have by Darcy's law: whereq x is the x-component of the averaged flux,q x H is the total flow through the domain measured by the flow through the outflow boundary, √ 2H = L 2H . Fig. 6 shows the stationary states of these configurations. www.gamm-proceedings.com

Conclusion
This work aims to propose an alternative simple approach to the determination of the effective moduli of flow fields. We investigated a stochastic model based on the Brownian motion. The stationary state of the resulting stochastic model solves the Darcy equation and can be iteratively solved by a Monte Carlo approach applied in form of a particle simulation. Improved numerical efficiency could be achieved by the usage of the related transition matrix in form of a forward iteration. Further improvements of the convergence by a cascadic approach and alternatively by applying the multi-grid method are described.