A Numerical Study on Hydraulic Fracturing Problems via the Proper Generalized Decomposition Method

The hydraulic fracturing is a nonlinear, fluid-solid coupling and transient problem, in most cases it is always time-consuming to simulate this process numerically. In recent years, although many numerical methods were proposed to settle this problem, most of them still require a large amount of computer resources. Thus it is a high demand to develop more efficient numerical approaches to achieve the real-time monitoring of the fracture geometry during the hydraulic fracturing treatment. In this study, a reduced order modeling technique namely Proper Generalized Decomposition (PGD), is applied to accelerate the simulations of the transient, non-linear coupled system of hydraulic fracturing problem, to match this extremely tight response time constraint. The separability of the solution in space and time dimensions is studied for a simplified model problem. The solid and fluid equations are coupled explicitly by inverting the solid discrete problem, and a simple iterative procedure to handle the non-linear characteristic of the hydraulic fracturing problem is proposed in this work. Numeral validation illustrates that the results of PGD match well with these of standard finite element method in terms o f fracture opening and fluid pressure in the hydro-fracture. Moreover, after the off-line calculations, the numerical results can be obtained in real time.


Introduction
Hydraulic fracturing, also termed as hydro-fracking, is a very important technology in the petroleum industry because it can obviously stimulate more production of oil and gas wells [Yew and Weng (2014); Wang, Shahvali and Su (2015)]. Recently, with advances in these technologies such as horizontal well drilling and segregated completion, hydraulic fracturing has become an essential part of commercially developed unconventional hydrocarbon reservoirs, such as shale gas, coal-bed methane, deep geothermal resources and tight hydrocarbon reservoirs. All the time, the determination of the fracture geometry and propagation path in hydraulic fracturing has been a research hotpot in petroleum engineering [Nolte and Economides (2000)]. However, the numerical simulation of hydrofracking problem is a time-consuming task due to the strongly non-linear coupling between the equations for the fluid flow in the hydro-fracture and rock deformation. To handle this hydro-mechanical coupled problem, various numerical methods have been successfully put forward in recent years [Hunsweck, Shen and Lew (2013); Yuan, Zheng, Moghanloo et al. (2017); Kumar, Camilleri and Brewer (2016); McClure, Babazadeh, Shiozawa et al. (2016)]. These numerical methods have been successfully applied in the simulation of hydraulic fracturing, which are involved in finite difference method (FDM), finite element method (FEM), extended finite element method (XFEM), discrete element method (DEM), and numerical manifold method (NMM), etc. For the first two methods, their solutions will depend on the mesh size near the hydraulically-driven fractures because of the stress singularity in the vicinity of crack tips [Pogacnik, Elsworth, O'Sullivan et al. (2016); Yan and Zheng (2017); Shimizu, Murata and Ishida (2011);Douillet-Grellier, Pramanik, Pan et al. (2016)]. DEM offers more advantages than other numerical methods in the simulation of discontinuous problems, but most of the simulation time is occupied in the explicit iteration of small time step and contact judgment between two adjacent elements [Pogacnik, Elsworth, O'Sullivan et al. (2016)]. For XFEM, the calculated crack can freely propagate that is independent on mesh sizes because additional enrichment functions are introduced into the standard shape functions. However, the enrichment process will take up most of central processing unit (CPU) time, and thus it is difficult to simulate complex fracture networks [Shimizu, Murata and Ishida (2011);Haddad, Du and Vidal-Gilbert (2017)]. Therefore, these mentioned methods are time-consuming in the simulation of hydrofracking problem because we have to deal with the same above-mentioned difficulties: non-linear, transient and coupled. It will take so much time, normally several hours to days, to simulate the fracture propagation process on a computer, therefore it still remains a great challenge to predict the fracture geometry in real time during the fracking operations [Shimizu, Murata and Ishida (2011);Douillet-Grellier, Pramanik, Pan et al. (2016); Miehe and Mauthe (2016); Haddad and Sepehrnoori (2015); Adachi, Siebrits, Peirce et al. (2007);Ji, Settari and Sullivan (2009)]. In order to efficiently solve this difficulty, the new idea is to solve the non-linear, transient and coupled partial differential equations (PDEs) via model order reduction method (ROM) such as Proper Orthogonal Decomposition method (POD) or Reduced Basis (RB) method [Secchi and Schrefler (2012) ;Chinesta, Keunings and Leygue (2013); Chinesta, Ladeveze and Cueto (2011);González, Alfaro, Quesada et al. (2015); Henneron and Clenet (2015)]. Model order reduction aims to reduce the computational complexity by means of a reduction of the associated state space dimension or degrees of freedom. The Proper Generalized Decomposition (PGD) method, a kind of novel ROM method, is firstly introduced by French mathematician Pierre Ladevèze et al. [Chinesta, Ladeveze and Cueto (2011);Ladevèze, Passieux and Néron (2010)]. In the PGD framework, the solutions are represented by several independently separated functions. It means that PGD method can successfully separate the time and space domain for transient PDEs, and thus the original time-dependent problem is divided into 1D time problem and 2D or 3D space problem. And then alternating fixed strategies are adopted to solve the non-linear and coupled problem with a relatively low computational cost in an enrichment process. Recently, PGD method has been successfully applied in many challenging issues involved in various fields of science and engineering [Signorini, Zlotnik and Díez (2017); Aguado, Huerta, Chinesta et al. (2015)], and it is proved to be a more efficient and powerful algorithm [Aguado, Huerta, Chinesta et al. (2015); Ammar (2010); Zlotnik, Díez, Modesto et al. (2015); Tamellini, Le Maitre and Nouy (2014)]. In this paper, the PGD method is firstly applied in hydraulic fracturing problems to achieve the fast simulation of the fracture geometry. Firstly, the singular valued decomposition (SVD) method is used to linearize the cubic term of PDEs, then the fast PGD solver is constructed to separate the time domain and space domain to reduce the associated dimensions, and finally get the fracture geometry in a real-time mode. Verification test is carried out by the comparison of results between PGD and FEM. This new PGD solver will be much more useful and attractive for simulating the hydraulic fracturing problem within the oil and gas industry in the future.

Problem statement
To simply the construction of PGD-ROM, the following assumptions are made without loss of generality: (1) the hydraulic fracturing treatment is considered to be a two-dimensional plane strain problem; (2) there is almost no fluid leak-off normal to the fracture wall because of the low matrix permeability; (3) the injection fluid is Newtonian fluid, and thus the viscosity is assumed as a constant in hydraulic fracturing; (4) it is an isothermal process during the hydraulic fracturing, and thus the fracturing fluid viscosity is independent on the temperature.

Governing equations
The governing equations of hydro-fracking problem are composed of two coupled equations: the classical stress equilibrium equation representing rock deformation, and the transient fluid flow equation in the hydraulically-driven fracture [Yew and Weng (2014)].

Stress equilibrium equation
Based on the classical theory of linear elasticity, the solid part of hydro-fracking PDEs consists in stress equilibrium equation and proper boundary conditions [Yew and Weng (2014); Nolte and Economides (2000)] finding a displacement function u(x, y) taking values in a 2D domain of Ω , as shown in Fig. 1, such that (4) u = 0 at the four corner points (5) (2)-(5) denote its boundary conditions; σ denotes the second-order Cauchy stress tensor; u denotes the displacement tensor vector; p(s, t) denotes the fluid pressure scalar acted on the fracture wall Γ f ; σ H and σ h respectively denote the maximum and minimum horizontal principal stress in the far field; n t and n f respectively denote the normal unit vector, which is perpendicular to the outer boundary Γ out = Γ L Γ R Γ T Γ B and inner boundary Γ f ; s denotes the space variable representing crack positions; t denotes the injection time. Under the condition of elastic deformation, the constitutive equation between stress tensor and strain tensor can be represented as [Yew and Weng (2014); Hunsweck, Shen and Lew (2013)]: where σ e denotes the effective stress tensor; ε denotes the second-order strain tensor; D denotes a fourth-order elastic stiffness tensor. If isotropic rock is under plane strain conditions, D is a 3 × 3 square matrix: where E denotes Young's modulus, and ν denotes Poisson's ratio. Under small-deformation assumption, there is a linear relationship between strain tensor and displacement vector [Yew and Weng (2014); Fjar, Holt, Raaen et al. (2008)], where ∇ s denotes the symmetric operator. According to the Terzaghi-Biot effective stress principle, Cauchy stress can be divided into two parts [Fjar, Holt, Raaen et al. (2008)]: where p p denotes the pore pressure acted on rock skeletons; α denotes Biot elastic coefficient, α ∈ [0, 1].

Fluid flow equation
According to the lubrication theory in fluid mechanics, the continuity equation of a planar Poiseuille flow between two parallel plates satisfies [Yew and Weng (2014); Nolte and Economides (2000)]: w 3 12µ ∂p ∂n = 0 on two crack tips (13) where w denotes the fracture width; µ denotes the viscosity of fracturing fluid; Q 0 denotes the source term, i.e., injection point; k denotes the fracture permeability; Ω f denotes the crack surface. It should be mentioned that Eq. (10) is a non-linear equation because of the cubic term w 3 , i.e., the coefficient of the second derivative. According to Eq. (13), it shows that the fluid flow equation satisfies the Neumann boundary condition at the two crack tips. In Eq. (12), the superscripts "+" and "-" denote the two fracture surfaces of Ω f . According to Eq. (12), there is a certain sparse geometry matrix A, such that For each row in matrix A, the only two non-zero entities are equal to 1 or -1, and others are all equal to zero. Eq. (14) smoothly builds the connection between displacement field u in domain Ω and fracture width field w in domain Ω f .

Weak forms
According to the variational principles for stress equilibrium equation and fluid flow equation, the trial solution spaces of displacement and fluid pressure fields are respectively defined as follows: where H 1 denotes Hilbert space.
In the above definitions, the two sets (15) and (16) satisfy with Dirichlet boundary conditions of solid part Eq. (1) and fluid part Eq. (10). After multiplying the weighted functions δ u and δ p on two sides of Eqs. (1) and (10), and integral by parts, their corresponding weak forms can be expressed as Eqs. (17) and (18). Its statement is as follows: find (u, p) ∈ u × p such that for all (δu, δp) where σ 0 is equal to σ H or σ h on the outer boundary Γ out [Hunsweck, Shen and Lew (2013); Adachi, Siebrits, Peirce et al. (2007)].
3 Numerical methods

FEM formulation
Assume that N u i and N p i respectively denote the nodal displacement and fluid pressure shape functions, u i and p i respectively denote the nodal displacement and fluid pressure value, N u and N p respectively denote the vector of nodal displacement and fluid shape functions, u and p respectively denote the vector of nodal displacement and fluid pressure. Then the displacement and fluid pressure field are respectively expressed as:

The discretization of stress field equation
According to Galerkin FEM method, the discrete form of Eq. (17) is represented as follows: where K s denotes the stiffness matrix, K s = Similar to the above-mentioned process of solid part, substituting δp = N p , p = N pp , u = N uũ into the above weak form, and then the discretization scheme of its weak form is obtained as follows: whereu denotes the first order time derivative ofũ; u n+1 is the nodal displacement in domain Ω at time n + 1; To discretize the discretization of time derivative term in Eq. (21), backward finite difference is adopted to approximate the time derivative of displacement u in Eq. (21).
where n denotes the n-th time layer. The cubic term w 3 in the integral of K f (w 3 ) makes Eq. (24) to be a nonlinear equation.
Finally the fully coupled solid-fluid forms between Eq. (19) and Eq. (23) can be expressed in the block matrix form: where u n+1 and u n respectively denotes the displacement at the time n + 1 and n.

Preliminary treatment
According to the discrete form of solid part Eq. (19), the displacement function u can be expressed as: where b = K −1 s f s , and B = K −1 s M s . Substituting the above equation into Eq. (14) and Eq. (23), we obtain For simplicity, the subscript "f " is ignored in the next section without any ambiguity. Eq. (26) indicates us that once the separated forms of pressure solution of Eq. (28) is obtained, it is easy to get the PGD solutions of displacement and crack opening according to Eqs.
(25) and (26). Here, it should be mentioned that the vector b and matrix B in Eq. (25) are only the function of space variable. Thus, the next step is only to seek for the PGD solution of fluid part by means of inverting the discrete form of solid part.

The separated representation of the coupled system
The separated representations of PGD of fluid problem are as follows [Chinesta, Keunings and Leygue (2013); Chinesta, Ladeveze and Cueto (2011)]: where m denotes the modes of separated terms; X and T are unknowns at current enrichment step m; p m−1 is the first m − 1 terms of PGD form, which is already computed.
Here again for simplicity, the subscript symbol "m" is neglected instead of X, T without any ambiguity. Substituting Eq. (29) into Eq. (28) yields, Because of the non-linear term K w 3 , here SVD method is used to separate the matrix into space domain and time domain [Chinesta, Keunings and Leygue (2013); Chinesta, Ladeveze and Cueto (2011)].
where [K l ] ij = Ωf N i θ l (x)N j dx. Now our expected separated forms of fluid problem are formulated as:

Alternating direction iteration
At enrichment step m, both functions X and T are unknown at the current enrichment step. Their product X · T makes the equation nonlinear. Thus herein alternating direction iteration is utilized to handle the non-linear problem. This kind of iterative scheme includes two steps for each iteration step p: (1) given T p−1 m (t), compute X p m (x); (2) given the justcomputed X p m (x), compute T p m (t). At the beginning of iteration, a random vector guess T 0 m is specified. The stopping criterion of alternating direction iteration is defined as follows [Chinesta, Keunings and Leygue (2013); Chinesta, Ladeveze and Cueto (2011); Signorini, Zlotnik and Díez (2017); Aguado, Huerta, ]: CMES, vol.122, no.2, pp.703-720, 2020 where · is the L2-norm, and ε is a specified tolerance. Once the computed X p m (x) and T p m (t) satisfy the criterion, both of them will be assigned to X m (x) and T m (t) respectively. In following, the superscript of iterative step p is also again neglected for simplicity without any ambiguity.
• Given T , Calculate X Since all functions of t are known, T is multiplied on both sides of Eq. (34) and then carry out the integration in time domain. It means that given T , find X such that In the above expression, the following coefficients associated with 1D integral over Ω time are defined as: Therefore the compact form of our final seeking function X is expressed as: The strong form of the Boundary Value Problem (BVP) is an algebra equation of X, it can be numerically solved by means of any appropriable numerical methods.
• Given X, Calculate T Similar to the above procedure, since all functions of x are known, X T is multiplied on both sides of Eq. (34). It means that given X, find T such that In the above expression, the following coefficients dependent on space variable x are defined as: Thus the compact form of Eq. (44) is formulated as: The strong form of the Initial Value Problem (IVP) is a first-order ordinary differential equation (ODE) of T, it can be numerically solved by means of the fourth-order Runge-Kutta (RK4) method. Assume initial condition T (0) = 0, and use RK4 scheme, we get:

Linearized strategy
In the previous section, the matrix K w 3 is non-linear, therefore Picard iterative method (i.e., fixed point method) is used to solve and linearize the discrete system of Eq. (27). The corresponding iterative scheme is [Yew and Weng (2014)]: where δ denotes the Picard iteration number. When it satisfies the stopping criterion: p δ+1 − p δ < , the final p will be obtained. strain) + 1D fluid equation. Without loss of generality, the same method can be extended into the case of 3D. Because the fracture width is very small, the fluid equation can be viewed as a 1D problem along the fracture length direction. Because the strain in the zdirection (fracture height) is very small, the stress equilibrium equation can be viewed as a 2D problem (Plane strain) [Yew and Weng (2014); Zlotnik, Díez, Modesto et al. (2015)].  results of PGD match well with the results of FEM, which verify that PGD method has good numerical accuracy and is a reliable method, as shown in Fig. 4. Further, the calculation time via PGD algorithm is much less than that via FEM algorithm, shown in Tab. 3 (a computer with 16 GB memory and Intel 3.7 GHz CPU). It is observed that the average computational cost of PGD is approximate 10% of that FEM, which means that the computational efficiency of PGD is approximate ten times of that FEM. Therefore, PGD method is much faster than FEM method in simulating the hydraulic fracturing problem [González, Alfaro, Quesada et al. (2015);

Convergence analysis
The errors of fluid pressure in the hydro-fracture and fracture opening at different mesh sizes are respectively calculated. The corresponding results are shown in Figs. 5 and 6, respectively. The error is defined as follows: where the reference solution is defined as the result that calculated on the finest mesh size; X h denotes the pressure or fracture opening solution with different mesh sizes; and X ref denotes the reference solution. As shown in the two figures, their solutions are stable because of the linear relationship between the errors and mesh size. Hence, Results show a very good convergence and accuracy towards the full finite element solution of the problem in terms of fracture opening and fluid pressure in the fracture.  [Yew and Weng (2014); Nolte and Economides (2000); Modesto, Zlotnik and Huerta (2015); Wang, Zhou, Ding et al. (2015)]. In future work, some key engineering and mechanical parameters, such as fluid viscosity, injection rate, Young's modulus, Poisson's ratio, and far-field stress, will be considered into the current PGD formulation, and try to solve the parametrized hydro-fracking equations.

Conclusions
1. In this work, PGD technique is applied to efficiently solve the transient, non-linear and coupled hydro-fracking problem in petroleum engineering. 2. Several numerical examples of hydro-fracking problem show that PGD solver can be smoothly used in most cases with good accuracy. Importantly, the PGD solution is much faster and more efficient than the full finite element solution. 3. The PGD strategy in this paper is quite different from the standard incremental time integration schemes. The number of non-linear iterations of the alternating direction algorithm rarely exceeds ten, while the modes N is often a few tens.