Modeling of Placement of Immiscible Fluids of Different Rheology into

ARTICLE pubs.acs.org/IECR

Modeling of Placement of Immiscible Fluids of Different Rheology into a Hydraulic Fracture Alexander Lakhtychkin,† Oleg Vinogradov,*,† and Dmitry Eskin‡ †

Department of Mechanical and Manufacturing Engineering, University of Calgary, 2500 University Drive NW, Calgary, Alberta T2N 1N4, Canada ‡ DBR Technology Center, Schlumberger, 450 17th Avenue, Edmonton, Alberta T6N 1M9, Canada ABSTRACT: A model of consecutive placements of two different immiscible fluids, possessing different densities and rheological properties, into an expanding fracture is presented. Hydraulic fracturing technology, based on such a fluid placement, is expected to allow better fracture shape control. Mathematically, the problem is reduced to solving a two-dimensional second-order nonlinear partial differential equation describing pressure distribution in a semiellipsoid domain with moving boundaries. The level set method is used to track the interface between the fluids in time. Numerical examples of the placement of power-law fluids are presented, and the accuracy of the computations is investigated.

1. INTRODUCTION Hydraulic fracturing is one of the most popular and efficient methods of oil and gas reservoir stimulation technologies, allowing a significant increase in well production. This technique is based on viscous fluid injection under high pressure into a perforated well that leads to reservoir fracturing. By the end of the treatment, the fracture length can reach several hundred meters and its height tens of meters, while the fracture width does not usually exceed about 1020 mm. Besides fracturing a rock, a fracturing fluid serves as a proppant (solid spherical particles, often sand) carrier. The proppant function is to prevent fracture closure after slurry pumping is stopped. A proppant pack, formed in a closed fracture, represents a high-conductivity channel in a reservoir that can increase oil production by several times compared to production of a nonfractured well. A detailed description of hydraulic fracturing technology is presented.1 Conventionally, the permeability of a fracture and its efficiency are estimated by the total mass of the proppant placed inside the fracture. However, the proppant distribution within the fracture may play an important role. There are at least a couple of potential benefits for the oil and gas industry from a technology allowing proppant placement control. The first one is the creation of large empty channels in a proppant pack filling the fracture. This kind of proppant placement may increase the total fracture conductivity, while reducing the total amount of proppant placed in the fracture. Another possibility is associated with possible fracture shape control. For example, providing a high concentration of small-size proppant in a certain fracture segment may significantly increase the hydrodynamic resistance of that segment, causing suppression of fracture propagation in this direction. This is important, for example, for a thin low-permeable reservoir when a long fracture of a small height is required. The placement of immiscible fluids of different densities and rheologies into a fracture is one of the possible methods for the fracture shape and proppant placement control. The fluids are considered to be pumped in a sequence through the perforated casing, one after another. Every newly pumped fluid pushes the r 2011 American Chemical Society

previous ones deeper into the expanding fracture. Because these fluids have different densities, they also move relative to each other carrying proppant, as is shown schematically in Figure 1. In practice, by variation of the fluid composition, it is possible to vary the fluids densities and viscosities in wide ranges; therefore, theoretically, it is possible to provide the desired fluid properties. For example, if the second fluid is heavier than the first one, then the second settles on the fracture bottom, while the first (lighter) flows to the top. Also, if the second fluid is more viscous than the first one, then the pressure drops faster toward the bottom than toward the top of the fracture, causing a slower fracture growth toward the bottom. Because the fluid properties affect the proppant distribution inside the fracture, the interface between the fluids has to be tracked with significant accuracy. Physically, the problems of the fluid (slurry) flow in the fracture and the fracture propagation are coupled. However, we focus our efforts only on the hydrodynamic part of the problem, i.e., on solving the problem of slurry transport, which by itself represents a big challenge. This decoupling of the real phenomenon can be justified by the following arguments. The first one is that even the solution of the hydrodynamic part of the problem already allows some quantitative estimations of the fracture height growth.2 In this paper, vertical fracture growth was studied when a top or bottom portion of the fracture was filled with a different fluid. The second reason is that the solution of the hydrodynamic problem can be incorporated into the existing codes solving the rock mechanical part of the problem. Decoupling between the hydrodynamic component of the problem and the problem of rock mechanics is achieved by assuming that evolution of the fracture geometry in time is known; i.e., some realistic fracture growth scenario is assumed beforehand. In the current work, we present the case of transport of two immiscible fluids only. However, the Received: September 16, 2010 Accepted: March 9, 2011 Revised: January 18, 2011 Published: March 28, 2011 5774

dx.doi.org/10.1021/ie101914d | Ind. Eng. Chem. Res. 2011, 50, 5774–5782

Industrial & Engineering Chemistry Research

ARTICLE

developed model can be generalized for calculation of three and more fluids. In the future, the developed model can be used as a component of a general model, in which a rock fracture part is incorporated.

2. PROBLEM DESCRIPTION The solution of the problem starts from some “small” preexisting fracture of a known geometry filled with a single fluid (see the left-hand side in Figure 1). At the initial time moment (t = 0), the second fluid starts to enter the fracture, displacing the first one and forcing it to flow deeper into the expanding fracture. The fracture is considered to be planar and vertical (see the illustration of the fracture shape in Figure 2). The fracture width can be represented as an arbitrary, smooth enough function of x and y. However, in the present analysis, the fracture is considered to be an ellipsoid with a cut edge at a certain thickness wmin determining the fracture tip width (see Figure 2b). Using the notations presented in Figure 2, the equation describing the shape of the frontal half of the fracture wall (see Figure 2) is written as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u 2 ! 2 2 ! u w x y min , þ wðx, yÞ ¼ w0 t1 1 L H w0 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi 2ffi x x 0 e x e L, H 1 e y e H 1 L L

where L and H are the fracture half-length and half-height, respectively, w0 is the fracture half-thickness at x = 0 and y = 0, and wmin is the fracture half-thickness at the boundary. In practice, the fracturing fluid does not reach the fracture tip. This effect is called the fluid lag. The value of this lag depends on the fluid rheology and rock mechanical properties. Because we do not model the fracture mechanics (the fracture shape evolution is prescribed), the fluid lag is not taken into account explicitly. Nevertheless, our model takes this effect into account through parameter wmin, which restricts the fracture domain filled with fluid. In principle, the developed model allows incorporation of the fluid lag phenomenon more accurately. The back part of the fracture wall is symmetric to the frontal one with respect to the xy plane. The fracture boundary is comprised of two parts, as shown in Figure 3. The first part, the wellbore-fracture contact line, is split into three sections. The first section is in the middle (of height Hi), where the fluids are injected into the formation through perforations in the wellbore at a known total flow rate Q0(t). Note that within our model we avoid detailed modeling of flowthrough perforations; therefore, the perforated section is represented by a flat linear source. We assume also that the diameter of the wellbore is large enough, so that the pressure along the perforated section changes only hydrostatically. The two other sections have no perforations and, therefore, are characterized by zero normal fluid flux. The second part, the curvilinear moving fracture boundary, is located inside the formation. To define the

ð1Þ

Figure 1. Schematic of fracture growth (the light shading identifies the fluid first pumped into the fracture; the second fluid is shown by the dark shading). It is assumed that 0 < t1 < t2.

Figure 3. Computational domain and boundary conditions.

Figure 2. Hydraulic fracture shape: (a) front view of the fracture; (b) top view of the fracture; (c) front view of the fracture (the contour lines define the fracture width in millimeters, while the fracture length and height are given in meters). 5775

dx.doi.org/10.1021/ie101914d |Ind. Eng. Chem. Res. 2011, 50, 5774–5782


ARTICLE

the times t and t þ Δt, using the constraints defined by eq 2, is approximated as Z Z Lf ðx, yÞ dS V ½LðtÞ, HðtÞ, w0 ðtÞ þ Q0 ðtÞ Δt þ 2Δt S

¼ V ½LðtÞ þ ΔL, HðtÞ þ RΔL, w0 ðtÞ þ βΔL

Figure 4. Determination of the normal component of boundary flux q at point (x0, y0): (a) scheme of the fracture in the xy plane; (b) view of the fracture (x0, y0) (x1, y1) element in a plane normal to the xy plane.

boundary conditions for this part of the fracture, the following procedure has to be followed. The fracture growth scenario is prescribed by the introduction of two constants constraining the fracture height and width expansion rates relative to the fracture length expansion rate: 8 dH dL > < ¼R dt dt ð2Þ dw0 dL > : ¼β dt dt where dL/dt is the fracture length change rate that is found at every time step by using the fracture volume balance equation (see eq 4). The volume of the fracture, shown in Figure 2, is calculated by using eq 1 as 0 1 ffi Z L Z H pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðx=LÞ2 @ V ðL, H, w0 Þ ¼ 4 wðx, yÞ dyA dx ð3Þ 0

0

Once a new fragment of the fracture is open, the fracturing fluid leaks through the fracture faces into the rock. During this process, the polymer component of the fluid forms a thin layer of low permeability on the rock surfaces (the so-called filter cake). The time of formation of this layer is relatively short (on the order of tens of seconds). However, during the filter cake formation, the leak-off rate is significantly higher than that through the formed filter cake, and therefore, in many practical situations, this effect (spurt-loss effect) should be taken into account. Following Carter,3 the total leak-off, Lf, can be introduced as a function of time: Lf = 2CLt1/2 þ Sp, where CL is the leak-off coefficient and Sp is the spurt loss. The spurt loss is measured as the volume flux per unit of newly opened area and depends on both fracturing and the reservoir fluid properties and also on the rock permeability. Note that the spurt loss is inversely proportional to the fracturing fluid viscosity. In the field, values of CL and Sp can be obtained from a mini-frac data (e.g., Fan and Cheng4). The spurt loss value can also be obtained from laboratory tests,5 providing dependence of the leak-off rate versus time. Although we did not account for the spurt-loss effect in the numerical examples presented in this paper, any leak-off function can be easily introduced in our model on the basis of either a laboratory or a field test. For the finite small time increment Δt and the corresponding fracture length increment ΔL, the fluid volume balance between

ð4Þ

where S is the fracture wall surface area. The only unknown, ΔL, in eq 4 is calculated numerically at the beginning of every time step. Then, from eq 2, the values of H and w0 are calculated at time t þ Δt. In accordance to the fracture growth scenario described above, the newly appearing portion of the fracture has to be filled with the fluid coming from the internal part of the fracture, as shown in Figure 4. In this figure, q represents the fluid flux per unit fracture length; in other words, q is defined as vector (vx, vy) integrated along the fracture thickness. The flux qn is the normal component of vector q to the fracture boundary. The corresponding fluid volume balance at the boundary is Z qn ðx0 , y0 Þ 3 Δt ¼ 2 wðx, y, t þ ΔtÞ dΔΓ þ Sp ΔΓ ð5Þ ΔΓ

where ΔΓ is the normal displacement of the fracture boundary during the time step Δt. The leak-off term 2CLt1/2 was not included in eq 5 because our analysis showed that reduction of the time step Δt allows us to make this term small compared to other (included in eq 5) terms. The spurt-loss term can also be neglected in cases when the condition Sp/wmin , 1 is met. We emphasize that, according to Fan et al.,5 this inequality is violated in the case of a highly permeable brine-saturated rock. Note that for illustration purposes in all computational examples presented in this paper we assumed that Sp = 0. Using this approach, it is possible to determine the normal component of the fluid flux qn at the boundary. This flux component qn will be used as the boundary condition for eq 14 formulated below. Therefore, the given total fluid injection rate into the well, Q0, along with the assumption about the prescribed fracture geometry evolution, is transformed into the boundary condition of the Neumann type (imposed on the pressure derivatives) at the curvilinear boundary. The nonperforated section of the well is also characterized by the boundary condition of the Neumann type (∂p/∂x = 0). Because it is assumed that along the perforation line the pressure changes hydrostatically, the absolute pressure at the point (x, y) = (0, 0) can be set equal to any constant (which was zero in our case).

3. GOVERNING EQUATIONS Derivation of the governing equations starts by applying the NavierStokes equation to a fluid flow: Dv þ v 3 rv ¼ rp þ r 3 T þ Fg ð6Þ F Dr where F is the fluid density, v is the fluid velocity, p is the pressure, T is the deviatoric stress tensor with components Tij(x,y,x), and g is the gravity acceleration. Following the recommendations of different authors (e.g., Advani et al.6 and Clifton and Abou-Sayed7), we assumed that fracturing fluids possess power-law rheology.1 The power-law model is simple and provides a relatively good description of the rheological properties of the majority of fracturing fluids in the 5776



ARTICLE

range of shear rates specific for hydraulic fracturing. However, the power-law fluid rheology model predicts infinite fluid viscosity at zero shear rates, which is not true for the real fluids. This contradiction is overcome by the introduction of a threshold viscosity value, μn, such that if the effective viscosity exceeds this limit, the fluid behavior index n is set to unity and the flow consistency index, K, is set to μn. In the developed code, the value μn is set to 2 Pa 3 s. Power-law rheology can be described as μ ¼ K γ_ n 1

inertial Fυ2 n w1 þ n ∼ viscous HK

ð8Þ

where υ is the absolute value of the velocity vector. Substituting regular hydraulic fracturing parameters (F ∼ 1500 kg/m3, H ∼ 50 m, w ∼ 0.01 m, v ∼ 0.3 m/s, K ∼ 3.6 kg/(m 3 s2n), and n ∼ 0.36) into eq 8, we obtain R ∼ 6 103

ð9Þ

To assess the minimum viscosity value needed to justify neglect of the inertial term, a single Newtonian fluid (n = 1 and K equals the fluid viscosity) flowing into the fracture can be considered. Let us assume that the ratio R should be less than 0.1; then using the same parameters (the velocity and fracture geometry) as those previously used, we obtain that the viscosity must be higher than 0.035 Pa 3 s. This criterion is usually satisfied in practice, except in applications in which water is employed as the fracturing fluid. After the simplifications mentioned above are applied to eq 6, we obtain rp Fg ¼ r 3 T

n ! d dvξ K Gξ ¼ dz dz

ð10Þ

The vector G = rp Fg is introduced for calculation convenience. Because we employed the lubrication approximation, which implies a small gradient of the fracture width, the vector G is the only vector variable that may influence the fluid velocity. Thus, the vector G is parallel to the fluid velocity vector v, and we can introduce a local coordinate system, in which these vectors are aligned with some axis ξ. Because power-law fluid rheology was assumed, the shear stress is calculated as n dvξ ð11Þ Tzξ ¼ K dz Thus, the problem is reduced to a one-dimensional flow in a channel of constant width, and eq 10 in this case is reduced to the

ð12Þ

Applying the no-slip condition at the fracture walls to eq 12, we obtain vξ as a function of z. Then the fluid flux qξ is calculated by integration of the equation for vξ along the fracture width: qξ ¼

ð7Þ

where K is the flow consistency index, n is the flow behavior index, γ_ is the shear rate, and μ is the fluid viscosity. To simplify eq 6, we employ the lubrication approximation.8 Because the fracture width is at least 3 orders of magnitude smaller than the other fracture dimensions, we can neglect minor viscous terms that do not contain velocity derivatives with respect to the fracture width. Considering that the fluid viscosity is high and the injection rate is relatively low, we can also neglect both the acceleration and inertial terms in the momentum equation. Moreover, it can be shown that the acceleration and inertial terms are of the same order of magnitude and their ratio to the major viscous term, R, is evaluated as R ¼

ordinary differential equation

2n 2 þ 1=n 1=n 1=n w K Gξ 2n þ 1

ð13Þ

The flux q is written in vector form in the xy coordinate system as " 2 #ð1 nÞ=2n 2n Dp 2 Dp 2 þ 1=n 1=n w þ Fg q ¼ K þ 2n þ 1 Dx Dy Dp Dp ex þ þ Fg ey ð14.1Þ Dx Dy where ex and ey are the unit vectors for the x and y coordinates, respectively, and g is the absolute value of the gravity acceleration. The mass balance equation in terms of the fluid flux is written as Dw þ Lf ð14.2Þ r3q ¼ 2 Dr This equation shows that the fluid is accumulated in a calculation cell because of both the width increase and the fluid filtration into the rock. Note that eqs 14.1 and 14.2 have been successfully used by a number of investigators (see, e.g., Ouyang et al.9). To make the set of eqs 14.1 and 14.2 applicable to the two-fluid system, we assume that the density F and the fluid rheology parameters K and n are functions of the spatial coordinates. These parameters have to be different in the domains filled with different fluids. Substituting eq 14.1 into eq 14.2 and performing routine math, we obtained the governing equation for the pressure distribution p (not presented here because it does not add anything in substance but requires substantial space). This governing equation is classified as a second-order nonlinear partial differential equation (PDE). Because of the assumption of the previously prescribed fracture geometry, the right-hand side of that resultant equation is a known function of time and coordinates. The boundary condition is also a known function of time as discussed above; see eq 5. The obtained formulation of both the equation and the boundary conditions allows resolution of the steady-state problem at any current time step and, as a result, obtainment of the pressure distribution, and thus the fluids velocities at this time. Then, using the calculated velocities, the positions of the fluids at the beginning of the next time step are determined.

4. NUMERICAL TECHNIQUE Because of the complexity of the governing equation and its boundary conditions, the problem has to be solved numerically. The steady-state solution of the nonlinear eq 14 is obtained by using the NewtonRaphson method. Linearization of eq 14 is achieved by representing the pressure distribution as pðx, yÞ ¼ p0 ðx, yÞ þ p1 ðx, yÞ, where p1 ðx, yÞ , p0 ðx, yÞ ð15Þ where p0(x,y) is the assumed initial pressure, which is supposed 5777



ARTICLE

Figure 5. Numerical patterns applied in the problem.

to be close to the real solution, and p1(x,y) is the correction term, which is supposed to be significantly smaller than p0(x,y). By substituting eq 15 into eq 14 and omitting the arguments (x, y), we obtain the following linearized governing equation for p1: Ap1xx þ 2Bp1xy þ Cp1yy þ Dp1x þ Ep1y ¼ F

ð16Þ

where all of the coefficients are functions of p0 and the coordinates. Equation 16 is solved at each iteration step, and then the obtained values of p1 are added to the corresponding values of p0 to obtain a better approximation to the solution; at the next iteration, this summation is used as p0. This iterative process is repeated until the numerical values of p1 reach the machine accuracy. Once this plateau is reached, one can state that convergence is achieved at the current time step. Numerical experiments showed that less than 10 iterations are usually required to reach the machine accuracy. It is clear that a proper initial guess of p0 at t = 0 is important to obtain a rapid convergence of the iterative process. In our case, the initial pressure distribution at the first time step is defined as a polynomial function of x and y, approximately satisfying the boundary condition, namely, at the fluid inlet zone, where the pressure is prescribed, and at three other points: at the upper-left and bottom-left corners of the computational domain and at the rightmost point of the domain (see Figure 3). At the first two corner points, the x component of the fluid flux is zero. By using the boundary condition expressed by eq 5, the y pressure derivatives can be found and then matched during selection of the initial guess of p0. At the rightmost point, the y component of the flux is zero if a fracture is filled with a single fluid, and thus the x pressure derivative can be determined and also used during selection of the initial guess of p0. It is important to emphasize that our calculations showed that the converged solution does not depend on the initial guess. This gives us confidence that the

iterative process converges to the true solution. For any sequential time steps, t > 0, the initial values of p0 are taken as the resultant values of p0 at the previous time step because the pressure distributions within a small time interval are relatively close. Under some conditions, for example, for a small value of the fluid behavior index n, the nonlinear terms of the governing equation reach large values. This requires an improved initial guess for p0 at the first time step. If the initial guess is not sufficiently close to the solution, the iterations may not converge, and thus the solution cannot be determined. In this case, the solution can be found in two steps. For example, if the flow behavior index n = 0.4 and the solution do not converge, we can initially find a solution of the same problem but at n = 0.6 and then, using the obtained solution as the initial guess, solve the original problem at n = 0.4. The numerical experiments performed have shown the robustness of this approach. To solve eq 16, the finite difference method (FDM) of second order with square mesh was employed. The computational domain, grid, and nodes used to determine the numerical pressure derivatives at some special nodes are shown in Figure 5. The first column of nodes is located along the wellbore, while the others are inside the computational domain. To make our numerical scheme more symmetric, we use symmetric patterns at the top and bottom parts of the domain, as indicated by nodes “1” and “2” in Figure 5. The perforated near-wellbore section of the fracture is characterized by the prescribed pressure (the Dirichlet-type boundary conditions) and, therefore, p1 should be constantly zero at this part of the boundary. The nonperforated sections of the wellbore boundary are characterized by zero flux (the Neumanntype boundary conditions). For these nodes, we apply the same scheme as that for the internal nodes of the domain but assuming that the left-hand-side value of p1 equals that on the right-hand side (this procedure automatically satisfies the zero normal flux boundary condition). If one or more nodes, required by the numerical pattern, are located outside the domain, we find the corresponding intersection with the boundary (the A point in Figure 5). The value of p1 at this point p1A cannot be directly included into the computational matrix. However, p1A can be found by using p1 values at five neighboring nodes inside the domain (the “stars” in Figure 5), which are already included in the set of equations. On the one hand, by using these values for the Taylor expansion of the pressure in the vicinity of the A point, we determine the pressure derivatives and the corresponding fluid flux at the A point. On the other hand, the fluid flux at the point A is known from the boundary conditions in eq 5. By equating these two fluxes, we determine p1A as a function of the p1 values at the “stars”, and the value of the normal boundary flux at the A point. Then the obtained expression for p1A can be included in the computational matrix for p1 values. Because the FDM employed for resolution of the governing PDE is based on the second-order Taylor expansion, the pressure gradient should be a continuous function of at least second order. However, because we consider the interface between two immiscible fluids with different properties, the pressure gradient is discontinuous. Formally, this makes a solution based on the Taylor expansion inapplicable. To overcome this problem, a formal transition (mixing) layer of constant thickness at the fluidfluid interface was introduced. Its size was selected to be equal to double the size of the mesh used for the numerical solution. This technique is standard and is routinely used for 5778



ARTICLE

Figure 6. Results of the numerical tests. (a) Numerical solution convergence (the maximum value of p1 in the domain is plotted versus the number of iterations, nI (different lines correspond to different values of the fluid behavior index n). (b) Average deviation of the numerical solution from the analytical one versus the mesh size. The solid line represents the fitted line.

these types of problems.10 The mixing zone does not significantly affect the solution, while enables us to work with the pressure as a smooth function. For example, the fluid density can be calculated as follows: 8 F2 , d < δ > > ! 2 δ 2 δ 2 2 > : F1 , d > δ δ e d e δ

ð17Þ

where δ is the half-thickness of the mixing layer and d = d(x, y) is the signed distance function (a scalar function) defined in the entire domain. The absolute value of the signed distance function at any point equals the distance between the point (x, y) and the interface between the fluids. The sign of this function is positive if the point is located in the first fluid and negative if the point is located in the second fluid. The described numerical algorithm was realized by using Mathematica 7 software.

5. VALIDATION OF THE SPATIAL SOLUTION Because there are no experimental data available to verify the numerical solution, it is important to validate the solution indirectly. First of all, it was important to test the reliability of implementation of the numerical method by investigating the solution convergence, and especially the convergence rate. We expected that the closer we approach the real solution, the better the linear approximation should become and the faster the solution should converge. The maximum value of p1 in the domain was detected and plotted after each performed iteration for the following parameters of the system: F1 = 1000 kg/m3, F2 = 1500 kg/m3, K1 = 0.5 Pa 3 sn, and K2 = 1 Pa 3 sn. The values of n were varied but maintained the same for both fluids. The results presented in Figure 6a confirm the expectations formulated above. It may also be noted that for the Newtonian fluids, i.e., n1 = n2 = 1, the solution should converge in a single step because of the linearity of the governing equation in this case. In Figure 6a, one can see that the corresponding curve reaches the machine accuracy at the second iteration, thus proving the validity of the code. Another possible way to validate the applied numerical scheme is to compare the numerical solution against an analytical

one for a particular case. We were able to find an analytical solution for a test problem after making a number of simplifying assumptions: both fluids are identical and Newtonian and the fracture has a penny-shaped form with radius R0 and constant width. In this case, the governing equation (eq 14.114) is reduced to the Laplace equation, which in polar coordinates (r, θ) is written as 1 D Dp 1 D2 p r ð18Þ þ 2 2 ¼0 r Dr Dr r Dθ The next step is to prescribe zero pressure at the entire nearwellbore boundary and set the fluid flux at the circular boundary as qðr ¼ R0 , θÞ ¼

nHarm

∑

j¼1

aj cos½ð2j 1Þθ

ð19Þ

where nHarm is the number of harmonics in the series and Qj are arbitrary coefficients. Under these conditions, the solution can be found analytically by using the method of separation of variables.11 The solution is given as 2j 1 3 KR0 nHarm aj r pðr, θÞ ¼ cos½ð2j 1Þθ 2 w3 j ¼ 1 2j 1 R0

∑

ð20Þ To test the numerical solution, nHarm was set to 2, Q1 = 1, and Q2 = 0.2. The fracture size was L = H = R0 = 20 m and w0 = 1 cm. It has been found that the discrepancy between the analytical and numerical solutions asymptotically approaches zero with mesh refinement. The average deviation between the numerical solution and the analytical one, given by eq 21, versus the mesh size is shown in Figure 6b. The average deviation is calculated as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u i i 2 u ti 1 ðpnum pan Þ ð21Þ average deviation ¼ N

∑

where pinum is the numerically calculated value of the pressure at the ith domain node, pinum is the analytically obtained value of the pressure at the ith node, and N is the total number of nodes in the domain. Note that the fitted curve in Figure 6b is close to the quadratic function. It is also possible to see that, even for a relatively large 5779



ARTICLE

Table 1. Discrepancy between the Total Incoming Flow through the Perforations and Total Boundary Flow at Different Time Steps time, min

0

5

10

15

20

25

30

35

40

45

50

55

60

total flux error, %

0.93

1.38

2.68

2.19

2.55

1.53

2.43

2.03

3.41

2.40

0.55

1.06

1.64

where velocity v is defined as v(x,y) = q(x,y)/2w(x,y).

rj ¼ 1 )

6. TIME-DEPENDENT SOLUTION After obtaining the pressure distribution at a given time step, we determine the fluid velocities inside the fracture and use them within the level set method (LSM), described below, to calculate the interface position for the next time step. Besides this, the computational mesh is updated because of the domain expansion and all functions defined in the domain are extrapolated to the new nodes. After the mentioned updates are completed, eq 16 is solved again at a new time step. Special care must be taken to track an interface between fluids. According to the literature, the interface is usually tracked by using one of three methods: the marker method12,13 the volumeof-fluid (VOF) method,14,15 and the LSM.16,17 As input, all of these methods require a known velocity distribution in a simulation domain, and as output, these methods provide the location of the interface between the two fluids at the next time step. In this work, we employed the LSM because it is considered to be optimal for such kinds of problems. In contrast to Marker-like methods, this method is based on tracking of the fluids themselves rather than the interface between them. This allows one avoid topological complexities, such as interface self-intersection, interface interaction with the domain boundary, and problems caused by excessive grouping or spreading markers. In comparison with the VOF method of the same order of accuracy, the LSM is computationally more efficient and easier in implementation. In addition to this, the LSM allows one to find the distance to the interface from any point of the domain in a simple and fast manner. The idea of the LSM is to set and then propagate in time the so-called “level function” j. This is a continuous function determined at the nodes of the computational domain, and it is positive in one fluid and negative in the other. Therefore, the zero of this function corresponds to the interface between the fluids. In our case, as a level set function, it is convenient to use the signed distance function (see eq 17). This allows one to avoid determination of the distance to the interface for every node, which is a very costly computational procedure. In our problem, the initial setting of the level function is done in a trivial manner because the initial condition is defined as a fracture filled with the first fluid only, while the second fluid is at the fracture entrance. Therefore, setting the level set function at the initial time moment is based on calculation of the distance from each point of the initial domain to the vertical boundary section representing the fracture entry. Further propagation of the level function is based on resolution of the following transport equation: FDM Dj þ v 3 rj ¼ 0 w jni, jþ 1 ¼ ðjni, j v 3 rjni, j ÞΔt ð22Þ Dr

The high computational speed of the LSM in our case is achieved by using the same numerical scheme for determining the derivatives of the level function j as that employed for solving eq 16. It is well-known18,19 that with time the level function loses its initial property of being the distance function. Therefore, at every time step, the level function has to be corrected. The correction is based on the solution of the following equation: )

mesh size of 2 m, the deviation of the numerical solution from the analytical one is about 500 Pa, which is significantly smaller than the average value of the absolute pressure (about 107 Pa). Another possibility of testing the validity of the linearized solution is to substitute the solution p0, found after a sufficient number of iterations, into the nonlinear equation and to check whether the equation is satisfied. This test was also successfully accomplished.

ð23Þ

Equation 23 is a typical Eikonal equation. This equation is solved by the fast marching method, which allows one to find the solution in O(N log N), the total number of domain nodes. A detailed description of the fast marching method is given in the literature.18

7. COMPUTATIONAL RESULTS AND DISCUSSION The developed code allows calculation of the fluid motion in a fracture. To test the accuracy of the developed code, we solved a test problem at the two different time steps (30 and 60 s). Computer simulation of the fluid displacement performed for 20 min showed that there was almost no visual difference between the fluid interface locations obtained in these two cases. Another test verifying the solution was made by a comparison of the fluxes flowing into and out of the fracture. Setting the fracture width to a constant and the leak-off to zero eliminates the right-hand side of eq 14.2, and therefore at any time the flux entering the fracture through the perforations has to be equal to the flux leaving the fracture through the expanding edge of the fracture boundary. The computational test was performed for the constant incoming flux of 102 m3/s. The mesh size was set to 1.4 m and the time step to 20 s. The fluid parameters used in the calculations were F1 = 1050 kg/m3, K1 = 2 Pa sn, n1 = 0.7, F2 = 1000 kg/m3, K2 = 1.5 Pa sn, and n2 = 0.9, where subscript 1 indicates the first fluid, which initially was inside the fracture, and subscript 2 indicates the fluid entering the fracture. The relative differences between the incoming and outgoing fluxes are given in Table 1 for different values of the fracture propagation time. The table shows that the inlet flux is rather close (the deviations are below a couple of percents) to the fracture volume increase rate. Thus, the fluid volume is accurately conserved within the numerical code developed. To demonstrate the numerical code performance, the results of the following simulations are presented in Figure 7. The fluid parameters used in the calculations were F1 = 1100 kg/m3, K1 = 1.5 m 3 s2n, n1 = 0.8, F2 = 1000 kg/m3, K2 = 3.6 Pa sn, and n2 = 0.36. The fracture parameters were the initial half-length L = 20 m, the initial half-height H = 20 m, the injection height Hi = 20 m, w0 = 1 cm, and wmin = 0.5 cm, the incoming flux was 102 m3/s, and the coefficients in eq 2 were: R = 0.2 and β = 0.0001. The mesh size h was 1.3 m, and the time step Δt was 15 s. The fluid interface boundaries at different time moments are shown in Figure 7. The results show that the developed technique allows one to compute the flows of two non-Newtonian fluids in a fracture. 5780



ARTICLE

Figure 7. Simulation of two non-Newtonian fluid propagations in a fracture. The pictures are taken at the times 0, 20, 40, and 60 min (the fracture dimensions are in meters).

The model of placement of two fluids into a fracture demonstrated a good performance. Thus, it can be applied to modeling of the very practically important problem of placement of NonNewtonian slurries. The density and viscosity of the slurry can be easily evaluated if the local volume concentration of the particles is known.20 However, the solid concentration distribution in a fracture continuously changes because of particle settling caused by gravity. Particles may also migrate through the fluid interface boundary. To solve the problem of placement of two slurries, the set of equations formulated above should be supplemented with the mass conservation equation for the solid phase.1 The equations of fluid motion in this case are coupled with the mass conservation equation for solids through the proppant concentration distribution over the fracture determining distribution of the slurry viscosity and density.

8. CONCLUSIONS A mathematical model of a flow of two non-Newtonian immiscible fluids inside an expanding fracture has been developed. The accuracy of the numerical implementation of the model based on the FDM has been proven by investigation of the numerical scheme convergence, by a comparison of the numerical results with an analytical solution available for a simplified geometry, and by investigation of the calculated flux balance. Evolution of the fluidfluid interface in a growing time fracture has been illustrated by numerical examples. The developed model is intended to be used as a basis for the development of a simulator for placement of multiple slurries containing proppant in a fracture. ’ AUTHOR INFORMATION Corresponding Author

*Tel.: 1 403 220 7187. E-mail: [email protected].

9. ACKNOWLEDGMENT The financial assistance provided by the Schlumberger Canada Ltd. and MITACS is gratefully acknowledged. We are also thankful to Dr. Tony Settari (University of Calgary) for valuable discussions. We express our gratitude to Dr. J. R. A. Pearson (Schlumberger Cambridge Research) for helpful critical remarks and constructive advice. ’ NOMENCLATURE Abbreviations

FDM LSM

finite difference method level set method

PDE VOF

partial differential equation volume-of-fluid method

Symbols

A, B, C, D, E, F coefficients of eq 16, (Pa2/m2, Pa2/m2, Pa2/m2, Pa2/m3, Pa2/m3, Pa3/m4) arbitrary coefficients, m2/s aj fluid loss coefficient, m/s C2 d signed distance function, m unit vectors parallel to coordinate system axes, dimenex, ey sionless G vector of modified pressure gradient, Pa/m g vector of gravity acceleration, m/s2 g gravity acceleration value, m/s2 H fracture half-height, m height of the perforated interval of the well, m Hi i domain node number, dimensionless K flow consistency index, Pa 3 sn L fracture half-length, m fluid leak-off through the porous fracture faces, m/s Lf N number of nodes in the simulation domain, dimensionless n flow behavior index, dimensionless nHarm number of harmonics in a series, dimensionless p pressure at a point (x, y), Pa p0(x,y) guess pressure distribution in the fracture, Pa p1(x,y) correcting pressure distribution in the fracture, Pa value of the correcting pressure at the fracture p1A boundary, Pa Q0(t) total flux into the considered half of the fracture, m2/s q vector of the linear flux of the fluid in the fracture, m2/s normal component of vector q to the fracture boundary qn at the boundary, m2/s R ratio of the inertial and viscous terms, dimensionless penny-shapes fracture radius, dimensionless R0 r radius in the polar coordinates, dimensionless spurt-loss coefficient, m Sp T deviatoric stress tensor with the components Tij(i,j = x, y,x), Pa t time since the second fluid started to enter the domain, s V volume of the considered half-fracture, m3 v vector of the fluid velocity equal to (vx, vy), m/s w fracture half-thickness at a point (x, y), m fracture half-thickness at point x = 0 and y = 0, m w0 fracture half-thickness at the fracture boundary, m wmin (x, y, z) spatial coordinates, m (x0, y0) some point at the fracture boundary, m 5781


Industrial & Engineering Chemistry Research Greek Symbols

R β γ_ δ

ΔΓ Δt ΔL θ F μ μh ξ v j

fracture height growth constraint, dimensionless fracture thickness growth constraint, dimensionless shear rate, 1/s half-thickness of the mixing layer between the different fluids, m normal displacement of the fracture boundary during the time step Δt, m time increment, s fracture length increment, m angle in the polar coordinates, dimensionless fluid density at a point (x, y), kg/m3 fluid viscosity, Pa 3 s threshold viscosity value, Pa 3 s coordinate axis aligned with v, m absolute value of the velocity vector, m/s level set function, m

ARTICLE

(18) Sethian, J. Level Set Methods and Fast Marching Method; Cambridge University Press: Cambridge, U.K., 1999. (19) Sethian, J. A.; Vladimirsky, A. Fast methods for the Eikonal and related HamiltonJacobi equations on unstructured meshes. Proc. Natl. Acad. Sci. U.S.A. 2000, 97, 5699–5703. (20) Schook, C.; Roco, M. Slurry Flow; Butterworth-Heinemann: Boston, 1991.

’ REFERENCES (1) Economides, M. J.; Nolte, K. G. Reservoir Stimulation, 3rd ed.; John Wiley & Sons: New York, 2000. (2) Settari, A. Quantitative Analysis of Factors Influencing Vertical and Lateral Fracture Growth, SPE 13862. SPE Prod. Eng. 1988, 310– 322. (3) Carter, R. D.Derivation of the General Equation for Estimating the Extent of the Fractured Area. Optimum Fluid Characteristics for Fracture Extension; Drilling and Production Practice; American Petroleum Institute: New York, 1957; pp 261269 (4) Gadiyar, B.; Morales, R.; Norman., W. Is Spurt Loss a Reality during Frac/Packing in High Permeability Formaitons? SPE 73767. SPE Prod. Eng. 2002. (5) Fan, Y.; Chen, Z.; Rapid, A. Method for Determining FracturingFluid Leakoff Coefficient and Spurt Loss; SPE 3830; SPE: Long Beach, CA, 1997. (6) Advani, S. H.; Lee, T. S.; Lee, J. K. Three-dimensional modeling of hydraulic fractures in layered media. Part I. Finite element formulations. J. Energy Resour. Technol., Trans. ASME 1990, 112, 1–9. (7) Clifton, R. J.; Abou-Sayed, A. S. A Variational Approach to the Prediction of the Three-Dimensional Geometry of Hydraulic Fractures, SPE 9879. SPE/DOE Low Permeability Gas 1981. (8) Pearson, J. R. On suspension transport in a fracture: framework for a global model. J. Non-Newtonian Fluid Mech. 1994, 54, 503–513. (9) Ouyang, S.; Carey, G. F.; Yew, C. H. An adaptive finite element scheme for hydraulic fracturing with proppant transport. Int. J. Numer. Methods Fluids 1997, 24, 645–670. (10) Hammond, P. S. Settling and slumping in a Newtonian slurry, and implications for proppant placement during hydraulic fracturing of gas wells. Chem. Eng. Sci. 1995, 50, 3247–3260. (11) Evans, L. C. Partial Differential Equations. Am. Math. Soc. 1998. (12) Zabusky, N. J.; Overman, E. A. Regularization of Contour Dynamical Algorithms. I. Tangential Regularization. J. Comput. Phys. 1983, 52, 351–373. (13) Harlow, F. H.; Welch, J. E. Numerical calculation of timedependent viscous incompressible flow of fluid with free surface. Phys. Fluids 1965, 8, 2182–2189. (14) Pilliod, J. E.; Puckett, E. G. Second-order accurate volume-offluid algorithms for tracking material interfaces. J. Comput. Phys. 2004, 199, 465–502. (15) Hirt, C. W.; Nichols, B. D. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 1981, 39, 201–25. (16) Osher, S.; Fedkiw, R. P. Level set methods: an overview and some recent results. J. Comput. Phys. 2001, 169, 463–502. (17) Osher, S.; Sethian, J. A. Fronts Propagating with CurvatureDependent Speed: Algorithms Based on HamiltonJacobi Formulations. J. Comput. Phys. 1988, 79, 12–49. 5782


Modeling of Placement of Immiscible Fluids of Different Rheology into

Recommend Documents