The fluid flow connecting the hydraulic fracture and associated unconventional gas or oil reservoir is of great importance to explore such unconventional resource. The deformation of unconventional reservoir caused by heat transport and pore pressure fluctuation may change the stress field of surrounding layer. In this paper, the stress distribution around a penny-shaped reservoir, whose shape is more versatile to cover a wide variety of special case, is investigated via the numerical equivalent inclusion method. Fluid production or hydraulic injection in a subsurface resource caused by the change of pore pressure and temperature within the reservoir may be simulated with the help of the Eshelby inclusion model. By employing the approach of classical eigenstrain, a computational scheme for solving the disturbance produced by the thermally and pressure induced unconventional reservoir is coded to study the effect of Biot coefficient and some other important factors. Moreover, thermo-poro transformation strain and arbitrarily orientated reservoir existing within the surrounding layer are also considered.

The International Energy Agency projections indicate that an increasing supply gap for liquid fuels will arise if discovery and investment in new fields are unable to meet the urgent requirements with the rapidly increasing consumption of resources. Although the unconventional sources will lower energy costs and allow the substitution of natural gas for coal for electrical power production, the lack of equipment and experience for exploring the complicated geological structure may inhibit or slow down development. The prospect of higher prices accompanying a supply crunch may lead to an acceleration of investment in the unconventional resources, such as shale gas and other frontier oils [13]. Hydraulic fracturing, whose process involves the high-pressure injection of fracking fluid into a wellbore to create cracks in the deep-rock formation, is a well stimulation technique commonly utilized in unconventional reservoirs.

Solving the thermally and pressure induced reservoir is of great importance to petroleum industry and the associated applications [46]. The fluid flow connecting the fracture and surrounding rocks, heat transport, and temperature changes within the resource may significantly disturb the stress fields and deformation of the unconventional oil reservoirs including the shale and tight sand reservoir. The well-known Eshelby solution has been widely utilized to determine the effects of inhomogeneities which differ from the regional material [7, 8]. Eshelby’s technique can be also applied to deal with the elastic fields of stress, strain, displacement, deformation gradient around the intrusions, reservoirs, and underground structures due to the pore pressure or temperature changes [9, 10]. Zhou et al. investigate the background, applications, and future research directions of inclusion problem, as well as the classical methods including the analytical and numerical approach to treat associated problems [11].

Extensively studies on the geostructures and applications of hydraulic fracture have been explored over the decades. An analytical solution considering the porosity and permeability of gas storage is presented and may be applied to study the corresponding properties of shale gas reservoirs with noncircular nanopores [12]. Cai et al. explore the mechanism and effect of shale bedding of supercritical carbon dioxide hydro-jet fracturing using the experimental method [1315]. The problems of two interacting ellipsoidal inclusions and inhomogeneities with arbitrary orientation are solved to determine the elastic stress and strain components using Eshelby’s technique, and such method can be extended to discuss the pressure induced fracture [16, 17]. The failure mechanism, crack initiation, and propagation of brittle granite are evaluated on the basis of discrete element method [18], as well as the thermally induced microcracking [19]. The interpolations of the experimental data are proposed to obtain the characterization of fractured medium and corresponding crack porosity according to Eshelby-Cheng model [20]. Kanin et al. adopt the Carter’s leak-off model to assess the effect of pore fluid leak-in, near-tip cavity, and the penny-shaped crack [21]. Mauludin et al. investigate the effect of mismatch fracture properties and simulate the circular capsule with different types surrounding by the mortar media utilizing the cohesive-zone model [22]. The hydrodynamic flow model is proposed to investigate the properties and distributions of heterogeneous fluid, the interaction between the liquid and wall, and liquid–liquid slip [23]. Eghbalian et al. use the approach of representative elementary volume to estimate the failure mechanism of brittle rocks containing the microcracks and nanopores [24]. Pang et al. study the mechanism of kaolinite swelling produced by CO2 and water with the help of molecular simulation method [25]. Cai et al. develop a numerical method to evaluate hydraulic fracture propagation with the help of the flow-coupled cohesive interface coupled method [26].

The problem of deformation and stress disturbance around a reservoir due to the heat transport and fluid migration may be solved on the basis of the Eshelby model, which involves the concept of eigenstrain to handle the inelastic issue [2729]. In this study, the stress fields around a penny-shaped reservoir with thermo-porous eigenstrains in a full space are investigated via the numerical equivalent inclusion method. Section 2 introduces the thermal and porous eigenstrains and the classical equivalent inclusion method, as well as the corresponding Eshelby tensor. The computational scheme to treat the arbitrarily orientated reservoir surrounding by the diffident soil and rock materials is proposed in Section 2. Parametric studies and examples on the Biot coefficient, distribution of a reservoir, effect of linear thermal expansion, temperature changes, pore pressure fluctuation, and elastic stress disturbance are provided. The applications of proposed solution and computational scheme are given in Section 3. Conclusions are presented in Section 4.

2.1. Thermo-Poro Transformation Strain

The formulation of poroelasticity is a branch of the classical elastic theory and can be extended to study the thermal problem for a continuous porous material. Mura [30] extends the celebrated work by Eshelby [7, 8] and uses the term “eigenstrain” to represent the inelastic strain, including the thermal expansion, pore pressure change, misfit, and dislocation. The distribution of uniform or nonuniform eigenstrains in a material will lead to the corresponding eigenstresses and also affect the mechanical properties of the surroundings. The well-known eigenstrain method can be widely utilized to handle the void, crack, and other material defects as strain released during geometry changes. Therefore, changes due to heat expansion or cold contraction in a porous medium may be computed with choosing the proper eigenstrain via the Eshelby model [31]. The eigenstrain caused by the pore pressure fluctuation and temperature change may be expressed as

The terms μ and ν refer to the shear modulus and Poisson’s ratio, Kronecker delta δij is 1 if the variables i,j are equal, and 0 otherwise. The Biot parameter α=1Cs/C, and the compressibility of the solid phase or called the unjacketed coefficient of compressibility is denoted as Cs, which is usually considered identical to the compressibility of the grains [32]. The corresponding jacketed compressibility coefficient C means the drained compressibility of the porous medium. Furthermore, the constitutive equation of thermo-poro eigenstrain (1) also depends on the unstrained state of pore pressure P0 and temperature T0 and the strained state of pressure P and temperature T. The term λ represents the coefficient of linear thermal expansion, which characterizes the ability of an object to expand under the effect of temperature elevation. It is worth noting that if the porous reservoir is compliant, which is compared to the solid constituents, α=1Cs/C~1, whereas if the resource is stiff α=1Cs/C~0.

2.2. Formulation

Eshelby’s work [7, 8] has shaped the fields of micromechanics and provided some outstanding and fundamental solutions to address the problems in solid with defects, especially for the fracture mechanics, plasticity and engineered composite materials, etc. It is noted that Eshelby’s equivalent inclusion method in elasticity is one of the most famous and well-known micromechanical approach to analyze the deformation of arbitrarily shaped inhomogeneity, which is restrained by its surroundings with different elastic constants. The general ellipsoidally shaped inhomogeneity is versatile enough to cover a wide variety of special shape including the penny-shaped reservoir with thermal and porous transformation strain. In this case, an isotropic, linear, and full-space medium applied with remote external load σij0 contains an ellipsoidal inhomogeneity, whose elastic constants Cijkl differ from those of the rest of the surrounding matrix Cijkl. In this case, when letting three axes of the ellipsoid a1a2a3, the general shape of ellipsoid transforms into the penny-shaped reservoir, as shown in Figure 1.

The nonlinear inhomogeneity problem may be converted into linear inclusion problem via the equivalent inclusion method. The constitutive equation of equivalent inclusion method is
wherein the term εij0 denotes the remote strain field, which can be calculated in absence of the inhomogeneity according to Hook’s law, σij0=Cijklεkl0. In addition, εijP is the preliminary eigenstrain distributed inside the reservoir and is given in Equation (1). The equivalent eigenstrain εij may be expressed as
The fictitious eigenstrain εij can be solved according to the constitutive equation in Equation (2). Moreover, Equation (2) can be rewritten in terms of the classical Eshelby strain tensors Sijkl. The detailed formulation of Eshelby tensors is listed in [30].

Solving above equation can obtain the fictitious eigenstrain εij and equivalent eigenstrain εij. Equation (4) is enforced at a set of particularly chosen points, resulting in a system of algebraic equations, which can be numerically solved. Thus, the stress field of thermally and pressure induced unconventional reservoir may be determined by adopting Eshelby’s equivalent inclusion method, as shown in Figure 2.

When the penny-shaped reservoir and surrounding rock matrix are both isotropic materials, Equation (4) may be expressed as [30]
wherein the parameters λ and μ refer to the Lame constants of the penny-shaped reservoir and λ and μ denote the corresponding Lame constants for the rock matrix. In order to obtain the shear components of εij, the deviatoric strains are introduced.
By adopting the above expressions, Equation (5) arrives at
where bulk modulus, K, can be expressed in terms of Lame constants and Young’s modulus.
After the algebraic operation, the shear parts of εij may be calculated as
In this special case, the reservoir may be assumed be a penny-shaped inhomogeneity, whose geometric shape is a1=a2 and a3=0. Then, Equation (5) becomes
and the equivalent eigenstrain can be solved as

To make the case a bit more general, the following examples and discussions of the geometry for the resource inhomogeneity are set to be a1a2a3. Due to the arbitrarily orientated reservoir existing within the surrounding rock after hydraulic injection, it is necessary to apply the Euler angles to evaluate this problem. According to Euler’s rotation theorem, any rotation can be measured by utilizing three angles. Thus, in conjunction with the equivalent inclusion method, the stress distributions of thermally and pressure induced reservoir with arbitrary orientation can be numerically solved.

In order to calculate the stress distributions of thermally and pressure induced reservoir surrounding by the unconventional oil reservoir, a computational code is programmed to investigate some fundamental parameters affecting the elastic stress fields. In this section, dimensionless factors are employed to study the penny-shaped reservoir caused by thermo-porous eigenstrains. The practical parameters including the Young’s modulus, Poisson’s ratio, remote applied loads, coefficient of linear thermal expansion, Biot’s coefficient of the unconventional reservoir, and rock matrix are considered. Due to the randomness of reservoir, the Euler angles to represent its orientation are also involved to determine the Von Mises Stress of reservoir caused by the heat elevation or cold contraction and pressure raising or dropping before and after the hydraulic injection. To make the case a bit more general, the values of Biot’s coefficient are on the basis of Bakken shale samples from Williston Basin; the text measurements of samples are detailed in [33]. Thus, the practical factors of Young’s modulus, Poisson’s ratio, pore pressure change, temperature fluctuation, coefficients of linear thermal expansion, remote applied loads, and rotation angles of reservoir can be input into the codes to evaluate this problem. For the purpose of investigating each factor affecting the final results, the general parameters are set as follows. Reservoir is embedded in an infinitely extended shale matrix with remote loads applied σ110=σ220=σ330=1 at infinity, the size of reservoir is a1=1,a2=0.75,a3=0.01, Young’s modulus and Poisson’s ratio for the reservoir matrix and surrounding rock E=30000Mpa and ν=0.3, the corresponding coefficient of linear thermal expansion λ=3×105, Biot’s coefficient for the unconventional reservoir α=1Cs/C=0.7, pore pressure change ΔP=15, and temperature fluctuation ΔT=50. The target line is from -2 to 2, and the Von Mises stress distribution and the each stress components along this line is obtained.

The effect of parameter on the Young’s modulus for the surrounding rock matrix with respect to those of reservoir Ereservoir=30000Mpa may be computed according to the formulation detailed in Section 2. The calculated Von Mises stress for a penny-shaped reservoir starting from x1/a1=2 to 2 is shown in Figure 3(a). Results are displayed as the ratio of Ematrix/Ereservoir; the Von Mises stress reaches the maximum points at x1/a1=1,1 when Ematrix=5/3Ereservoir and minimum points at x1/a11,or1. That means the jump between the matrix and reservoir has its maximum value when Ematrix is larger than Ereservoir. The effect of Poisson’s ratio on the Von Mises stress is plotted in Figure 3(b). The Poisson’s ratio for the reservoir is set to be 0.3; the corresponding Poisson’s ratio for the matrix varies from 0.2, 0.25, 0.3, 0.35 to 0.4.

Results are given for two parameters based on the equation of eigenstrain caused by the pore pressure fluctuation and temperature change, which is related to the coefficients of Biot and linear thermal expansion. Figures 4(a) and 4(b) show the normal stress components for Biot’s coefficient α=0.7 and the coefficient of linear thermal expansion λ=1×105,3×105, respectively. The difference between 4(a) and 4(b) is that the stresses σ11,σ22 when λ=1×105 show smaller than for the λ=3×105. It is noted that the stress σ33 exhibits almost the same tendency inside and outside the reservoir due to the target route lying on the x-axis. In reality, the stress distribution inside the inhomogeneity would like to be disturbed by surrounding geostructures such as fault and cracked layer. In this case, the reservoir is assumed to be perfectly bounded by the rock matrix, without being disturbed by any other geostructures or neighbor resources. The total tendency of Figures 4(a) and 4(b) is seen to be almost equal; only the maximum and minimum stress values differ.

As expected, before or after the hydraulic fracturing, there is a temperature difference between the injected fluid and surrounding matrix underground. Figure 5 displays the predicted normal stress distributions of the penny-shaped reservoir due to temperature changes produced by hydraulic injection. Figure 5(a) presents the normal stresses after injection leading to the temperature decreasing of reservoir matrix caused by the injected cold water. In like manner, Figure 5(b) describes the normal stress distributions when the temperature elevates. It can be seen that stress σ33 keeps a constant as the temperature heating or decreasing as opposite to the stresses σ11 and σ22. The results tell that the stress disturbance may be influenced by the injected cold fluid, and the associated stress field of reservoir would like to be strongly disturbed leading to the stress concentration. The stresses inside the reservoir show a big difference among each stress component if the shape of resource is irregular.

As we know, the pressure of reservoir formation would like to increase with the injection of fracturing fluid until it peaks at the breakdown pressure. In order to compare the pressure dropping and increasing, the last example describes the pore pressure fluctuation affecting the total stress field after or before the fracturing work. Figures 6(a) and 6(b), respectively, exhibit the predicted normal stress parts of the reservoir produced by the change of porous eigenstrains, which means the pressure dropping ΔP=PP0=30 and pressure elevating ΔP=PP0=30. In Figure 6(a), the stresses σ11 and σ22 reach their maximum value at the interface between the resource and formation when ΔP=30, while the stresses σ11 and σ22 when ΔP=30 in Figure 6(b) are lower than those of Figure 6(a). It is concluded that the pore pressure of reservoir matrix dropping or elevating may both influence the stress components σ11 and σ22 at this special case, but stress σ33 in Figures 6(a) and 6(b) shows no change due to any pore pressure fluctuation along the target line from x1/a1=2 to 2 because of the special shape of reservoir.

Fluid production or hydraulic injection in a subsurface resource would like to change the pore pressure and temperature within the reservoir, which generally induces the alteration of stress distribution within and around the reservoir. The general ellipsoid is more versatile to cover a wide variety of special case, including the penny-shaped inclusion. Using the classical Eshelby inclusion method to solve the induced stress changes within an ellipsoidal resource in a full space or half space is an effective way to deal with the geostructure problem. In this paper, a computational scheme for treating the stress distributions a penny-shaped reservoir with thermo-porous eigenstrains in an infinitely extended formation is explored with the help of numerical equivalent inclusion method. Parametric studies on the Biot coefficient, effect of linear thermal expansion, temperature changes, pore pressure fluctuation, Young’s modulus, and Poisson’s ratio of the unconventional reservoir are performed. Results are displayed as the ratio of Ematrix/Ereservoir, the Von Mises stress reaches the maximum points at x1/a1=1,1 when Ematrix=5/3Ereservoir, and minimum points at x1/a11,or1. That means the jump between the matrix and reservoir has its maximum value when Ematrix is larger than Ereservoir. Some other results and discussions are detailed in Section 3.

The processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

The author(s) declare(s) that there is no conflict of interest regarding the publication of this paper.

This work is financially supported by the college programs of Chongqing Industry Polytechnic College (Nos. GZY2020-GGRC-39 and GZY2021-ZX-10). The authors would also like to acknowledge the financial support from research funds for the doctoral program of Chongqing Industry Polytechnic College.

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).