Numerical Simulation of the Propagation of Hydraulic and Natural Fracture Using Dijkstra ’ s Algorithm

Utilization of hydraulic-fracturing technology is dramatically increasing in exploitation of natural gas extraction. However the prediction of the configuration of propagated hydraulic fracture is extremely challenging. This paper presents a numerical method of obtaining the configuration of the propagated hydraulic fracture into discrete natural fracture network system. The method is developed on the basis of weighted fracture which is derived in combination of Dijkstra’s algorithm energy theory and vector method. Numerical results along with experimental data demonstrated that proposed method is capable of predicting the propagated hydraulic fracture configuration reasonably with high computation efficiency. Sensitivity analysis reveals a number of interesting observation results: the shortest path weight value decreases with increasing of fracture density and length, and increases with increasing of the angle between fractures to the maximum principal stress direction. Our method is helpful for evaluating the complexity of the discrete fracture network, to obtain the extension direction of the fracture.


Introduction
Hydraulic fracturing is one of the key technologies in exploitation of natural gas, which plays a major role in enhancing petroleum reserves and daily production especially from unconventional reservoirs [1].The hydraulic fracture and branch fractures could be generated during Hydraulic fracturing, as shown in Figure 1 [2].The hydraulic fracture can provide a high-conductivity path for methane migration in the low permeability reservoir.In recent years, to investigate propagation and re-orientation of hydraulic fracture, lots of research has been conducted [2][3][4][5].
The Finite Element Method (FEM) and Displacement Discontinuity Method (DDM) are two popular methods.Mohammadnejad and Khoei used a fully coupled extended finite element method for hydraulic fracture propagation of porous media [6].Ferté et al. combined extended the finite element method and cohesive zone models to predict crack paths and to compute the crack deflection angle [7].Fu et al. simulated hydraulic fracturing in arbitrary discrete fracture networks by using an explicit coupling simulation strategy [8].However, due to the multiple dynamic physical processes with characteristic length-scales across several orders of magnitude, the time steps of the Finite Element Method (FEM) for simulating fracture propagation are very small.Furthermore, the matrix in physical simulation is not only complex, but also ill-conditioned.Hence, FEM is not suitable to be used to simulate hydraulic fracture propagation.Verde and Ghassemi implemented Fast Multipole Displacement Discontinuity Method to accelerate the solution of large scale coupled fluid flow-geomechanical problems [9].Dong and De Pater used a general displacement discontinuity formulation to deal with the problems of curved and straight-line cracks [10].Behnia et al. studied the boundary element method based on the displacement discontinuity formulation for mixed-mode crack tip propagation of pressurized fractures [11].However, the influence of the coefficient matrix of the boundary element method is not sparse.In addition, with the increase of gird, the time consumed increases rapidly.Hence, the application of DDM in simulating complex fracture network is also limited.In this paper, due to the weaknesses of these two methods for simulating hydraulic fracture propagation in a large scale, a new numerical model based on graph theory, energy theory and vector methods was studied to obtain the final configuration of fracture propagation.A graph is an important mathematical method to describe the special relationship between variables [12].A graph includes some vertices and arcs.The vertices represent the members, and the arc represents the unique relationship between them [12].The shortest path problem is one of the most fundamental problems with widespread applications [13], which is to find the shortest path between two vertices.The shortest path problem's solution can be used in a number of fields.Thus, more researchers have focused on this problem.Shirinivas et al. [14] have given an overview of the graph theory in various fields to some extent and mainly focused on the computer science applications.Klampfer et al. [15] have used mathematical analysis methodologies based on circular graphs to solve a shortest path routing problem in Routing Information Protocol.Sarbu and Valea [16] have proposed that an optimal solution could be obtained based on graph theory for selection of source location and the main path of water transmission.All of the above applications are given with complete graph information, that is, all vertices and arcs are given in this graph.However, parts of vertices and arcs in graph are not given in some special fields.For example, new vertices and arcs are generated with the fracture propagation in discrete fracture networks (DFN) during hydraulic fracturing.This paper conducts exploratory research based on this special field.
In this paper, a new method is presented to obtain the final configuration of fracture propagation with high computation efficiency, which combines Dijkstra's algorithm, energy theory and vector method.Changes in hydraulic fracture propagation are generally attributed to various distribution properties of fractures (such as fracture length, fracture density and dip angle) and simulation conditions (like elastic modulus, normal stress pressure, Poisson's ratio and so on).Hence, different parameters are analyzed, which can show how these parameters influence the propagation of hydraulic fracture.Some preliminaries are given in Section 2. Section 3 discusses how to obtain the configuration of propagated hydraulic fracture using Dijkstra's Algorithm.Parameter analysis is done in Section 4. Section 5 states a discussion and some conclusions are given in Section 6.

Preliminaries
Graph G(V, E) is a collection, which is composed of a vertex set V(G) and edge set E(G).In this paper, the element of E(G) denotes fracture.Definition 1. [12] Make G = (V, E, W) be a weighted graph, where V is the vertex set, E is the edge set, and W is the weight matrix for E. W represents the cost from one vertex u to other vertex v (indicated by ω(u, v)).If e = (u, v), replace ω(u, v) with ω(e).Definition 2. [12] The length of a path p = (e 0 , e 1 , ..., e k ) is the sum of the weights of its constituent edges: (1) Definition 3. [12] The distance from u to v, shown as δ(u, v), is the shortest path with the minimum sum of the weights on the edges in a u − v path.u, v denoted start vertex and sink vertex, respectively.
The following example is given to illustrate the meaning of Definition 1-3, which is shown in  From Definition 2, it is easy to obtain that length(A, D) = ω(A, D) = 9 or length(A, D) = ω(A, B) + ω(B, D) = 5.Hence, the distance from A to D is 5.

Dijkstra's Algorithm
The shortest path problem can be divided into four situations in actual problems, and there are always four kinds of the corresponding algorithms [17]: (1) Dijkstra's algorithm for Single Source Shortest Path, the weight of edge is non-negative [18]; (2) Bellman-Ford's algorithm, it is also for the Single Source Shortest Path, and the weight of edge can be negative.However, the negative weight circuit cannot exist [19]; (3) Shortest Path Fastest Algorithm for Single Source Shortest Path-it is an improved algorithm based on the Bellman-Ford algorithm [20]; (4) Floyd algorithm for Shortest Path among all vertices-the weight of edge can be negative [21].However, the negative weight circuit cannot be existed.For hydraulic fracturing, the weight of fracture is non-negative and the hydraulic fracture propagation belongs to the Single Source Shortest Path problem.Hence, Dijkstra's algorithm is the most suitable method to simulate the hydraulic fracture during propagation of hydraulic fracturing.
Dijkstra's algorithm can be used to find the shortest path between the node and every other [18].Here, V = {v 0 , v 1 , v 2 , ..., v n } is divided into two vertex sets: S = {v 0 } and T = {v 1 , v 2 , ..., v n } .Assume v 0 is the initial vertex, and dist[v i ] is used to store the upper bound of the shortest path distance from v 0 to v i .If the final shortest path from v 0 to v i is determined, denoted S = S ∪ v i , T = T \ v i .The algorithm repeatedly selects the node v i from T with the minimum dist[v i ], and puts it in set S. The main steps for Dijkstra's algorithm can be described as follows: Step Step 3.
Step 3. Set S = S ∪ v i , T = T \ v i , and S[v i ] = 1.If T = ∅, stop, else go to Step 4. Step An example is given to illustrate the working of the Dijsktra's algorithm, which is shown in Figure 3.To begin with, the weight matrix should be computed.Assume that A is the start vertex, the weight of the minimum-weight paths from A to every other vertex is +∞, and A to A is 0.

Model Description
In this section, Dijsktra's algorithm is used to simulate the hydraulic fracture propagation in the DFN.DFN is given under the meter scale, which has two groups of fractures.During the hydraulic fracture propagation, two kinds of energies are used to describe different fracture propagation: mechanical energy and interface energy.The mechanical energy is used to keep natural fracture open, and interface energy is used to create new fractures.Combined the mechanical energy and interface energy, the fracture weighted formula is successfully deduced.

Discrete Fracture Network
An example of DFN is a two-dimensional horizontal section of a rock with 4 m length, 4 m width, which is shown in Figure 4. DFN is composed of two groups of random fractures, and the blue lines represent natural fractures.The fracture number is 250, the length of a fracture is 0.2 m.The angle between the two groups of fractures and the dip angle (which also means the angle between one set fractures to the maximum principal stress direction) is 60 • .Figure 4 also shows that the DFN is not fully connected from one side to the other.During hydraulic fracturing, new fractures are the only way to link natural fractures.Therefore, new fractures should be an important factor in Dijsktra's algorithm.

The Fracture Intersection Point
There are many fracture intersection points in Figure 4.During hydraulic fracturing, when fracture fluid flows through a fracture intersection, fracture fluid will enter into two fractures.Hence, the new vertices generated from fracture intersections have a significant influence on DFN connection, and these new vertices must be found out.Many mathematic methods to detect the fracture intersection have been studied, such as parameter equation and vector method [22].Vector method is a great approach to this problem.e 1 = (u 1 , v 1 ) and e 2 = (u 2 , v 2 ) represent two fractures, respectively.Defining the two-dimensional vector cross product: Vector quantities r = (u 3 , v 3 ) and s = (u 4 , v 4 ) represent the other two fractures, respectively.Supposing the two line segments run from e 1 to e 1 + r and from e 2 to e 2 + s.The main steps are as follows [23]: (3) If r × s = 0 and (e 2 − e 1 ) × r = 0, then the two lines are parallel and non-intersecting.(4) If r × s = 0 and 0 ≤ t, u ≤ 1, the two line segments meet at the point e 1 + tr = e 2 + us.
(5) Otherwise, the two line segments are not parallel but intersect.
Using this method, all of the fracture intersection points can be found in DFN, which are shown in Figure 5, and marked by the red points.

The Fracture Weighted Formula
During hydraulic fracturing, energy can well explain the fracture propagation, which includes the elastic modulus, normal stress, Poisson's ratio and so on [1].Hence, the fracture weighted formula with energy will be innovatively deduced in this section.Energy can be divided into two categories: mechanical energy and interface energy.The mechanical energy is for keeping natural fractures open, and it is related to the fracture length.Interface energy is for creating new fractures, which is related to the fracture length [24].The energies during fracture propagation are described by weighted value in the shortest path problem.Considering fracture length, fracture type and so on, the weighted values of fracture are different from each other.Mechanical energy should be considered in natural fractures.However, for new fractures, both mechanical energy and interface energy should be considered during computing the weight values of new fractures [24].The derivation processes of the fracture weighted formulas with energy are described as follows: Assuming that fracture propagation length is c and the normal stress is P. The expression is [24]: where S[σ xx , σ xy , σ yy ] is in situ stress, and σ xx and σ yy are the x and y direction stress.σ xy is shear stress and β is the dip angle, and the schematic describing the fracture and stresses is shown in Figure 6.Mechanical energy U M is the sum of the bearing system potential energy U A and elastic strain energy U E [25]: where G is the energy per unit area or length [25], c is the length of fracture, λ(c) is elastic compliance, and ν is Poisson's ratio.Generating Equations ( 4) and ( 5) into Equation ( 3), it is easy to obtain In another case, new fracture length is c, then, interface energy U S is: where γ is the surface energy per unit area or length.
During fracture propagation, the weighted value of natural fracture is computed only by considering mechanical energy U M , which can be better understood by the energy of keeping the fractures open during gas extraction.The weighted value of new fractures is computed by combining mechanical energy U M and interface energy U S .

Determination of the Hydraulic Fracture
Assume that the direction of the maximum principal stress is horizontal.Predetermine that the new fracture in DFN has connected its downstream vertex with its upstream vertex.For example, vector (u, v) represents the start vertex (indicated by u(u x , u y )) and the final point (indicated by v(v x , v y )) of a branch.If u x is greater than v x , vertex u is marked as an upstream vertex.In addition, vertex v is made as a downstream vertex.On the contrary, vertex v is marked as upstream vertex.All the vertices in Figure 4 and intersections in Figure 5 are the elements of set V.
Select N vertices from the set V as source points.If the abscissa of these N vertices is smaller than other vertices', these N vertices are denoted as initial vertices.In addition, found N vertices are also regarded as terminal points, if the abscissa of these N vertices are larger than other vertices'.Assume N = 12.The initial vertices and terminal points can be obtained in this DFN by comparing their abscissa values, which are shown in Figure 7.The vertices marked by red '*' are the start vertices (indicated by v 1 , v 2 , ..., v 12 ), and the terminal points are marked by a blue '*' (indicated by u 1 , u 2 , ..., u 12 ).Basing on Definition 3 of Section 2, the DFN is a collection of three sets, indicated by G(V, E, W).All endpoints of the fracture in Figure 3 and intersection points in Figure 4 are the element of set V = v 1 , v 2 , .., v n , where v j+12 = u j (j = 1, 2, ..., 12), and V 1 = v 2 5, v 2 6, .., vn describe the intersection points.The element of set E is composed of the fracture in DFN.Set W is computed by Equations ( 6) and ( 7) in Section 3.3.Take in situ stress S[σ xx , σ xy , σ yy ] = [6.0× 10 7 Pa, 1.0 × 10 7 Pa, 5.0 × 10 7 Pa] in Equation ( 2), Poisson's ratio ν = 0.35 in Equations ( 5) and ( 6), γ = 1.0 × 10 6 J/m in Equation (7).Based on Definition 2.2-2.3 in Section 2, p ij = (v i , u j ) presents the weight from v i to u j (i, j = 1, 2, ..., 12), and δ(u, v) represents the shortest path with the minimum sum of the fracture weights.Using Dijkstra's algorithm with energy, the hydraulic fracture of DFN can be obtained.The main steps of algorithm are as follows: Step 1. Set x = 1.
Step 3. Select the vertex v i in T, such that dist Step 4.
In addition, the algorithm flow chart is shown in Figure 8. From the above steps, it can be shown that the potential distances are now recorded as updated edge weights so that no additional memory is required.Through combining the vector method and geological parameter with Dijkstra's algorithm, the hydraulic fracture of DFN is obtained, which is marked by a red color in Figure 9.

Experiment and Comparison
The sample is sampled in the Jiaoshiba shale gas field, Sichuan Basin, and the sample's detail is shown in Figure 10.We drilled a hole with diameter 10 mm and depth 50 mm in the center of sample.During the hydraulic fracturing test, the axial stress and confining stress are 30 MPa and 20 MPa, respectively, and the injection rate is 30 ml/min.For comparison, the parameters in Section 4.1 are equal to the experimental conditions.The failure morphology of the sample after the experiment is shown in Figure 11a.It is worth noting that the sample has obvious horizontal bedding faces.The direction of natural fractures in the numerical model is along these horizontal bedding faces.As demonstrated in Figure 11a, it can be shown that the hydraulic fracture propagates through the whole sample along the direction of maximum horizontal in situ stress.Figure 11b shows the simulation results using the method above in Section 4.1.It can be shown that the hydraulic fracture, obtained by the method in this paper, is similar to the the failure morphology of the sample.

Parameter Sensitivity Analysis
During fracture propagation, several factors can influence the propagation of hydraulic fracture in a discrete fracture network system, such as intensity of discrete fracture (fracture number, and marked with n), fracture length (marked with l), and dip angle (marked with α).In this section, the weighted value of hydraulic fracture is used to explore how these parameters influence the hydraulic fracture propagation.Here, the DFN is a two-dimensional horizontal section of a rock with 4 m length and 4 m width.Assuming that the DFN is static, partially connected and randomly distributed.The maximum principal stress direction keeps in line with the x-axis, and many groups of fractures are repeated computing for eliminating the influence of randomness.Figures 12-15 show quantitative comparison of the behavior of the models based on different parameters.Firstly, the influence of fracture number on the weight value of hydraulic fracture is shown in Figure 12.The trends of scatter distribution and average distribution are consistent with each other, and the weight value of hydraulic fracture decreases with the increasing of fracture number.That also means that the greater the fracture, the better for the hydraulic fracture propagation.Then, Figure 13 shows that the weight value of hydraulic fracture decreases with increasing of fracture length.It also indicates that as the fracture length is growing, the weight value of the hydraulic fracture is reducing.Longer fracture length will benefit the hydraulic fracture propagation.Finally, as demonstrated in Figure 14, it has little influence on the weight value of hydraulic fracture for different values of α when the angle between two groups of fracture keeps orthogonal.Hence, the influence of different values of α on the weight value of hydraulic fracture will also be studied under two groups of fractures are not orthogonal.Figure 15 shows the weight value of the hydraulic fracture firstly decreases and then increases with the increase of α when the angle between two groups of fractures are 20 • .That being said, the farther away from the maximum stress direction, the larger the weight value of the hydraulic fracture.

Discussion
The operating parameters are effective only when the geological parameters are favorable for fracturing treatment.Therefore, the evaluation of the geological parameters before hydraulic fracturing is crucial to hydraulic fracturing treatment in shale formation.The effects of geological parameters are investigated, including fracture density, fracture length, and dip angle (the angle between one group of fracture and the maximum principal stress direction).By analyzing the sensitivity of DFN with various parameters, it can be concluded that the hydraulic fracture weight value decreases with increase in the intensity of fracture and fracture length, and increases with increasing of the angle between fractures to the maximum principal stress direction.It means that the fracture density and fracture length are significant to the DFN connectivity.

Conclusions
Increases in natural gas extraction are being driven by rising energy demands, mandates for cleaner burning fuels, and the economics of energy use [26].Directional drilling and hydraulic fracturing technologies are allowing expanded natural gas extraction from organic-rich shales [27], while the propagation of the hydraulic fracture is the crux of hydraulic fracturing, and can provide a high-conductivity path for methane migration in the low permeability reservoir.Accompanying the benefits of such extraction [28] are public concerns about hydraulic fracturing that are ubiquitous but lack a strong scientific foundation.In this paper, Dijkstra's algorithm, energy theory and vector methods were combined to simulate the hydraulic fracture propagation in DFN during hydraulic fracturing.We successfully give the shortest path weighted value formula of fracture propagation, and obtained the hydraulic fracture in DFN during hydraulic fracturing.The direction of the hydraulic fracture is also simulated, which plays a key role in identifying the spatial position of the well bore during hydraulic fracturing.Numerical results with experimental data show that the hydraulic fracture obtained by our method is similar to the failure morphology of the sample (see Figure 11).
This study has played a very important role in the hydraulic fracturing field, including the following aspects: (1) obtaining the extension direction of the fracture; through comparing the energies of fracture propagation in different directions, the direction of the minimum energy fracture propagation is the most likely direction of the hydraulic fracture; (2) evaluating the complexity of the DFN.Through Dijkstra's algorithm, several paths can be obtained whose energies are smaller.If the positions of these paths are very near, it is shown that the hydraulic fracture is a great advantage and it is harder to get a complex fracture network.On the contrary, if the positions of these paths are very far away, it is shown that the hydraulic fracture lacks strength and it is easy to get a complex fracture network; and (3) through computing the weight value of hydraulic fracture, the injected energy during hydraulic fracturing can be obtained, which is strongly related with pore pressure and fluid viscosity.During hydraulic fracturing treatment, the pore pressure decreases with the distance far away from the well bore, and this factor has been considered while computing fracture weight value.

Figure 2 .
Figure 2. The flowchart of a graph.

Figure 4 .
Figure 4.An example of discrete fracture network.

Figure 5 .
Figure 5.The fracture intersection points of discrete fracture networks.

Figure 7 .
Figure 7.The start vertex and terminal vertex of discrete fracture network.

Figure 8 .
Figure 8.The flow chart of the completed algorithm.

Figure 9 .
Figure 9.The hydraulic fracture of the discrete fracture network.

Figure 10 .
Figure 10.The dimensions of the sample.

Figure 11 .
Figure 11.(a) the failure morphology of the sample; (b) the simulated fracture.

Figure 12 .Figure 13 .Figure 14 .
Figure 12.The influence of fracture number (n) on the hydraulic fracture weighted value.Here, α = 30 • , l = 0.2 m, and the angle between two groups fracture remains orthogonal.Eight kinds of fracture numbers (175, 200, 225, 250, 275, 300, 325 and 350) have been analyzed, with 150 groups of randomly generated fractures for each kind of fracture number.(a) scatter distribution of the the hydraulic fracture weighted value with different values of n; (b) average of the hydraulic fracture weighted value with different values of n.

Figure 15 .
Figure 15.The influence of dip angle (α) on the hydraulic fracture weighted value when the angle between two groups of fracture are 20 • .Here, n = 250 and l = 0.2 m.Ten kinds of α (0 • , 9 • , 18 • , 27 • , 36 • , 45 • , 54 • , 63 • , 72 • and 81 • ) were simulated, with 150 groups of fractures randomly generated for each kind of α.(a) scatter distribution of the hydraulic fracture weighted value with different values of α; (b) average of the hydraulic fracture weighted value with different values of α.