Creep deformation of fracture surfaces analysis in a hydraulically fractured reservoir using the finite element method

In this study, we develop a Galerkin-type finite element (FE) solution to assess the creep deformation of fracture surfaces. These surfaces are treated as comparatively compliant rocks which are separated by a comparatively stiff proppant pack. Viscoelastic properties of fracture surfaces are assessed by the temperature and deformation dependent modulus. The key challenge is to appropriately develop the creep and the relaxation functions which are continuous spectrum material functions and are unique for each model. Time dependent mechanical properties are connected through an instantaneous elastic response and coupled by way of nonlinear and convolution integral relationships. Coupling is achieved by way of convolution integrals computed using the full-time history of the material functions. As a novel numerical procedure, we apply a more computationally efficient algorithm to use discretespectrum approximation to the creep and relaxation functions, thus, allowing for adjustment of the energy dissipativity of the model in a consistently low computation-time. The proposed approach is novel and is the only approach in the literature (to the best of authors’ knowledge) that considers the conductivity in an intended reservoir at the level of proppants. It will be demonstrated that the proposed model can be used for the assessment of fracture conductivity under short or long-term loading conditions when rock material properties can be determined.


INTRODUCTION
The economic development of unconventional petroleum reservoirs is frequently dependent on the creation of an effective hydraulically induced fracture system that can carry hydrocarbons to the wellbore.The petroleum industry spends significant effort and cost on hydraulically fracturing these reservoirs and on placing proppant to keep the induced and reactivated fractures effectively propped open during production.Low viscosity waterbased fluids (eg.Slickwater) with relatively low proppant volumes are commonly used in this process.Although field evidence is limited it is generally acknowledged that many such fracture treatments result in significant fracture area that is minimally propped (single monolayer of proppant), partially propped or not propped at all and yet is in hydraulic communication with the wellbore.This concept is supported by evidence from mineback experiments, offset well behavior, microseismic and tiltmeter fracture mapping, decline curve analysis, mass balance approaches to pumping and flowback volumes, and laboratory tests (Warpinski et al., 2008;Cipolla et al., 2009;Morales et al., 2011).Under the right conditions this minimally propped region may have sufficient conductivity to play a significant role in the long term performance of many fracture-stimulated unconventional reservoirs (Palisch et al., 2007;Vincent, 2012).Many factors may contribute to the reduction of initial fracture conductivity and we will present some work related to the time dependent mechanical behavior of fractures that may be relevant at the scale of the lifecycle of a well.
As fluid pressures stabilize following a fracture treatment a component of the far-field stress acts perpendicular to the fracture surfaces.In a newly propped fracture the relative stiffness and shape of a proppant grain may result in some embedment of the grain in the surface of the rock fracture.Alramahi and Sundberg (2012) have published laboratory results that successfully demonstrated a correlation between the Young's Modulus of the rock perpendicular to an ideal fracture surface and the depth of proppant embedment in that fracture surface under given loading and boundary conditions.Our goal in this work was to contribute to a better characterization of the performance of propped fractures over the lifespan of a well by incorporating both short and long term effects on embedment and fracture conductivity.Several time dependent effects are potential factors in this process, including physical and chemical degradation of the proppant (McDaniel, 1986;Duenckel et al., 2011;LaFollette et al., 2011), time dependent chemical reactions between the rock and the fracture fluid (Horsrud et al., 1994;Akrad et al., 2011), and the time dependent mechanical behavior of the propped fracture surface (Akrad et al., 2011;Sone and Zoback, 2011;Yawei and Ghassemi, 2012).
In this paper we will present the results of a numerical model that can be used to simulate the time-dependent behavior of propped fractures.Some rocks from producing reservoirs exhibit significant time-dependent strain behavior which over the lifespan of a producing well would be expected to lead to a loss of fracture width due to time-dependent proppant embedment.This process is being observed in the controlled environment of a lab and while lab results are being obtained a numerical model was developed to simulate the observed results.Günter et al. (2010) added to the large body of work on time dependent behavior of salt rock by describing a constitutive model appropriate for describing the elastic and visco-plastic behavior of rock salt pillars via finite difference methods.Guo and Liu (2012) described an analytical method for applying elastic and visco-elastic equations to simulate proppant embedment by considering the process to involve an evenly distributed load along the fracture surface.Here we present a model that considers the dissipation of strain energy within the rock as stress is applied at the fracture surface resulting in a loss of fracture width through time under given loading conditions.The model uses a finite element solution that reduces the computation time required when simulating more complex geometries and scales of investigation.

METHOD
The methodology provided in this section is adapted for a viscoelastic/viscoplastic material model and is then transformed into a Galerkin-type finite element formulation.

The viscoelastic model
The model applied in this study relates the internal pressure, , and the area between the fracture surfaces, .Assuming the medium is purely elastic, this area is denoted as ., and are related using a convolute integral such that , in which and represent nonlinear relationship and convolution integral relationship, respectively.The convolution integral is a relationship defined between the applied pressure and the measured strain field: (1) is the zero pressure area, is the zero pressure fracture width, is the static modulus of elasticity of the medium, is the effective fracture wall thickness, and proppants are assumed rigid.and are, also, linked to and in the form of convolute integrals. (2) is the normalized creep function and elastic model is, in fact, the zero frequency limit of the viscoelastic model.Equation (2) can be written as: (3) is the normalized relaxation function.Both and are material functions which are continuous spectrum of which is defined as:  2) and (3) using full-time history of variables for the case of continuous-spectrum material functions.For the sake of ease, here we define discrete-spectrum functions for relaxation and creep behavior such that: (5) We note that for , and is the normalized static modulus of elasticity so as .Also, is the number of discrete time constants ( with their corresponding .At the onset of the simulation the fracture is assumed to be under nonzero and hydrostatic pressure (preload).For this case the immediate elastic response will take the form: is the area of reference pressure.The discrete form of the above equation for time takes the form: The variable saves the history of deformation for the j th component of and at each time-increment it is updated as such: (8) The immediate elastic response and the pressure maintain the relationship as such: (9) We note that if all are zero.

The finite element procedure
A discrete Galerkin-type finite element is applied in this study.Each element is considered a control volume in which the internal stress is defined based on the strains in the control volumes.This strain accounts for the flow motion (temporal large strains) of the corresponding control volume.The equilibrium equation in tensors is defined: (10) Applying the weighted residual formulation in the weak form yields: (11) Due to the divergence theory: (12) and the governing equation takes the form: (13) The isotropic viscoelastic material based on the small strain theory is defined: is the shear relaxation function and is the bulk relaxation function.
As mentioned, viscoelastic properties of fracture surfaces are assessed by the frequency dependent modulus, which is also termed the complex shear modulus.The relaxation function can be broken down into elastic and viscous terms in order to appropriately define the intended viscoelastic model.While developing the governing equation, to define the stiffness matrix, the stiffness takes the form .In the stiffness matrix formulation, the integration of and is handled in the same way such that: After the integration of the weak form the nodal solutions of stress and strain fields are available.

RESULTS AND DISCUSSION
We apply the formulation developed in this study to a few generic models with somewhat extreme material properties and conditions in order to verify the robustness of the formulation.Then, the transient and steady-state stress and strain for a viscoelastic material model of fracture surfaces is analyzed.The time-variables and spatial positions of each node are coupled to allow for any relaxation modals.For the steady-state viscoelastic analysis, the equilibrium equations are built by coupling the control volumes at the nodes.It should be noted that in this stage thedirect validation was not possible, however a tangible validation will be performed when the proposed model is further improved.
For the viscoelastic model, we assume that the

Example 1
We apply a torsional test on a cylindrical shell model at the frequency of 1 Hz and temperatures of -20 to 20°C with an temperature increments of 10°C for which a different torsional angle, , is applied (Figure 1).
The relationship of and is , r is the radius and H is the height of the model.To proceed with the solution, we consider 2 nd order polynomials to define the from that, that is, small and large deformations, tends to decrease to provide more stiffness to the viscoelastic model.
In this study we perform the developed formulation on a simulated fracture reservoir.To assess the impact of fracture deformation on conductivity two options have been proposed: (1) mechanically model the propped fracture to obtain the geometry and properties of the fracture surfaces and proppants.Perform a fluid/solid interactions analysis and apply appropriate pressures in a multiscale scheme to define the conductivity of the propped fracture by considering its physical dimensions.This can be done by applying Darcy's law or Forchheimer equations depending on the ratio of pressure drop to fluid velocity.(2) Fracture conductivity can be assessed indirectly by estimating the geometry of the fracture, that is, fracture width and the proppant size using fracture mechanics analysis.This analysis can be compared to measurements of flow conductivity under similar laboratory loading conditions that match the model case.The latter method was considered in this study in order to maintain control over the many variables at play in the

Contact mechanics
To model contact, we apply the method of interactive computation of a push back force in which contact surfaces and the corresponding force is determined (ref).To characterize the kinematics of the contact surfaces, that is, the proppant and fracture surfaces, two spherical and flat surfaces are defined and denoted as js 1 and js 2 in a 3D space.The distance of the two surfaces (node to node) is defined by a scalar function, g(x), which is denoted as the "gap" function.We assume a frictionless contact so that the contact force is normal to the contacting surfaces.This is defined as a partial derivative of g(x) with respect to the spatial coordinates as such (16):

Pressure follower force
Force vectors need to be defined such that their directions always stay normal to the surface in the current configuration each time (Mohammadi et al., 2011(Mohammadi et al., , 2013)).A face of a 3D element where the pressure is applied ( ) on the current configuration (unit normal n) is considered.The traction force vector (t) is expressed as n and the corresponding virtual external work in the current configuration such that: is the displacement vector of the surface, n is the normal vector to the surface, is the pressure and is the external work done by the external pressure.To assess the robustness of the proposed model, we develop a model in which one single spherical proppant is placed between the fracture surfaces.A 3D model in the context of large deformation theory with a linear (elastic) material model is applied.We allow 0.1 mm penetration of proppant into the fracture surface.The modulus of elasticity of the fracture surfaces and proppants is considered 20 and 40 GPa, respectively, and the radius of the proppant is 1 mm.The model is shown in Figure 4.
The mesh model (Figure 5) should be precise to ensure the stress field close to the contact is calculated irrespective to the element size and topology implemented in the study.
After performing the FE solution the stress distribution around the contact zone (Figures 6 and 7) is calculated.The calculation of the stress field around the contact surface is of particular importance because the failure phase and how and when it is initiated and propagated is the key indicator of what type of stress causes failure in mechanical interaction between the proppants and the fracture surfaces.

Failure analysis
After the areas with high stress are detected we can determine whether failure is initiated in those areas.We apply the following procedure to proceed with the failure analysis as shown in Figure 8.The idea is that the elements are deleted from the mesh model if the calculated von Mises stress in those elements exceeds the ultimate compressive stress (UCS).
The final geometry is based on the remaining elements after iterative model runs.
For dramatic demonstration purposes we consider a closure stress of 100 MPa and a UCS of 70 MPa for proppants and 50 MPa for the fracture surfaces.The failure analysis in 3 iterations is shown in Figure 9.
It is acknowledged that this algorithm is flawed in that the "deletion" of material is not what happens when rock material fails under proppant loading, as such the elimination of an element entirely is not strictly the right assumption and that approach will lead to overly pessimistic predictions of failure and embedment damage.
One way to address this issue is to maintain some portion of the volume of the failed elements and add it back to the new geometry (obtained from the previous iteration).This approach allows for the use of more realistic post failure geometries.As shown in Figure 10 this will improve the failure analysis, for instance, in  iteration #3 and the penetration of the proppant into the fracture surfaces reduces to 2.41 μm and the modeled post loading fracture width changes from 0.670 to 0.951 mm.It should be noted that the material properties used this demonstration were chosen for demonstration purposes and not to approximate the material properties of any particular rock-proppant-stress system.

Viscoelastic model
In the last step we apply the proposed viscoelastic model to an array of 3x3 proppant grains placed in between of the two fracture surfaces.Fracture surfaces are considered viscoelastic and the proppants are assumed to be rigid with infinitive stiffness.The FE mesh model of the viscoelastic model is shown in Figure 11.
A generic viscoelastic model is applied to the proposed finite model since the actual viscoelastic properties of the fracture surfaces were not yet available.However, such generic functions for viscoelasticity of fracture surfaces could help to determine whether or not the model developed in this study is robust enough for its intended use.The fracture surfaces modeled in this demonstration are highly viscoelastic and tend to flow (permanent deformation) over time and this will clearly affect the conductivity of the fracture (Figure 12).
As shown in Figure 13, the deformation in our demonstration case due to the viscoelastic properties of the model rock are much larger than the elastic (initial) deformation at the fracture surfaces.
In this study we developed a novel formulation using the nonlinear finite element method to assess the viscoelasticity of fracture surfaces for the purpose of estimating the conductivity of a propped hydraulic fracture.Nonlinear geometry, contact mechanics between the proppants and the fracture surfaces, and failure analysis were all developed within the framework of the finite element method.Unfortunately, the experimental data measuring the viscoelastic properties of rocks was not available at the time of publication so we have applied generic functions for relaxation and creep functions in order to evaluate the robustness of the developed model.
The model seems to be capable of simulating a variety of processes where viscoelastic behavior is involved and seems suitable for rocks mechanics and fracture deformation applications.However, there are a few additional considerations that could improve the model.The model described here is applied at the micro level and can only consider a few proppant grains on a fracture surface.Therefore, periodic boundary conditions are proposed to apply this technique to a macro-model with more realistic complexity.In addition, computational efficiency needs to be improved in order to reduce CPU time.
There are some additional improvements that could be made to the setup of the model itself.Firstly, in the FE analysis the fracture surfaces were subjected to a uniform pressure load rather than the increasing fracture closure pressure experienced in a real reservoir as pore pressure is depleted.Secondly, the material properties of the rock used in the model were homogeneous and isotropic throughout the entire model which is an oversimplification of realistic rock properties.Quantitative validation of the fracture surface deformation in a dynamic analysis will require the material property data from the biaxial/triaxial tests for each model of a particular viscoelasticity.This must be followed by a comparison of the predicted local strain distribution in the model with experimentally measured strain fields from the same model tested under dynamic loading.
As it stands, the current model formulation promises to allow an improvement in prediction of propped fracture performance by using laboratory and/or log derived material properties.

Conclusion
In this study, a more computationally effective algorithm was applied to use discrete-spectrum approximation to the creep and relaxation functions which allows for adjustment of the energy dissipativity of the model in a significantly low computation-time.As a preliminary demonstration of this model, this numerical technique was applied to an idealized rock fracture model with high energy dissipativity in the rock material.This proposed technique can be used for the evaluation of fracture conductivity under short or long-term loading conditions when rock material properties are considered nonlinear e.g., visco-elastic or visco-plastic.One of the methods to estimate the precise efficiency of a reservoir is the application of multiscale modeling and the present study is one step forward towards this goal.
the exponential integral function and basically .The objective is to solve Equations (

Figure 1 .
Figure 1.The mesh model of the hollow cylindrical model considered in example 1 which consists of 4872 triangular 2 nd order elements.
The model robustly shows the variation of and with respect to the percentage of deformation at different T for the given generic viscoelastic model.Based on the artificial material constants chosen in this study decrease as U increases, however results show that are optimal at moderate deformation and away

Figure 2 .Figure 3 .
Figure 2. The variation of with respect to U at different T, -20, -10, 0, 10 and 20°C, G is in Pa, and U is in %.

Figure 4 .
Figure 4.The simplified model of one proppant placed between fracture surfaces.The material model considered here is linear poroelastic for both proppants and fracture surfaces and the proppant is stiffer twice as much.A three-dimensional model is designed and the contact between the proppants and the proppants and the fracture surfaces are considered nodal (node to node)(Mohammadi et al., 2011).

Figure 5 .
Figure 5.The FE mesh model of the contact problem.For a 3D contact model, extra attention must be paid to provide the most appropriate mesh model specially around the contact area.

Figure 6 .
Figure 6.The von Mises contact stress around the contact area using a nonlinear FE solution.As observed in the diagram the maximum von Mises stress occurs not on the contact surface however in the depth of solid fractures meaning that the fracture and failure starts from the area which is located underneath of the contact area.

Figure 7 .
Figure 7.The shear stress distribution around the contact area using a nonlinear FE solution.As observed in the diagram the maximum shear stress occurs not on the contact surface and not on the midline of the model.There are two identical twins (or eyes) located at depth and offcenter in which the maximum shear stress occurs.

Figure 8 .
Figure 8.The algorithm applied for the failure analysis in this study.

Figure 9 .
Figure 9.An initial failure analysis algorithm applied in this study.The failure initiation is in the area in which contact stresses exceed the UCS of the fracture surfaces.The next area to see failure is at the corner of the interface as indicated in the diagram.The iteration stops when the von Mises stresses are less than the UCS.

Figure 10 .
Figure 10.The new algorithm to model the failure analysis by maintaining the 80% of volumetric fraction of the failed elements.

Figure 11 .
Figure 11.A sketch of the initial FE mesh model and the initial geometry before loading.

Figure 12 .
Figure 12.Creep deformation of fracture surfaces due to constant external force in 3 time-steps.The graphic at far right is the initial geometry and the far left graphic shows the steady state deformation after several creep time steps.

Figure 13 .
Figure 13.A cross section of the model showing the portion of proppants already embedded into the fracture surfaces.