Mathematical Problems in Engineering

Volume 2016 (2016), Article ID 4052483, 12 pages

http://dx.doi.org/10.1155/2016/4052483

## Finite Element Optimised Back Analysis of In Situ Stress Field and Stability Analysis of Shaft Wall in the Underground Gas Storage

^{1}College of Mechanical and Electronic Engineering, China University of Petroleum, Qingdao 266580, China^{2}College of Storage and Transportation and Architectural Engineering, China University of Petroleum, Qingdao 266580, China

Received 14 August 2015; Revised 11 January 2016; Accepted 13 January 2016

Academic Editor: Masoud Hajarian

Copyright © 2016 Yifei Yan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

A novel optimised back analysis method is proposed in this paper. The in situ stress field of an underground gas storage (UGS) reservoir in a Turkey salt cavern is analysed by the basic theory of elastic mechanics. A finite element method is implemented to optimise and approximate the objective function by systematically adjusting boundary loads. Optimising calculation is performed based on a novel method to reduce the error between measurement and calculation as much as possible. Compared with common back analysis methods such as regression method, the method proposed can further improve the calculation precision. By constructing a large circular geometric model, the effect of stress concentration is eliminated and a minimum difference between computed and measured stress can be guaranteed in the rectangular objective region. The efficiency of the proposed method is investigated and confirmed by its capability on restoring in situ stress field, which agrees well with experimental results. The characteristics of stress distribution of chosen UGS wells are obtained based on the back analysis results and by applying the corresponding fracture criterion, the shaft walls are proven safe.

#### 1. Introduction

In situ stress is an important factor that should be considered during the construction of underground gas storage (UGS). A direct way to get the in situ stress is to measure the principal stress by test methods, such as micro-fracturing test [1].

As the engineering technology develops, kinds of techniques have been selectively accepted for obtaining accurate in situ stress orientation and magnitude. Hydraulic fracturing, when poroelastic effects are considered, seems to be the best option to derive the horizontal in situ stress magnitude and direction [2]. If a vertical fracture is created in the process, then the minimum in situ stress must be in the horizontal direction. Although fracturing has been the most commonly used method in in situ stress measuring, it may not be available or even fail sometimes, for example, in inclined test holes where the fracture direction is not normal to the minimum horizontal stress. On the other hand, the fracturing method will also be ineffective when the measurement requires a higher injection pressure and the terrain, with UGS located, is bearing a very large stress. Some researchers implement methods such as leak-off tests and breakout orientation to the standard well testing and logging; however, these methods are not always reliable. Moreover, the stresses obtained from breakout analysis are confusing as breakout can be induced by different kinds of in situ stresses.

When in situ stress data from test wells is obtained, some techniques can be used to identify the magnitude and direction of the whole in situ stress field [3]. These techniques are organised as an overdetermined system of equations, and the in situ stress state can be solved by using several methods mentioned in some articles [4–6]. Not only can the magnitude of the two principal stresses be solved, but also the direction of the field can be determined. However, the methods assume that the shear stress is negligible as it is zero when the normal stress is a principal stress. Djurhuus and Aadnøy extended the study on in situ stress and presented a new theory [4]. They obtained experiment data from multiple fracturing and induced fractures from imaging logs to calculate the in situ stress state. The magnitude information is mainly calculated from fracturing data, while the orientation information can be obtained from image logs.

However, measurement is difficult to carry out because of high cost on an injection-production well. Therefore, several in situ stress values can only be obtained through experiments. The stress distribution of the test zone can be calculated quickly with the help of finite element optimisation and the results can be more accurate with some optimisation methods implemented. To carry out the inversion calculation, we need to construct the finite element model with proper boundary conditions and loads. Some minor modifications should also be incorporated into the model to eliminate the effect of stress concentration, which will be discussed below.

Boundary conditions and loading methods cannot be input directly because the tectonic movement is unknown and the geologic structure is complicated. These circumstances make in situ stress field analysis difficult to carry out. So boundary load inversion comes to be the crucial way to conduct research on the in situ stress field. From the viewpoint of geomechanics, the tectonic stress field is the in situ stress field forming the tectonic system and type, including its distribution area and inside stress distribution. During the analysis, rock density could be determined relatively accurately in a small variation range with different precision levels of gravity stress and tectonic stress. Instead of being inversed, based on terrain conditions, the finite element method could be used to calculate gravity stress to ensure the precision can satisfy the engineering requirement. However, tectonic stress field calculation is different [7]. During the construction of UGS, the distribution of the in situ stress field along the plane at the top of the salt cavern is needful while the in situ stress information of only a fraction of wells in a certain area is available. Then, inversion and extrapolation need to be performed to obtain the stress field of the whole area. Among the present numerical simulations and analysis methods of the in situ stress field, the most widely used approach is cut-and-try and adjustment of boundary loads [8–13]. Given the complexity of the geological structure and the limited reference of the in situ stress, boundary conditions and loading methods become the difficult problem in calculation and analysis part of a numerical simulation. The major issue is how to improve the precision of the inversion and reliability of the data. Simulation is an effective way to take terrain and geological features into consideration. And only by this can relevant theoretical or numerical analyses such as back analysis, back calculation, and simulation be performed to predict the tectonic stress and its mode based on the measured data [12]. The boundary loads identification and optimisation method presented by this paper adopts damped least squares algorithm to approximate the loads by setting parameters from stress testing references. The initial value for iteration is got by the unit load method to avoid the limits from the structure geometric features and boundary conditions. By controlling the error of observing point, the identification on uncertain and unknown boundary loads of the in situ stress field is achieved. A few ways are attempted and compared to improve the data reliability and the precision of in situ stress field inversion calculation in this paper.

#### 2. Building the Finite Element Model

It is difficult to achieve an analytical solution of boundary loads for a formation due to the complexity of the geological structure and boundary conditions. The finite element method is employed here for mechanical analysis to make sure that the research method can have a general application.

Though salt is viscoelastic, a UGS reservoir stress field can be treated as an elastic problem to simplify the analysis. This study focuses on the original rock salt, without considering the effect of creep and relax. The linear elastic assumption may bring some deviation with the actual situation, but it is high needed for figuring out a specific engineering problem.

The applied loads generate displacement field, stress field, and strain field inside for a certain area. The basic equations of small deformation elastic mechanics are as follows: equilibrium differential equation: geometric equation: physical equation: where is the differential operator and is the elastic matrix depending on the elastic modulus of the formation and Poisson’s ratio .

Generally, the finite element formula iswhere is the global stiffness matrix, is the node displacement matrix, is the boundary load array, and is the parameter vector used to describe the boundary loads of UGS reservoir.

Equation (4) is transformed to

Equation (2) is substituted in (3) to get the following equation:

Then, (5) is substituted in (6); that is,

From (7), the nodal stress and the boundary loads satisfy the following equation:where . is the transfer matrix between the nodal stress matrix and the boundary loads array. For the linear elastic body, could be determined when the finite element model of UGS reservoir is established. is uncertain because the loads are unknown, and only a part of is certain, since the stress values of only a few wells in the reservoir to be studied are known. Therefore, in general, the boundary loads parameter vector cannot be calculated from (8) and is waiting to be discovered.

#### 3. Optimisation Method

In situ stress field inversion is used to find the model or parameters that can fit actual data under certain principles. It is not only a method of data fitting, but also a way of optimising. The inversion result is closely related to the optimal criterions, which is used to evaluate the degree of data fitting. Sometimes, we have prior information about some areas that have already been researched. Based on this prior information, we could restrain the optimal criterion and parameters. Generally, the data going to be used in stress inversion calculation is obtained from the reservoir stress field and contains the stress of logging, fracturing, borehole breakout, or core testing of certain wells. The principle of stress inversion (or objective function) could be described as the following models by using the least squares method as the optimal criterion.

*(1) Constraint Model*. If the stress value of the test well is , that is, the simulation value of the stress field from the finite element calculation is known, the model could be described aswhere is derived from the finite element calculation, is the mean squared error between the measurement and calculation, is the unknown boundary load vector, and and are the boundary load parameter functions of the upper bound and lower bound, respectively, determined by prior information. and are the boundary load parameters of upper and lower bound, respectively.

*(2) Unconstraint Model*. An unconstraint model exists when the condition in (9) is not met; that is,where symbols are the same as those in (9).

#### 4. Inversion Algorithm

According to the principle of inversion calculation, back analysis is used to calculate the input data of a system by analysing its output data [13]. Therefore, if the system is well constructed, it is possible to obtain the input data from its output data, which is usually achieved through experiments. The back analysis ensures that a unique solution can be generated, as long as the mechanical model is real and the data is sufficient.

The general way of inversion is to present the parameter vector of boundary load and then calculate the reservoir’s stress response value under . During this process, the finite element method is used to approach the actual data of some well sites based on the least squares method. Usually, the Gauss-Newton method is adopted to solve the least squares problem of (9) or (10). However, ill-conditioned coefficient matrix of the equation, low iterative speed, and spurious convergence may occur sometimes. In this paper, a damped least squares method [Levenberg-Marquardt (LM) method] [14] is implemented to approach the parameter vector of boundary load . The LM method, which has a high optimisation speed and precision, possesses both the advantages of the gradient method and the Newton method. When the initial value of damping factor is small, the step size is set to be equal to that of the Newton method, while it is equal to that of the gradient descent when is large.

When the iteration reaches step , the column vector of stress response and boundary load and are set. Then, series expansion around is performed and the trace above the second-order polynomial is ignored to get the following equation [15]:where is the Jacobi matrix of in , and it is defined from the finite difference method aswhere could be calculated from (7).

The search direction used in LM method is from the solution of a set of linear equalities. The iterative formula iswhere is the search direction vector of and is the unit matrix. Scalar is the damping factor controlling the direction and value of . When is 0, the direction of is the same as the result of the Gauss-Newton method. When tends to infinity, we could derive the following equation from (13):Thus,where is in the direction of the zero vector and has a steep descending slope for , which is large enough; that is,From (13), we derivedThen we got

The boundary load vector could be derived from iterative solving of (13)–(18). The key step is how to choose iterative initial value. If the range of the boundary load is known, we can set the iterative initial value of as the average of the upper and lower bounds. If the range is an unbounded domain, the iterative initial value of could be determined after experiential trials.

From the above discussion, if the measured data of reservoir stress is available, the boundary load could be determined. The load parameter vector and reservoir’s horizontal stress could satisfywhere denotes the transfer matrix, is the point number for the measurement, and is the dimension of the boundary load vector. If , a set of boundary load vectors could be clear based on the principle of controlling the comprehensive error to be minimum for the nonuniqueness of the solution.

The specific steps of finite element inversion of the in situ stress field go as follows: Firstly, choose a proper objective function. Secondly, adjust and search parameters by the optimal method and forward-model the computed value by the finite element numerical algorithm. Thirdly, substitute the measured and calculated values into the objective function and judge the result from the objective function. If the result did not reach the minimum, the adjusting and searching should be continued. Repeat this process until the result from objective function reaches the minimum value and finally, the boundary load parameter is generated. Afterwards, the forward modelling is carried out to acquire the reservoir’s stress field by the boundary load derived from the inversion. (Figure 1 shows the inversion process.) The applied boundary condition is shown in Figure 2. The calculation is performed in the entire circle area to avoid stress concentration. Assuming the strata are boundless and stress is continuously changed, if the stress applied on the circle can make the inversion results close enough to the test results, the stress applied on the rectangular must be right. The stress distribution in the object area is the result we want from the inversion.