The problem of optimal control of the hydraulic fracturing process is set as the optimization problem of finding the input parameters vector for the fracture propagation model that provides the minimization of objective functions. Various objective functions based on the output data of the fracture propagation model are considered. The problem is formulated within the framework of multivariate and multiobjective optimization method, which is based on the combined features of fracture propagation model and Genetic Algorithm. Plane radial model of fracture propagation caused by Herschel-Bulkley fluid injection is used to establish relationships between input parameters and fracture growth. The potential of the proposed method for control of hydraulic fracturing process is demonstrated by application to hydrocarbon reservoirs of shallow bedding. Results show that the proposed methods for optimal control of hydraulic fracturing process play the important role in maximization of the volume of mined hydrocarbon with significant decrease of the costs for hydraulic fracturing execution.
The goal of the optimization of the hydraulic fracturing process is the maximization of gas and oil production by maximization of the fractured reservoir rock volume. The optimization of hydraulic fracturing is based on the model of fracture propagation in elastic media caused by viscous fluid injection. Input parameters of fracture propagation model are: the surface of the cavity in infinite elastic media; pumping pressure of fluid that causes fracture initiation and propagation (or pumping schedule for fluid of given rheology); elastic media parameters. Output characteristics of the model are: the fracture surface; fracture width distribution; speed of fracture front propagation. Calculation of output characteristics using the input parameters is the direct problem of fracture propagation. By solving it one can predict the geometry of forming fracture, volume of hydrocarbon mined from the fracture, calculate costs of this process, etc. The inverse problem is to find the vector of input parameters of fracture propagation model that satisfies the given performance criteria of fracture formation and propagation processes. Optimal control of hydraulic fracturing process consists in the solution of the inverse problem. In this problem it is required to find the parameters of rheological laws for fluid, pumping schedule (e.g. time dependence of fluid injection rate), conditions of fracture initiation (geometry of the cavity, its orientation against in-situ stresses of elastic media) that satisfy the needed location of incipient fracture, linearity of fracture propagation trajectory, uniformity of fracture width distribution along the trajectory, no sharp bends along the fracture trajectory, minimal costs for hydraulic fracturing, maximal volume of mined hydrocarbon.
In [1] the optimization of hydraulic fracturing process is carried out for the PKN model. The model describes straight fracture propagation from the linear source. It is supposed that the fracture is of constant height much lesser than the fracture length. On this assumption the change of the parameters along fracture height is negligible and the rock deformation can be considered as plain strain state independently in each vertical section. In PKN model fracture toughness is not taken into account. It is considered that the hydraulic fracturing fluid fills the whole of the fracture, fracture tip mechanics is not considered. Filtration leakoff to the rock is taken into account. Fracture geometry is considered as a function of the following parameters that influence the hydraulic fracturing process: hydraulic fracturing fluid viscosity , injection rate of fracturing fluid in( ), injection time , proppant concentration and fracture half-length frac. Fracture height ℎ frac and width are calculated solving coupling relationships based on fracture geometry and material balance. Optimization problem for hydraulic fracturing process consists in maximization of total production over 10 years with bound and design constraints. To solve the optimization problem the algorithm INTEMOB (INTElligent Moving Object) is used. The main concept of INTEMOB is based on the concepts of Genetic Algorithm, simplex-method and EVOP-algorithm (EVolutionary OPeration).
In [2] the optimization of hydraulic fracturing process is carried out using Pseudo-3D (P-3D) model of fracture propagation. P-3D model differs from PKN model by taking into account variability of the parameters along fracture height that influences fracture width. Problem statement and solution method for optimization problem are the same as in [1]. The results of [1] and [2] show that the suggested method for solving of the stated optimization problems plays the important role in hydraulic fracturing process improvement. Only a 12 % compromise with the production over 10 years saves about 44 % of the treatment cost.
In [3] the optimization of hydraulic fracturing process using 3D model is presented. Barnett Shale gas field is considered as the object of investigation. It is characterized by 7 rock layers and 4 sets of joints per rock layer. The goal of the optimization of the hydraulic fracturing procedure is maximizing the fractured reservoir rock volume which results in the maximization of gas production. For three dimensional modeling, analysis and post processing FEM simulator ANSYS is used. For fracturing process modeling the program multiPlas is used. Optimization tool optiSLang is used for calibration of more than 200 parameters of the reservoir production simulator. The main result of modeling is 3D fissured reservoir rock. In the calibration process, the physical parameters are updated until time and location of seismic fracture measurements show reasonable agreement with the simulation results. Correlation analysis of optiSLang identifies the main reservoir parameter to additionally calibrate the mechanism of how the fracturing design parameter effects the fracture growth. Among the wide range of the varied parameters the most important are layer thickness and its location, elastic properties of the rock, hydraulic parameters, bedding plane of the shale, the location and frequency of the fractures in this plane. These parameters are compared (correlated) with reservoir and well test data. The advanced functionality of the software supports taking into account some additional physical effects such as thermal effects. Calibrated model was used to predict the gas production rate. The predicted gas production rate from the calibrated model showed very good agreement to the real production rate and much better agreement than the estimated production rates with the help of seismic fracture measurements only. With help of the optimization design an increase of gas production of 25 % was possible with just an optimized well position in the reservoir.
In [4] optimization pump flow rates scenarios has been performed for the particular wellbore situated in Iranian Sand Stone Reservoir. Pump flow rates scenario is characterized by two parameters in this study: pump rate that is assumed to be constant and time of pumping that is replaced by the fracture half length. These two parameters are varied to maximize wellbore productivity over a period of one year. The parameters are varied independently and only one parametric problem of minimization arises. So no optimization algorithm is needed. But the feature of the paper is the fact that a few different models have been applied: pseudo three dimensional model, PKN, KGD and radial ones.
In [5] 2-D hydraulic fracturing model GEOS-2D, is used to simulate dynamic fracture propagation within a pre-existing facture network. Instead of integrating physical models and economic models to maximize net present value as the objective function, or estimating the total production of the wellbore over any time period the authors of [5] focus on physical criteria, that is, the optimal hydraulic fracture propagation under uncertain natural conditions. The fractal dimension of the connected fractures can be derived from the post-fracturing network simulated by GEOS-2D to represent the network density and connectivity. Therefore, the fractal dimension is chosen as the objective function to optimize the hydraulic fracturing well design. BOBYQA, a derivative-free nonlinear optimization algorithm, is applied in [5] to drive a global search on the modeled response surface.
In the present paper the problem of optimal control of the hydraulic fracturing process is formulated as the optimization problem of finding the input parameters vector for the fracture propagation model that provides the minimization of several objective functions. Various objective functions based on the output data of the fracture propagation model are considered. The problem is formulated within the framework of multivariate and multiobjective optimization method, which is based on the combined features of fracture propagation model and Genetic Algorithm. Plane radial model of fracture propagation caused by Herschel-Bulkley fluid injection is used to establish relationships between input parameters and fracture growth. The capability of the proposed method for control of hydraulic fracturing process is demonstrated by application to hydrocarbon reservoirs of shallow bedding. Results show that the proposed methods for optimal control of hydraulic fracturing process play the important role in maximization of the volume of mined hydrocarbon with significant decrease of the costs for hydraulic fracturing execution. 2
Direct problem of hydraulic fracture propagation The construction of the methods for optimal control of hydraulic fracturing is carried out using the plane radial model of fracture propagation caused by Herschel-Bulkley fluid injection. 2.1
Hydraulic fracture propagation is simulated by the classical penny shaped model [6,7,8] improved by more complex model of slurry flow inside the fracture. Newtonian fluid model is replaced by more general Hershel-Bulkley model [9] and the variation of proppant concentration is taken into account by convection equation added into the model. The geometrical concept of the penny-shaped fracture model is presented in Fig. 1. Rock deformation under axisymmetric pressure load in fracture is given by the integral relation ( ) =
8 ′ ︁∫ frac⎛ ∫︁
⎝ 0 ︀√ 2 net( ) − 2√︀ 2 − 2 ⎞
1 − 2 ⎠ , ′ = , net = − min, where is absolute pressure of the slurry of fracturing fluid and proppant flow, frac = frac( ) is fracture front position defined from the well-known criterion of brittle fracture propagation
Slurry flow in the fracture is described by continuity equation of fluid phase where ︁∫ frac 0 = √ 2 frac
net( ) ︀√ f2rac − 2
= . ( + 2 1 ( ) + 1 2
,
)
(
= 0, 2 +1 net = −2 ︂( 2 + 1 )︂ + ︂( 4 + 2 )︂ 0 + 1 , where , and 0 are slurry parameters. In the considered model the rheology of Herschel-Bulkley fluid is taken into account. In Herschel-Bulkley fluid shear stress is expressed from yield stress 0, power law fluid rheology consistency coefficient , power law fluid rheology behavior index and the shear rate ˙ by the formula
Volume concentrations of the fluid and the proppant for all possible
and are connected by the relation
Maron-Pierce [10] relationship
(
(, ) + (, ) = 1. ( ) = (0) 1 − * ︂( )︂ −2
, ( w, ) = in( ) ( w, ) = in( ) frac − fluid > 0. where * is the critical concentration of proppant. According to the experimental study performed by Mueller [11], this relationship can be used to calculate viscosity of suspensions of solid particles. Yield stress 0 and behaviour index correspond to the used hydraulic fracture fluid.
On the wellbore the boundary condition for the slurry pumping rate and for one of the two concentrations, e.g. for proppant, is set. The fluid front fluid is supposed to lag the fracture front frac On fluid front fluid Stefan condition fluid = ( fluid , ) =
( fluid ) 2 fluid ( fluid ) and the condition of zero absolute pressure
net( fluid , ) = − min are set. Net pressure net on the interval from fluid front fluid to fracture front frac is also supposed to be − derive the equation of slurry balance at any time moment
min. From (
(, 0) = 0, w ≤ ≤ 0. 2.2
The hydraulic fracturing process lasts seconds. In Fig. 2 the flowchart of numerical solution of incorporated to the model subproblems is presented. Let it be the fracture with front is defined by the increment frac at the timestep . Fracture front at ( + 1) timestep
frac to fracture front frac. Fracture increment
frac has the fixed value for all timesteps of the problem. The fluid front
position is set flui+d1 < fr+ac1. For the given fluid front position the equations (
min. Note that the time needed for the given fracture increment
(
m
s = 0: ps(r) = p0
W s(r) =W (ps(r)): Eq.(
k = 0: Δtk = Δt0
Qk(r) =Q (W n(r),W s(r), QLn(r), Qinn,δinn,Δtk ), uk : Eq.(
no ps+1(r) =ω p%s(r) +(1−ω ) ps(r) yes yes Δts = Δtk, Qs(r) = Qk(r)
Eq.(
no Rfmlu+id1 = R (Rfmluid )
yes Rfnlu+i1d = Rfluid m The primary goal of hydraulic fracturing is to increase the productivity of a well by superimposing a highly conductive structure onto the formation. One of the objective functions that are maximized in present paper is productivity index (PI) . This index comes from the linear relation between the production rate oil and the driving force (pressure drawdown) [12] oil = .
To find PI the problem of non-stationary oil filtration from the reservoir to the fracture and through the fracture to the wellbore is solved. The scheme of this problem is presented in Fig. 3.
The process of non-stationary oil filtration in plane radial fracture with faces closed on the proppant at the moment p and then having the fracture opening distribution
p( ) = (, p) (, p), is described by the convection-diffusion partial differential equation p
p − * p (︂
︂) p p + * p 2 res = 0 and the equation of Darcy’s law p = − p p
.
* = p( oil + p), In (19), (20) p(, ) is pressure in the fracture; p is fracture permeability; p is Darcy flux in the fracture; res is Darcy flux of the oil from the reservoir to the fracture; total system compressibility * is defined as where p is the porosity of the proppant in the fracture, oil and p are the coefficients of oil and proppant compressibility correspondingly. It is supposed that oil is a newtonian fluid with the coefficient of dynamic viscosity .
The following boundary conditions are set on the wellbore (17) (18) (19) (20) (21) (22) (23) (24) and at the fracture front Here oil( ) is debit of the wellbore. The initial condition for (19) is established
The process of non-stationary oil filtration from the reservoir to the fracture is described by the convection-diffusion equation [13], [14] and at the fracture the following condition is set
Initial data for (25) is established res(, , 0) = 0, 0 ≤ ≤ res. and the equation of Darcy’s law In (25), (26) res(, , ) is reservoir pressure on the cylindrical surface of radius (see. Fig. 3); res is permeability of the reservoir; res is Darcy flux of the oil along the cylinder surface; total compressibility ** is defined by the formula ** = res( oil + res), where res is reservoir porosity, res is reservoir compressibility.
At the distance res from the fracture the condition is set res
res 2 res = 0 − ** 2
The stated two subproblems of the production model are solved iteratively by means of interchange of parameters res(, 0, ) and p(, ). Differential equations (19) and (25) are solved numerically using the finite difference methods. 3
Inverse problem of hydraulic fracture propagation
The vector of free design variables x (which have been independently varied in this study to find an optimum design) includes injection rate of fracturing fluid in( ), injection proppant concentration (volume fraction) in( ), fracturing fluid parameters , , 0, the Young’s modulus and the Poisson’s ratio of the rock, the Carter’s leak-off coefficient . In that way x = ( in( ), in( ), , , 0, , , ).
(31)
The vector of free variables y includes debit of the wellbore oil( ), oil viscosity , fracture p and reservoir res permeability, fracture p and reservoir res porosity, compressibility coefficients of proppant p, oil oil and reservoir res. In that way y = ( oil( ), , p, res, p, res, p, oil, res). (32)
lows:
The design variables are constrained within lower and upper bounds as fol l ≤ ≤ u , = 1, . . . , .
To build the functional (36) one should solve all the problems formulated above: fracture propagation during the time caused by the pumping of the fluid-proppant mixture; non-stationary oil filtration from the reservoir to the fracture closed on the proppant and through the fracture to the wellbore during the period of time 0 ≤ ≤ prod. It is quite complicated even if the considered 1D problem statement is used. Therefore some simplifications are implemented. E.g., in [1], [2] instead of solving the filtration problem numerically its approximate closed-form solutions are used for transient period [13], [14] and pseudo-steadystate [15], derived with the assumption of the constant pressure in the wellbore and using the empirical relations between flow rate and pressure gradient in the near-wellbore zone. In practice even simpler formulation is often used. E.g., according to [12] in terms of productivity the optimal fracture closed on the proppant can be considered the one for which the following relation is fulfilled = 1.6 res .
p
Design constraints are formulated to prevent uncontrolled fracture growth, multiple secondary fracture initiation, excessive fluid loss, to ensure that the designed treatment program can be executed in the field using specified surface equipment and downhole tubing, to ensure adequate fracture width, fluid efficiency and desired geometric proportions. Design constraints are stated as the inequalities: For example, to constrain below the minimal fluid discharge to the fracture in (34) the following function is assigned
The objective of the hydraulic fracturing is the increase of reservoir productivity. Therefore in optimization problem it is reasonable to maximize the cumulative production over time prod years by means of minimization of the major design objective (x) ≤ 0, = 1, . . . , . 1(x) = imnin
− m in in( ). 1(x, y) = −
oil( ). prod ︁∫ 0 (33) (34) (35) (36) (37) Then the maximization of cumulative production during the time prod can be carried out by means of the minimization of the functional 2(x) = ⃒⃒⃒ p(0) ⃒ frac( ) − 1.6 res ⃒⃒ ,
⃒ p ⃒ for which one need to solve just the problem of fracture propagation during the time
caused by the pumping of the fluid-proppant mixture.
As far as treatment cost of hydraulic fracturing is quite high, it also should be considered during design work. One of the ways of reduction of the treatment cost is the increase of fluid efficiency, i.e. the relation between the final fracture volume and the volume of the fluid pumped into the fracture. The maximization of fluid efficiency corresponds to minimization of the fluid leakoff total volume that is achieved by using the functional 3(x) = (,
For low-permeable reservoirs with low value of res in (37), such as shale, it is urgent to create the fractures of large radius. The problem of radius maximization can be reformulated as the problem of minimization of fracture width after the pumping stops
4(x) = ( w, ), and the problem of average fracture front velocity maximization can be reformulated as the problem of averaged fracture width opening minimization. 1 ∫︁ 0 5(x) =
( w, ).
To demonstrate the possibilities developed methods for multiobjective design optimization of hydraulic fracturing in this study these functionals are considered.
It is necessary to find the parameter vector (38) (39) (40) (41) (42) (43) providing the minimums for the functionals
x = ( 1, . . . , ) ∈ X, x∈X min 1(x), . . . , min (x) x∈X in the presence of bound (33) and design (34) constraints. In (43) the numeration of the functionals is not linked with the numeration of the particular functionals defined above. It is supposed to consider ≥ 1 abstract functionals. 3.2
Since in general case it is impossible to find the vector x minimizing two or more functionals at the same time, the solution of the problem is Pareto front. In [1], [2] building Pareto front is carried out by solving the series of one-objective optimization problems with freezing one of the functionals by means of the design constraint. The technique used in our paper allows building Pareto front directly. In case of multi-objective optimization Pareto front allows choosing compromise between several performance criteria. Optimization problem (42), (43), (33), (34) was solved by Genetic Algorithm that was used earlier by the authors for multi-objective shape optimization of hydraulic turbines [16].
0.35 0.3 0.25 , 0.05 0 0 2 1 150 t, s 50 100 200 250 300 (45) with i0n = 0.022 m3 / s and in = 0.350 m3 / s, obtained from the solution For fracturing fluid discharge law variation during the time interval from 0 to we used its representation in the form of second order polynomial in( ) = 2 + + , 0 ≤ ≤ , which coefficients are found from the conditions in(0) = i0n, in( ) = in, in ( ) = i*n and have the following representation = (︀ −6 i*n + 3 in + 3 i0n)︀ / 2, = (︀ 6 i*n − 2 in − 4 i0n)︀ /, = i0n. ︁∫ 0 (44) (45) (46) and in
The parameter i*n in (45) defines the volume of the pumped hydraulic fracturing fluid at constant injection rate in( ) and is set as non-varied parameter. Parameters i0n and in are the first components of the vector x: 1 and 2 correspondingly. In Fig. 4 two fracturing fluid discharge laws at = 300 s are presented: in( ) = i*n with i*n = 0.1 m3 / s and (44), (45) with i0n = 0.022 = 0.35, obtained from the solution of the problem 4.2 with bound constraints 0.02 ≤ 1 ≤ 0.1, 0.1 ≤ 1 ≤ 0.5. providing the minimum for the functional constraints. The values of the rest parameters are shown in Tab. 1. x = ( 1, 2) = ( i0n, in), x∈X The solution of the problem is the vector providing minimum 3 equal to 9.5 m3 and corresponding to fracturing fluid discharge law shown in Fig. 4 under number 2. If one uses four parameters x = ( 1, 2, 3, 4) = ( i0n, in, , ) instead of two (47) with bound constraints then the better minimal value of 3 equal to 8.4 m3 for the solution vector is obtained. It should be noted that the fracturing fluid discharge law is the same for solutions (51) and 54. Also there is insignificant difference between the dependences frac( ) and (, ) in the solutions of two- and four-parameter problems shown in Fig. 5.
If one minimizes 4 instead of 3 in (48) min 4(x) x∈X with the same constraints then two-parameter problem statement gives x = (0.1, 0.1) with minimal value of 4 equal to 8.4 mm, and four-parameter statement gives x = (0.1, 0.1, 0.5, 0.8) with minimal value of 4 equal to 5.4 mm. The dependencies of the direct problem solution for this vectors x are shown in Fig. 6. Finally, the solutions of minimization problems for 5 are the vectors x = (0.022, 0.35) and x = (0.02, 0.348, 0.53, 0.81) in cases of two (47) and four (52) parameter optimization. In the first case the minimum is 6.1 mm, in the second case – 3.9 mm. Corresponding dependencies of the direct problem solution for this x are presented in Fig. 7. 4.2
Let us find the vector (47) providing the minimum for the functionals min 3(x), x∈X min 4(x) x∈X with bound (49) and design (50) constraints. The solution of this problem is Pareto front presented in Fig. 8 with extreme points 1 and 2 marked. This points correspond to the vectors x1 = (0.1, 0.1) and x2 = (0.022, 0.35). In Fig. 9 the solutions of direct problem for this vectors are shown.
Using four parameters (52) instead of two (47) with bound constraints (53) gives Pareto front presented in Fig. 10. Points 1 and 2 on the front correspond to (52) (53) (54) (55) (56) Fig. 5: Minimization of 3: two-parameter (curves 2) and four-parameter (curves 4) 150 t, s (b) 150 t, s (a) 50 100 200 250 300 50 100 200 250 300 50 100 1t,5s0 200 250 300 5 10 15 r, m20 25 30 35 (c) (d) Fig. 6: Minimization of 4: two-parameter (curves 2) and four-parameter (curves 4). 2 4 2 4 150 t, s (b) 4 2 50 100 200 250 300 50 100 200 250 300 150 t, s (a) 0.35 0.3 50 100 200 250 300 50 100 200 250
300 2 1 2 1 35 30 25 R,mfrca2105 10 5 00 14 12 10
20 r, m 1 2 50 100 200 250 300 5 10 15 25 30 35
14 2 12 m 10 m ,4 F 8 6 8
14
F3 , m3 10 12 16 18 1 Fig. 10: Pareto front in four-parameter optimization problem 4.2 with marked points for analysis. (b) 2 00 1 2 1 Fig. 12: Pareto fronts in two-parameter (solid) and four-parameter (dashed) twoobjective optimization problems with extreme points 1, 1’, 2, 2’ correspondingly. the solution vectors x1 = (0.1, 0.1, 0.5, 0.8) and x2 = (0.022, 0.36, 2, 1). In Fig. 11 the solutions of direct problem obtained with this two vectors are shown.
In Fig. 12 Pareto fronts obtained from the solutions of two- and fourparameter two-objective optimization problems are compared. In tab. 2 the summary data for two- and four-parameter optimization problems is brought together. The parameter vectors x and minimal values for the functionals from the solutions of one-objective optimization problems and from the extreme points of Pareto fronts are compared.
Note that the extreme points 1 and 2 on Pareto front of the absolute minimum for functionals 4 and 3 correspondingly give the values of this functionals close to the ones obtained in one-objective optimization problem in which only one functional is minimized while the other is not taken into account. Hence the solutions obtained in 4.1 are the particular cases of the solution of two-objective optimization problem or the extreme points of Pareto front. min 3(x)
x
building of Pareto front.
The methods for optimal control of hydraulic fracture are proposed. The procedure consists of the following parts.
– The simulation of hydraulic fracture propagation caused by the pumping of
Herschel-Bulkley fluid and proppant slurry. – The computation of productivity of the fracture filled with proppant based on the combination of the models of plane-parallel and plane-radial filtration for the simultaneous description of fluid filtration in the reservoir and the – Genetic algorithm for solution of multi-objective optimization problem and
The method capabilities are demonstrated on the proposed functionals that have practical value and are used in non-automatic hydraulic fracturing design. By means of choosing the discharge law of the fluid pumping and its rheological parameters the following problems have been solved.
– Leakoff minimization. – Fracture width on the wellbore minimization (or maximization of its radius). – Time-averaged fracture width minimization (or maximization of fracture propagation velocity).
It has been shown that while solving multi-objective optimization problem the obtained Pareto front includes the solutions of one-objective optimization problems as particular cases.
It has been shown that if only fluid discharge law is varied the it is impossible to reduce the fracture width and fluid leakoff simultaneously. But if the variation of fluid rheological parameters is allowed then it is possible to decrease fracture width and leakoff volume simultaneously opposing to the case of fixed fluid rheology.
Authors gratefully acknowledge the financial support of this research by the Russian Scientific Fund under Grant number 14-11-00234.