Patent Description:
The purpose of stimulation is to the enhance productivity of the well. A common stimulation method for carbonate reservoirs is acid stimulation whereby the selected acid is allowed to chemically react with the reservoir rock, which leads to dissolution and enhanced productivity.

For wells completed open-hole, a complicating factor is the acid placement, i.e. the ability to distribute acid across the entire reservoir section. Bull-heading acid from the surface typically results in a mediocre stimulation treatment because the majority of the acid is spent reacting at the heel of the well.

The solution to the acid placement challenge is addressed by the limited-entry liner (LEL) technique, also denoted as controlled-acid jetting (CAJ). The key concept is to distribute small holes of varying sizes and frequency in the liner. These holes act as flow restrictions, which leads to mechanical diversion of flow along the liner. An appropriate hole size distribution design is capable of ensuring that the entire reservoir section is treated with acid. Aspects of the hole size distribution has been addressed in a number of references.

A further complicating factor is to ensure maximum acid penetration into the reservoir rock. Acid is an expensive commodity and should not be spent on dissolving all the rock in the near-wellbore area. Rather, the stimulation programme should be designed in such a way that acid penetrates as furthest as possible into the formation because this situation leads to the highest negative skin and hence the highest productivity index.

Lab experiments by a number of authors clearly show that for any given rock, acid penetration depends on the interstitial velocity of acid. There exists an optimum velocity, which minimizes the amount of acid needed to generate dissolution patterns known as wormholes. This optimum velocity depends on the rock, and the acid system (type, concentration, temperature). In addition to ensuring uniform acid coverage, the hole-size distribution must also be designed in such a way that it maximizes the propagation of wormholes.

Problem pertains to acid stimulation of both vertical and horizontal wells. The challenge is to achieve uniform stimulation throughout the completed well trajectory. Some operators choose not to stimulate the wells, others bullhead from the wellhead, others stimulate through a coiled tubing. Segmented completions which allow acidization in stages and use of diverters is employed. A few operators make use of the Limited Entry Liner (LEL) concept, but does not describe a comprehensive workflow for the hole size design. The design of LEL in terms of varying hole sizes and frequency remain a subject matter of challenge because of multiplicity of considerations.

<CIT>, the authors describe the LEL concept (called controlled acid jet) for matrix-acid stimulation and develop a very simple steady-state model using polynomial approximation with orthogonal collocation. However, their model assumes constant friction factor and does not describe a workflow for design of the optimum hole size distribution. Their model does not estimate the maximum design rate, does not take into account the experimental wormhole curve, does not have a skin model and is unable to estimate the required acid coverage and the optimum distance between holes.

<CIT>, the authors model a limited-entry liner with a transient model. The flow equations are the same as in the current work, but their numerical solution scheme uses the finite-difference method to solve for pressure and rate in the liner whereas as in this work, polynomial approximation with orthogonal collocation is used. Their model is transient, whereas the current work describes a steady-state solution. Their workflow does not allow design of the actual hole size distribution but is used to analyse data from an existing stimulation job.

The teachings of the prior art does not solve the problem. For example the assumption of constant friction factor in <CIT> does not bode well for the actual phenomenon. They have also not laid out design basis for optimum hole size distribution. It cannot estimate maximum design rate and required acid coverage. The work did not incorporate experimental wormhole curves and skin model. <CIT> used a transient model. The teachings also did not describe optimum hole size distribution.

<CIT> discloses a controller for controlling an apparatus for performing a wellbore intervention or process, a corresponding processing device and method for simulating and/or controlling consecutive flow of a plurality of fluids in a wellbore of arbitrary geometry.

The invention, which will be described in further detail in subsequent paragraphs, consists of a newly developed method and system for stimulating a well, which addresses, at least partly, the afore-mentioned challenges in a novel and inventive way.

It is thus an object of the present invention to provide an accurate and efficient numerical solution strategy for providing an initial estimate of the number of holes per segment which honours the acid coverage per segment and the drop in pressure (dp) across the last hole, in particular in the context of acid stimulation of wells completed in a carbonate reservoir with a Limited-Entry-Liner or LEL liner. The invention is defined by the method of independent claim <NUM>, the data processing system of claim <NUM> and the computer program product of claim <NUM>. Other aspects of the invention are defined by the dependent claims.

According to a first aspect of the invention the accuracy of simulations of fluid transport in a system for stimulating a well in a material formation of a resource reservoir can significantly be improved by including a workflow for design of the optimum hole-size distribution. Therefore optimized hole-size distribution in the liner of a LEL liner system is modelled, which results in an improved modelling accuracy and providing an improved construction and operation of the stimulation system. In particular, providing an initial estimate of the number of holes per segment which honours the acid coverage per segment and the dp across the last hole, such that the initial estimate can be found from the relationship between interstitial velocity, pump rate, and total cross-sectional hole area for a particular discharge coefficient and liner configuration.

According to some embodiments, it is a further object to improve the accuracy of the simulation by ensuring that the annulus pressure remains below fracturing pressure. The maximum allowed pump rate is dictated by the permeability, the fluid viscosity, the length of the completed interval, the skin, and the difference between annulus pressure and reservoir pressure.

According to some embodiments, it is a further objective to improve the accuracy of the simulation by estimating wormholing characteristics to facilitate an optimal hole-size distribution.

The wormholing estimate includes a nodal analysis calculation performed to estimate the downhole temperature at the heel of the liner, and based on the choice of the acid, the permeability and the temperature, the optimum velocity for wormhole propagation is estimated together with the anticipated pore volume to breakthrough.

According to some embodiments, the model accuracy has been improved by providing a method comprising estimation of the total number of holes and dp in pressure across the last hole. Based on the optimum velocity and the calculated design pump rate, the total cross-sectional area of the holes is calculated, wherein the area is linearly correlated with the dp across the last hole.

According to another aspect, a data processing system is configured to perform the steps of the method described herein.

According to yet another aspect, the invention relates to a method of stimulating a well by means of a workflow system for adjusting the hole-size distribution which honours the acid coverage per segment and the drop in pressure (dp) across the last hole in the context of acid stimulation of wells completed in a carbonate reservoir with a LEL liner. The method comprises:.

In one embodiment the model describes a first order non-linear boundary value problem and consists of two first-order coupled non-linear differential equations with boundary conditions at each end of the interval of the independent variable, coupled with a set of algebraic equations, wherein polynomial approximation is used to solve the numerical model.

According to yet another aspect, the method of stimulating a well by means of a workflow system for adjusting the hole-size distribution which honours the acid coverage per segment and the drop in pressure (dp) across the last hole in the context of acid stimulation of wells completed in a carbonate reservoir with a LEL liner in which the constraints of;.

According to yet another aspect, the method of stimulating a well by means of a workflow system for adjusting the hole-size distribution which honours the acid coverage per segment and the drop in pressure (dp) across the last hole in the context of acid stimulation of wells completed in a carbonate reservoir with a LEL liner. The method comprises further:.

According to yet another aspect of the invention, a data processing system is configured to perform the steps of the method of stimulating a well as described herein.

The term data processing system includes any electronic system or device having a processor configured to perform the step of the method, and to communicate the outcome of those steps to a user of the system or device. Such system or device includes, but is not limited to, a computer, a laptop, a handheld electronic device, or electronic workstation.

These and other features of the invention will become more apparent by the following description of the embodiment, which is made by way of example, with reference to the accompanying drawings in which:.

The limited-entry liner consists of a number of unevenly spaced holes with the purpose to distribute fluid, in this case acid, evenly along the reservoir section to be stimulated. The concept was initially described in <NUM> by Shell for fracturing applications (Lagrone and Rasmussen, <NUM>) and is still widely applied. It was later adapted for matrix-acid stimulation and patented by Maersk Oil (known as controlled acid jetting or CAJ) and implemented in North Sea chalk reservoirs on a large scale, see Hansen (<NUM>) and Hansen and Nederveen (<NUM>). Since then, this novel stimulation concept has been tested by various operators such as ConocoPhilips (Furui et al. , 2010a,b), Petrobras (Fernandes et al. , <NUM>), ExxonMobil (Sau et al. , <NUM>; Troshko et al. , <NUM>), ZADCO (Issa et al. , <NUM>) among others (Mitchell et al. , <NUM>; van Domelen et al. , <NUM>, <NUM>). Rodrigues et al. (<NUM>) provided a good general overview of stimulation techniques for low-permeability reservoirs and Shokry (<NUM>) described the acid stimulation practice in ADNOC for offshore reservoirs.

<FIG> shows a schematic cross-sectional view of a well-bore <NUM>. The well-bore <NUM> is conventionally formed by techniques commonly known in the art, and includes a wall <NUM> created by the drilling process, a leading end <NUM>, which extends into the formation <NUM>, and a trailing end <NUM> for accessing the well-bore.

A limited-entry liner <NUM> is introduced into the well-bore <NUM>. The liner <NUM> has an open end <NUM> and opposed sealed end <NUM>. An annulus <NUM> is formed between the wall <NUM> and outer surface <NUM> of the liner.

The liner <NUM> is provided with a number of pre-formed holes <NUM> that form flow passages between the interior of the liner <NUM> and the annular space <NUM>. The holes <NUM> have a shape and location that comply with particular, pre-defined specifications.

Typically, the distances between adjacent holes <NUM> along the liner <NUM> decrease towards the end <NUM> of the liner.

The acid is pumped into the liner in the liner <NUM> and exits holes <NUM> at high velocities resulting in jetting into the formation <NUM>. By limiting the number and size of holes, a choke effect is obtained and a significant pressure drop occurs between the inside and the outside of the liner during stimulation. A non-uniform geometric distribution of the holes is used to compensate for the friction pressure drop along the liner section. This means that the average hole spacing decreases towards the bottom of the liner. The open annulus <NUM> outside the liner in combination with the overpressure on the inside of the liner (due to the choking over the holes) ensures that the acid eventually reaches the bottom of liner, and the well is thus stimulated along its full length.

Acid is bull-headed from the surface and enters the liner <NUM> in the direction of arrows <NUM>. The liner does not have to be horizontal but very often is. When acid reaches the first hole <NUM>, which has a size of <NUM>-<NUM>, the pressure drop across the hole is so high that only a small portion of the acid exits the liner through the hole; the remaining portion continues along the liner until it reaches the next hole where the same process is repeated. An appropriate hole-size design makes it possible to honour a specified acid coverage, defined as barrels of acid per feet of reservoir section. Prior to the stimulation, the mud can be circulated out so that only completion brine with the right density is found in the wellbore <NUM>.

The acid stimulation process is modelled by discretizing the wellbore <NUM> into a number of nodes <NUM>, typically <NUM>-<NUM>. The nodes do not need to have the same size. From a practical design point of view, the wellbore is split into a smaller number of segments <NUM>. These segments may be physically isolated from each other on the annulus side by hydraulic packers <NUM> (not shown) but do not have to. Nodes can overlap between two segments, as shown in <FIG>.

Displacement of brine by acid is considered to occur by single-phase plug flow with minimal dispersion. The negative excess mixing volume is not taken into account. The liner <NUM> is closed at sealed end <NUM> before stimulation and it is not cemented, which means that fluid can in principle flow in the annulus <NUM> along the well-bore trajectory before packers <NUM> are set. In practice, annulus flow occurs predominantly due to jetting of acid through the holes <NUM>, perpendicular to the wellbore. Annulus flow along the liner can be ignored for practical modelling purposes.

The completion design, and the associated modelling workflow covered in this document, allows for reservoir segmentation using packers and the resulting liner is hence referred to as a segmented limited-entry liner. The desired acid coverage can be specified per segment to take into account differences in porosity, permeability, initial water saturation, and reservoir pressure. The number of segments for modelling the process can be larger than the number of packer-isolated intervals.

Design of the hole-size distribution depends primarily on liner geometry and flow rate, which in turn is governed by reservoir properties, i.e. reservoir permeability. Acid stimulation is inherently transient in nature because the skin factor at any given position along the well changes with time from a positive value initially (caused by a mud filter cake) towards a negative value once the acid has reacted with the reservoir rock minerals. If the skin evolution over time is uniform along the well, it will not affect the flow distribution, which means that the overall process can be modelled based on steady-state principles.

The invention consists of a comprehensive algorithm for designing the hole-size distribution for limited-entry liners. The next sections describe the algorithm for designing a hole-size distribution which achieves a specified (often uniform) distribution of acid volume per interval length, also known as acid coverage.

The algorithm is shown schematically in <FIG> and <FIG>. <FIG> shows the overall algorithm whereas <FIG> shows a more detailed part of <FIG>. The algorithm is discussed by reference now to <FIG> and the first block, input data and constraints <NUM>.

As a starting point for implementation of the algorithm, input data constraints are entered into the system. The input data is made up of rock properties, completion data, fluid properties and other data, such as pump rate, number of nodes for the numerical algorithm, pressure drop across the last hole of the liner and annulus pressure. These inputs are either known, or may be sourced from historical data from the wellbore.

The algorithm defines certain constraints which must be adhered to in the functioning of the system. These constraints form part of the input data and constraints <NUM>. The constraints includes, but are not limited to; annulus pressure must exceed minimum reservoir pressure to avoid cross-flow inside wellbore; annulus pressure must not exceed fracturing pressure to avoid fracturing; wellhead pressure must not exceed maximum design pressure rating - in turn this impacts the design rate and/or the amount of friction reducer to be added; cross-sectional area of all LEL holes combined should be equal to or exceed a minimum cross-sectional area to avoid creating an additional pressure drop during normal production or injection of the well after stimulation - this impacts the number and size of the holes; average distance between two neighbouring LEL holes should equal twice the wormhole radius - this impacts the pressure drop across the last LEL hole which is a design variable; and, the liner ID cannot exceed the wellbore size.

Based on the input per segment, the maximum rate per segment is found by applying the transient inflow equation. Note that although the well is horizontal, it acts as a vertical well in the early injection phase because the boundaries have not been felt. Hence, the reservoir section length, L, replaces the reservoir thickness, H.

B is the acid formation volume factor, which is in the range <NUM> to <NUM>. In practice, it is assumed to be <NUM>. The viscosity is the maximum value of the oil or gas viscosity and the acid viscosity. In heavy oil reservoirs, the transient phase injectivity is initially controlled by the oil properties. Thus,
<MAT>.

The permeability will see a contribution from the two horizontal directions as well as the vertical direction:
<MAT>.

The vertical/horizontal permeability ratio may attain values in the range <NUM> to <NUM>. For the current application, the value is close to <NUM>, which makes the overall permeability equal to the horizontal permeability.

Where the total system compressibility is given as a contribution from the rock and the fluid phases present in the pore space.

rw refers to the wellbore radius. In gas reservoirs, co equals gas compressibility.

The maximum pump rate allowed is then the sum of the individual segment rates:
<MAT>.

However, any segments which must be left unstimulated and therefore require joints without holes, do not contribute to the calculation of the total rate. To start the design algorithm detailed later, the actual design rate is taken as a value <NUM>-<NUM>% lower than the maximum allowed rate. This value may be adjusted in a subsequent iteration.

T is the total pump time calculated from the acid coverage and length of all the segments
<MAT>.

It is noted that T depends on Q, which depends on T.

Research into matrix-acid stimulation fundamentals took off in the <NUM>'s with the pioneering work of Fogler and co-workers from the University of Michigan (Hoefner et al. , <NUM>; Hoefner and Fogler, <NUM>; Bernadiner et al. , <NUM>; Fredd and Fogler, <NUM>, <NUM>, <NUM>; Fredd et al. , <NUM>) who demonstrated that the acid reaction with the rock gives rise to different etching patterns depending on the type and concentration of acid as well as the velocity and the temperature. Key subsequent contributions in the literature to the current understanding includes work by Halliburton (Gdanski and Norman, <NUM>; Gdanski and van Domelen, <NUM>; Gdanski, <NUM>), Buijse and Glasbergen (<NUM>), and Hill and coworkers from Texas A&M University (Al-Ghamdi et al. , <NUM>; Dong et al. , <NUM>, <NUM>; Dubetz et al. , <NUM>; Etten et al. , <NUM>; Furui et al. , <NUM>, <NUM>, 2010a,b; Izgec et al. , <NUM>; Ndonhong et al. , <NUM>, <NUM>; Sasongko et al. , <NUM>; Schwalbert et al. , <NUM>; Shirley et al. , <NUM>; Shukla et al. Further references to experimental and theoretical studies on wormhole growth are listed within these references.

<FIG> shows the effect of rate on dissolution through a series of images <NUM>. A low rate leads to uniform dissolution and hence a very inefficient usage of the acid. This is shown by the image to the far left <NUM>. In the image the acid <NUM> has not permeated the formation <NUM> to any appreciable extent. At slightly higher rates (i.e. moving from left to right in the images), the acid creates wormholes <NUM> through the rock. In fact, any acid formulation has an optimum velocity at which the least volume of acid is required to etch a pattern from inlet to outlet. This volume is called the pore volume to breakthrough <NUM>. Note that <NUM>% HCl corresponds to <NUM>, hence the <NUM> concentration used in the experiment is quite low.

<FIG> illustrates the impact of interstitial velocity <NUM> on pore volume to breakthrough <NUM> at two different temperatures 204A (depicted by the dot-dash line) and 204B (depicted by the solid line). A temperature increase (i.e. from temperature 204A at 25oC to temperature 204B at 600C) leads to higher reaction rate and hence faster dissolution; optimum wormhole growth therefore requires a higher acid velocity to avoid spending all the acid near the wellbore. It is also clear that it is better to pump at a rate which is slightly above the optimal than below. In a low-permeability reservoir, the maximum pump rate is limited by the fracturing pressure, which may prevent the operator from reaching the optimum velocity. In such situations, it is necessary to select a different acid <NUM> formulation to shift the curve to the left and preferably also down.

The wormhole data can be reproduced with a model proposed by Buijse and Glasbergen (<NUM>) containing two fitting constants, αand β, which can be reformulated in terms of the lowest point on the curve (optimum interstitial velocity <NUM>, optimum pore volume to breakthrough <NUM>)
<MAT>.

Increased temperature <NUM> and increased HCl concentration both increase the optimum velocity <NUM> for wormholing. For low-permeability rocks where the optimum rate may be limited by the fracture propagation pressure, it may be beneficial to reduce the acid concentration, although the pore volume to breakthrough <NUM> increases and hence the volume of acid solution needed. If the acid concentration is halved then the volume must double to maintain the same number of moles. Several authors have investigated the effect of weaker acids, see Punnapala et al. (<NUM>) and Shirley et al. A friction reducer may shift the PV curve upwards, which means that more acid is required to achieve the same skin.

Talbot and Gdanski (<NUM>) proposed a general wormhole model where they correlate the two input parameters to the Buijse-Glasbergen model as a function of rock and acid properties as well as temperature. However, they do not specify the values of the constants in their correlation.

In this invention, we make use of a concept whereby we shift the default wormhole curve shown in <FIG> up, down, left, or right as a function of the temperature <NUM>, the permeability, and the acid type. Table <NUM> shows some rough rules-of-thumb when adjusting the optimum (lowest) point on the wormhole curve. Based on the default curve, the optimum point is shifted with the amount indicated. The optimum point cannot be lower than (<NUM>, <NUM>). Values in the table are only indicative and serve to illustrate a concept.

Acid reactivity increases with temperature <NUM>, which means that the optimum velocity <NUM> for wormhole growth also increases. For low-permeability reservoirs, it can be difficult to reach the optimum velocity without fracturing the formation. Therefore, it is important to evaluate the downhole temperature of the acid <NUM> when it reaches the formation <NUM>.

As shown in <FIG>, which illustrates the relation between the temperature of the acid at the entrance of the liner <NUM> under different pump rates <NUM> and wellhead temperatures <NUM>. It is an advantage to inject at high rate and at the lowest possible wellhead temperature to limit the in-situ acid reactivity. This is shown by line 304A. As the temperature increases, lines 304B, and 304C we can see an increase in the acid reactivity at the entrance of the liner <NUM>. Furthermore, any brine used to clean out the mud prior to the acid stimulation should be injected at the lowest possible temperature.

The temperature to be used for adjusting the wormhole curve is the temperature of the acid when it enters the reservoir, not the reservoir temperature.

Economides et al. (<NUM>) derived a formula to determine the volume of acid required to achieve a certain wormhole length <NUM> based on the pore volume to breakthrough <NUM> from core flood data:
<MAT>.

The formula is plotted in <FIG>. The ratio V/L is known as the acid coverage in bbl/ft <NUM>. The equivalent skin <NUM> is given as:
<MAT>.

The algorithm aims to achieve a given final skin factor and then calculates the equivalent wormhole radius and then the required acid coverage. However, for economic considerations, the maximum acid coverage is limited by the acid volume which can be pumped. For instance, in offshore wells, the volume is limited by acid boat capacity. In the current application, the acid coverage should not exceed <NUM> bbl/ft.

Alternatively, the acid stimulation can be fixed, which enables calculation of the maximum, final wormhole length <NUM> and consequently the final, negative skin <NUM>.

<FIG> illustrates the outcome of a larger sensitivity analysis involving the pump rate <NUM>, the dP across the last hole <NUM> (502A to 502E, respectively), the discharge coefficient <NUM> (504A to 504E, respectively) and the total hole cross-sectional hole area <NUM>. The linear relationship is very clear and it is therefore possible to predict the dP <NUM> required to obtain a certain cross-sectional hole area <NUM>. This constraint that the cross-sectional LEL hole area <NUM> must be equal to or larger than a minimum cross-sectional area to avoid imposing an addition pressure drop during production/injection after stimulation therefore results in a constraint on the dP across the last hole, which can be estimated based on the relationship provided by the sensitivity analysis. This is a novel concept.

At this stage, we are able to estimate the initial hole-size distribution for the starting point of the algorithm. This is depicted by block <NUM> in <FIG>.

The next step, block <NUM>, requires that the equations are set up. These are then solved as part of the following step, block <NUM>, dealt with later in this specification.

The equation of motion for isothermal one-dimensional pipe flow describes the pressure drop as a contribution from friction, gravity, and acceleration. The gravity term dominates in the vertical section of the wellbore, whereas friction losses become relatively more important in the horizontal section. The acceleration term is only required when velocity changes occur, such as when fluid enters the liner from the tubing (change in inner diameter), or whenever fluid exits through a hole in the liner. The contribution of the acceleration term to the total pressure drop is less than <NUM>% and can often be neglected. <MAT><MAT>.

□ is the angle relative to the z-axis and D is the pipe diameter. The acceleration term can be expressed in terms of volumetric flow Q instead of velocity v,
<MAT>.

The Fanning friction factor, f, is defined in terms of the wall shear stress
<MAT>.

Hence the friction pressure drop for Newtonian flow is:
<MAT>.

Friction pressure loss, in psi/ft, for pipe flow of Newtonian fluids becomes:
<MAT>.

For laminar flow the Fanning friction factor is linked to the Reynolds number,
<MAT>.

Friction pressure loss for annulus flow of Newtonian fluids in concentric pipes:
<MAT>.

For laminar flow in the annulus, the Fanning friction factor is defined as:
<MAT>.

The pressure difference due to the static head is found from:
<MAT>.

Annular flow can occur if acid is first injected through a coiled tubing inside the production string followed by the main acid treatment injected into the production tubing while the coil is still in hole.

The Fanning friction factor for pipe flow in smooth pipes is described by the Prandtl-Karman equation:
<MAT>.

For rough pipes, the friction factor depends on the relative pipe roughness, □/D, and is given as
<MAT>.

<FIG> illustrates the effect of Reynold's number <NUM> on friction factor <NUM> for different values of pipe roughness <NUM> (604A to <NUM>, respectively). A typical relative roughness for a new pipe is <NUM>-<NUM>.

Similarly, for annular flow in smooth pipes, the friction factor becomes
<MAT>.

Finally, for annular flow in roughness pipes, the friction factor is given as
<MAT>.

There is a potential discontinuity going from laminar to turbulent flow because the flow regime is poorly defined in the <NUM>-<NUM> Reynolds number region. This has no impact on the LEL hole design. <FIG> shows that roughness plays a role only if it exceeds <NUM>.

Typical pumping rates are <NUM>-<NUM> bbl/min, depending on reservoir permeability and liner length. Such rates may lead to high surface pressures and thus require the upper completion to be designed appropriately. There is often a need to reduce the friction pressure loss to stay within safe operating limits and this requirement may necessitate the use of drag reducing agents (DRA). Drag reducers are mostly dilute polymer solutions, which lower the frictional resistance to flow in the turbulent regime when added to a solvent, for instance water or acid. Very low concentrations (a few thousand ppm) may in some instances reduce friction by as much as <NUM>%. Friction reducers, may, however, cause reservoir damage, according to some studies.

When adding drag reducing agents a zone named the elastic sub-layer is formed between the viscous sub-layer and the Newtonian core. The extent of the elastic sub-layer will be governed by the amount and type of polymer, and by the flow rate. Within the elastic sub-layer the turbulence structure is significantly different from the Newtonian plug. The turbulent eddies are broken up, whereby the flow characteristics approach those of laminar flow. Therefore, it follows that the friction factor is independent of the relative roughness of the pipe in the presence of drag reducing agents.

Maximum drag reduction is achieved when the elastic sub-layer extends to occupy the entire pipe cross-section. Drag reduction by dilute polymer solutions in turbulent pipe flow is bounded between the two universal asymptotes described by the Prandtl-Karman law for Newtonian turbulent flow and a maximum drag reduction asymptote. In between is the so-called polymeric regime in which the friction factor relations are approximately linear in Prandtl-Karman coordinates, see <FIG>. The polymeric regime may be described by two parameters: The wall shear stress at the onset of the drag reduction, □w*, (or equivalent the onset wave number w*) and the slope increment,δ, by which the polymer solution slope exceeds Newtonian slope. The onset of drag reduction occurs at a well defined onset wave number. For a given polymer solution w* is essentially the same for different pipe diameters. For solutions of a given polymer-solvent combination w* is approximately independent of polymer concentration.

When modelling the effect of the drag reducer, it is assumed that the fluid friction factor is reduced and that the fluid viscosity remains the same. Acid viscosity, via the Reynolds number, has a minor impact on friction losses at typical operating conditions, as seen from <FIG>.

The following formula, developed by Virk (<NUM>, <NUM>) relates the friction factor to the concentration of the drag reducer for pipe flow:<MAT>.

The drag reduction model parameters are
<MAT>.

K and α are constants. The parameters are specific to the chemical used and must be fitted based on flow loop test data provided by the vendor.

The maximum drag reduction asymptote for pipe flow is described by:
<MAT>.

<FIG> shows the impact of drag reduction on the friction factor, in a Prandtl-Karman plot.

For the particular Drag Reducing Agent (DRA) <NUM> model constants used, the maximum asymptote <NUM> is only reached if the DRA concentration exceeds <NUM> ppm. Without the addition of a DRA is shown by <NUM>. Incrementally increasing the amount of DRA is shown by lines <NUM>, <NUM> and <NUM> respectively.

In a <NUM>" ID liner, a pumping rate of <NUM> bbl/min, equivalent of <NUM> bbl/d, leads to a Reynolds number of approximately <NUM>, which is well inside the turbulent flow regime.

<FIG> illustrates the influence of drag reduction on friction pressure <NUM> in a <NUM> ft long <NUM>" OD top completion string as a function of pump rate <NUM>. Friction is reduced to <NUM>/<NUM> by adding <NUM> ppm DRA. The concentration of DRA is similar to that as illustrated in <FIG>.

The limited-entry liner consists of a number of holes allowing fluid to exit the liner and enter the annulus and subsequently the reservoir. The holes are small compared to the liner dimensions, both in terms of length and diameter and can therefore be considered as an orifice. The pressure drop across N holes in the liner may be calculated as:
<MAT>.

Qhole is the flow rate in bbl/min through the holes. The positive direction is from the liner and into the annulus. Dhole is the inner diameter, in inches, of the holes in the liner. N is the total number of holes. CD is the dimensionless discharge coefficient, which accounts for the fact that the pressure loss is only partially recovered due to the short length of the hole (equal to the pipe thickness). Based on the work by Crump and Conway.

(<NUM>), a lower value of <NUM> is used for flow of water and gelled fluids in round sharp-edged drilled holes; values up to <NUM> are also possible, depending on fluid type and how the hole was actually drilled, see El-Rabba et al. (<NUM>) and McLemore et al. CD should be considered a sensitivity variable during the first LEL design jobs. Drilling the holes at a slight angle may reduce the splash-back of unspent acid hitting the formation and improve the jetting process.

The model for the friction factor in the presence of a drag reducer is combined with the model for the friction factor for Newtonian turbulent pipe flow in rough pipes.

If no drag reducers are used then δ = <NUM>. If drag reducers are used the roughness is set to zero.

Inserting the expression for Reynolds number:
<MAT>.

The hole distribution function, <MAT>, is defined as number of holes per foot along the liner. For the present application, it is a step function reflecting different geometrical distribution of holes in sections along the liner. It can also be specified on the basis of the actual liner tally in which case the length of the individual sections equals the length of the joints.

A mass balance across a unit liner section states that the rate lost across a liner segment equals the rate through the holes along the same segment:
<MAT>.

For steady state applications the annulus pressure is a user-provided input, which should not exceed the fracture propagation pressure.

The model for the liner without an inner string includes two differential equations to describe the pressure and flow rate profile within the liner. In addition, a set of algebraic equations defines the Fanning friction factor for Newtonian turbulent flow with and without drag reducer. The inner diameter, the hole diameter and the well inclination may change along the liner.

The flow rate change inside the liner, which is equal to the flow through the holes, is described by:
<MAT>.

It is clear that the liner pressure must exceed the annulus pressure at all times to avoid unwanted cross-flow from a segment with higher reservoir pressure. Zones with substantially higher or lower pressures should be isolated with packers to improve acid coverage control. Transformation into the dimensionless independent variable u=x/L, where L is the total length of the liner and introducing position-dependent constants gives:
<MAT>
<MAT>
<MAT>
<MAT>
<MAT>
<MAT>
<MAT>
<MAT>
<MAT>
<MAT>.

The stimulation rates in bbl/min/ft are found from:
<MAT>.

The flow rate through the holes in bbl/min/hole is found from:
<MAT>.

The algorithm enters into the inner loop <NUM>. This lead by block <NUM>, Solve equations.

The model describes a first order non-linear boundary value problem and consists of two first-order coupled non-linear differential equations with boundary conditions at each end of the interval of the independent variable, coupled with a set of algebraic equations. The problem can be solved in at least two ways, either using a finite-difference scheme, see Mogensen and Hansen (<NUM>), or with polynomial approximation. Following the approach outlined by Hansen and Nederveen (<NUM>), polynomial approximation is used to solve the numerical model.

The liner is spatially discretised with N nodes, as will be described in the subsequent paragraphs. The solution vector, Y, is defined as an array containing the pressure, P, in each of the N nodes, the flow rate, Q, in each of the N nodes and the friction factor, f, in each of the N nodes.

The iteration scheme to be outlined in the subsequent paragraphs is somewhat sensitive to the guess of the initial solution vector. A safe starting guess is to assume that the pressure in the liner decreases linearly with position, while remaining larger than the specified annulus pressure. The same approach is taken for the rate inside each liner segment, which must equal the total flowrate at the inlet and zero at the end. The initial estimate of the friction factor should be <NUM> or lower. Note that since the flow rate is zero at the end of the liner, the flow regime will go from being turbulent around the heel towards laminar at the toe.

The N collocation points are found as zeros of N'th order orthogonal polynomials <NUM>, against the length of the liner <NUM>; in total there are N-<NUM> interior points and <NUM> interval endpoints. Location of the zeros is shown in <FIG> using <NUM> collocation points.

Use of polynomial approximation of order N gives:
<MAT>.

Construction of Lagrange polynomials are defined by:
<MAT>.

The nodes of the polynomial are found from the orthogonality condition:
<MAT>.

Expansion of the four dependent variables into an N'th degree Legendre interpolation polynomial (α = β = <NUM>) gives:
<MAT>
<MAT>
<MAT>
<MAT>.

Based on experience, the number of collocation points is set to <NUM> and is independent of the number of physical holes in the liner. Rewriting the equations gives the discretised model: <MAT><MAT><MAT><MAT>.

The maximum drag reduction asymptote is not included in the model formulation. In practice, use of drag reduction concentrations in the order of <NUM>-<NUM> ppm does not exceed the asymptotic line. <MAT><MAT><MAT> <MAT><MAT><MAT><MAT>.

Q0 is the specified pump rate. The stimulation rates in bbl/min/ft are found from:
<MAT>.

The next step as part of the inner loop <NUM>, is to update the solution vectors, as shown by block <NUM>.

The 3N set of non-linear algebraic equations is solved using a Newton-Raphson iteration procedure. Given an initial estimate of the solution vector an improved estimate is obtained from:
<MAT>.

The Jacobian matrix is defined with the following elements:
<MAT>.

The elements of the Jacobian matrix are set up as described below. The collocation matrix is defined as:
<MAT>.

Function F (i = <NUM>. N):
<MAT>
<MAT>
<MAT>.

Function G (j = <NUM>. N-<NUM>):
<MAT>
<MAT>
<MAT>.

Function H (j = <NUM>. N-<NUM>, i = j):
<MAT>
<MAT>
<MAT>
<MAT>.

Boundary Conditions:
<MAT>
<MAT>
<MAT>
<MAT>
<MAT>.

The final step of the inner loop <NUM> is to determine whether the solution vector is constant, block <NUM>.

Typically, the Newton-Raphson technique will converge within <NUM> iterations using carefully selected relaxation parameters to guide the convergence during the first iterations. The method ensures that final convergence speed is quadratic.

The iterative inner loop will repeat by following arrow <NUM>, and restating resolving the equations, as set forth from block <NUM>.

This iterative inner loop <NUM> finishes when the absolute change to the solution vector is below a certain threshold value, typically 1E-<NUM>. To avoid the possibility of an infinite loop, the procedure stops after a pre-specified number of iterations has been reached, typically in the range <NUM>-<NUM>.

Once the solution vector is deemed constant, the next step is to follow arrow <NUM> to calculate the acid coverage, depicted by block <NUM>.

Once the stimulation flow rates are calculated from the solution, the acid coverage per liner segment is the product of segment flow rate and pumping time. If the overall pump rate changes during the job, the stimulation rate for each segment changes.

The transient period where the acid front moves through along the liner while displacing the brine must also be taken into account. However, this is compensated for when water displaces acid at the end of the job. The time it takes for the front to reach a given position i, is called the retention time, which is calculated recursively:
<MAT>.

Since the liner flow rate gradually decreases towards zero at the heel, it is clear that it takes gradually longer time for the acid front to displace the brine out of the liner. In other words, the inner part sees acid for longer time than the outer part. The hole-size distribution should compensate for this. The retention time is therefore also a measure of the minimum time needed for water to displace acid from the liner at the end of the stimulation.

The next step is shown in block <NUM>, to determine if the dp across the last hole is matched. This step goes in combination with the following block <NUM> which is to determine if design acid coverage matched?.

The dP across the last hole is calculated as the difference between the pressure in the last node of the liner and the annulus stimulation pressure (which is constant and user-specified):
<MAT>.

The difference between calculated and specified target acid coverage is given as
<MAT>.

This formulation ensures that the dCOV (acid coverage distribution) function is always positive. Hence, it must be minimized to obtain the best possible match. The relative acid coverage is determined as follows:
<MAT>.

Turning to <FIG>, the above two steps are combined into block <NUM>.

Whereas the inner loop <NUM> consists of solving the material balance for a given combination of LEL holes, pump rate and other variables, the first part of the outer loop <NUM> consists of adjusting the LEL hole size distribution to match both the desired pressure drop across the last hole, dP <NUM>, and the desired acid coverage for each segment <NUM>. The outer loop <NUM> serves to honour both constraints at the same time.

Therefore, the hole size-distribution must be satisfied, as shown in block <NUM>.

If the dP is too small <NUM>, then there are too many LEL holes and one LEL hole is then subtracted from the segment with the highest relative acid coverage <NUM> and the material balance inner loop <NUM>, via block <NUM>, is then reinvoked.

If the dP is too large (arrow <NUM>), there are too few holes, and one hole is then added to the segment with the lowest, non-zero relative acid coverage <NUM> and the material balance inner loop <NUM>, via block <NUM>, is reinvoked.

Segments with zero acid coverage are not adjusted.

If the dP is close to the target value within a certain tolerance, then the acid coverage distribution dCOV is calculated <NUM>. At this point, the total number of LEL holes is correct but the holes just need to be redistributed among segments. One LEL hole is added to the segment with the lowest, non-zero relative acid coverage, whereas one LEL hole is subtracted from the segment with the highest relative acid coverage <NUM>. Then the inner loop <NUM>, via block <NUM>, is reinvoked and the procedure is repeated until the dCOV function reaches a minimum. Since the algorithm adjusts integer values, i.e. number of LEL holes, it is not possible for the dCOV function to be exactly zero.

Once the minimum dCOV function is reached, it must be determined if the calculated wellhead pressure (WHP) is below the wellhead pressure maximum constraint, as shown in block <NUM>.

Every wellhead has a maximum pressure rating, such as <NUM> psia, <NUM> psia and higher. Similarly, every tubing has a maximum pressure rating. Therefore, if the reservoir pressure is high, the design rate may give rise to a wellhead pressure, which exceeds the pressure rating.

If the calculated wellhead pressure exceeds the maximum rating of the tubing (shown by arrow <NUM>), adjustment to the design (block <NUM>) requires the following steps:.

Next, the average hole distance constraint must be met, block <NUM>.

As described earlier, Economides et al. (<NUM>) derived a formula to determine the volume of acid required to achieve a certain wormhole length based on the pore volume to breakthrough from core flood data:
<MAT>.

The ratio V/L is known as the acid coverage in bbl/ft. The equivalent skin is given as
<MAT>.

Schwalbert et al. (<NUM>) defined the stimulation coverage as twice the wormhole radius relative to the length of the perforated interval, which for LEL completions equals the distance between LEL holes.

Turing to <FIG> which illustrates skin factor as a function of stimulation coverage.

Therefore the skin factor <NUM> becomes constant when the stimulation coverage <NUM> reaches <NUM>%.

<FIG> shows that an effective wormhole radius <NUM> of <NUM> ft would result in an equivalent negative skin factor <NUM> of -<NUM>, assuming that all the wormholes generated along the well have the same radius. Combining the two plots shows that the maximum distance between wormholes should not exceed twice the wormhole length. A skin of -<NUM> for the entire well, for instance, means that the holes should be drilled with a maximum distance of <NUM> ft.

This means that the average distance between LEL holes, defined as the length of the stimulate reservoir section divided by total number of holes, should not exceed twice the expected final wormhole radius. The following check is therefore performed:
<MAT>.

If the average hole constraint is not met it becomes necessary to adjust the hole size, block <NUM>.

Based on the evaluation of the above equation, the following possible actions are taken:.

If the distance between LEL holes is too small, the LEL hole size can be increased by <NUM> and the entire simulation is then repeated.

If the distance between LEL holes is too large, the LEL hole size can be decreased by <NUM> and the entire simulation is then repeated.

If the average distance between LEL holes is close to or equal to twice the wormhole radius, the algorithm has converged with a final design and proceeds with output of results, block <NUM>.

Node properties, including position, pressure, rate, friction factor, number of holes per foot, velocity, retention time, stimulation rate, cumulative volume of acid leaving the node through holes.

Segment properties, including segment number, segment interval, number of holes in segment, distance between holes, calculated and design acid coverage, acid coverage ratio, acid stimulation rate, acid velocity at the exit point of the holes, pore volume to breakthrough, final wormhole radius, and final skin factor.

Total liner volume, total tubing volume, displacement volume, retention time
A detailed tally list containing the number and size of LEL holes for each joint to be run in hole, as well as the order in which the joints must be run in hole. Furthermore, the total number of joints with a particular number and size of holes is summarised, such as number of joints with <NUM>, <NUM>, <NUM> or <NUM> LEL holes of size <NUM>, <NUM>, or <NUM> etc..

Wellhead pressure and bottom-hole pressure during pumping are calculated from the pressure at the first node and then subtracting hydrostatic pressure and adding friction up to the given gauge depth.

Wellhead and bottom-hole instantaneous shut-in pressures ISIP are calculated from the pressure at the first node and then subtracting hydrostatic pressure up to the given gauge depth. The friction is zero because the rate is zero during an ISIP.

Input to the numerical design model includes:.

The software will then estimate the design pump rate based on the standard transient inflow model (not the Darcy model, which is a steady-state assumption) while ensuring that the injection pressure remains below fracturing pressure. Key parameters include the permeability, the length of the completed interval, and the difference between annulus pressure and reservoir pressure. For the calculation, it is assumed that the skin can be reduced to zero. Thus, note that because the stimulation job typically takes less than <NUM> hours, the injectivity is higher than predicted by the Darcy formulation. The reason is that the boundaries are not yet felt by the pressure pulse emitted during stimulation. So, even though the flow inside the liner is a steady-state formulation, the inflow model used for design of the pump rate is transient.

A nodal analysis calculation must be performed to estimate the downhole temperature at the heel of the liner. Based on the choice of acid system, the permeability, and the temperature, the optimum velocity for wormhole propagation is estimated, together with the anticipated pore volume to breakthrough based on published literature data. The Buij se-Glasbergen model is used to characterise the wormholing at different velocities.

Based on the optimum velocity and the calculated design pump rate, it is straightforward to calculate the total cross-sectional area of the holes. This area is linearly correlated with the dP across the last hole, which is a key design parameter.

The stimulation design aims for a negative skin of -<NUM> or better, which requires the holes to be not more than <NUM>-<NUM> ft apart. The model by Economides et al. (<NUM>) is used to calculate the acid coverage required to achieve this skin. A higher acid coverage requires more acid and longer pumping time and hence higher cost. This must be weighed against aiming for a more negative skin.

Provide an initial estimate of the number of holes per segment and let the software find the solution which honours the acid coverage per segment and the dP across the last hole. The initial estimate can be found from the relationship between interstitial velocity, pump rate, and total cross-sectional hole area for a particular discharge coefficient and liner configuration.

To illustrate the design concept in more detail, an example is shown below. The well in question will have an approximate reservoir length of some <NUM> ft. Since a stand (<NUM> drill pipe lengths) is approximately <NUM> ft, the well is numerically split into <NUM> segments, each with a length of <NUM> ft, corresponding to <NUM> stands.

The initial design coverage is set to <NUM> bbl/ft. The transient inflow equation predicts that the maximum rate without fracturing the formation is <NUM> bpm, assuming that the skin is zero. As the stimulation progresses, the rate can be increased further. The resulting pumping time will be <NUM> hours, which leads to a slight adjustment of the design rate, but not much.

Although the reservoir temperature is <NUM> F or more, nodal analysis based on the design rate of <NUM> bpm predicts a BHT of <NUM> F at the first hole. This temperature is used for estimating the position of the optimum velocity for wormhole propagation based on a measured curve and the Buij se-Glasbergen model.

The final skin is initially assumed to be -<NUM>, which yields a maximum distance between adjacent holes of <NUM> ft. This corresponds to a pressure drop across the last hole of about <NUM> psia, which is then used as input for the design model.

The discharge coefficient is assumed to be <NUM>, which is mid-way between the theoretical minimum of <NUM> and a high value of <NUM>. Post-job analysis will help identify the pressure drop across the holes and hence the actual discharge coefficient.

A first estimate of the hole size distribution makes use of a linear relationship between hole cross-section area and pressure drop across the last hole. Based on this initial input, the actual optimum hole size distribution is calculated using the numerical algorithm outlined. In the inner loop, the flow equations are solved. In the outer loop, the number of holes is adjusted to match the pressure drop across the last hole as well as the acid coverage for each segment.

The results from the calculations are show in the four plots above. The distance between adjacent holes is in the range <NUM>-<NUM> ft, which yields optimum stimulation coverage (wormholes cover the entire well length). The distance is not uniform because the hole size is chosen to be constant at <NUM> to avoid complicating the pilot design.

Based on the wormhole growth model of PVbt versus interstitial rate, the minimum PVbt to be inserted into the skin model by Economides and based on the specified acid coverage of <NUM> bbl/ft. This yields a skin factor of -<NUM>, which is considered close enough to the initial estimate of -<NUM>. If a skin of -<NUM> is desired, we would need to increase the acid coverage, recalculate the pumping time, recalculate the flow rate, redesign the hole sizes and then check the resulting skin.

While the embodiment of the invention has been described above and discussed in detail, the invention is not deemed to be restricted to this particular embodiment. A person skilled in the art will appreciate that a number of variations may be made to the described embodiment or features thereof, without departing from the scope of the present invention.

Claim 1:
A method, performed by a data-processing system, of simulating fluid transport in a system for stimulating a well in a material formation (<NUM>), which system comprises a limited entry liner, LEL, (<NUM>), wherein the liner (<NUM>) is divided into segments (<NUM>), each of which having a length less than that of the total length of the liner, and which includes one or more holes (<NUM>) along a wall of the liner for discharging a fluid into the formation, characterized in that
the method provides (<NUM>) an initial estimate of the number of holes along the wall of the liner which honours the acid coverage per segment and the drop in pressure across the last hole, for an optimized hole-size distribution in construction and operation of the system,
adjusting (<NUM>) the initial estimate of the number of holes along the wall of the liner if the acid coverage per segment and the drop in pressure across the last hole is not honoured; and wherein the method further comprises:
performing a series of algebraic equations for an initial hole-size distribution guess;
calculating acid coverage and dp across the last hole;
comparing acid coverage and dp across the last hole against a design variable in a first iteration;
evenly decreasing the number of holes across a segment for a next iteration until dp across the last hole is honoured; or
evenly increasing the number of holes across a segment for a next iteration until dp across the last hole is honoured, as a first step; and
performing a second step which includes:
redistributing existing number of holes between various segments as a first iteration, wherein segments, where the calculated acid coverage is the furthest away from design values, exchange one hole;
performing the next iteration until acid coverage is honoured; and performing the first step and the second step until dp across the last hole and acid coverage is honoured; and
wherein the method further comprises:
operating the system for stimulating the well in the material formation (<NUM>) using the optimized hole-size distribution.