CN108779669A - The continuous fully implicit solution well model with tridiagonal matrix structure for reservoir simulation - Google Patents

Abstract

The subterranean hydrocarbon reservoir with horizontal well or multiple vertical wells is simulated by the continuous solution of reservoir and well equation, to simulate the flowing of the fluid in reservoir and well productivity.Well processing is defined bottom hole pressure well by the continuous solution of reservoir and well equation.This avoids the big matrix caused by the simultaneous solution of reservoir and well equation is solved, it can be computationally very expensive for a large amount of unknown number and need dedicated sparse matrix solver.This continuous solution is related to the reservoir systems solver for the rule supplemented by the minor matrix of the numerical solution of bottom pore pressure force.To the tridiagonal matrix of adjacent with well unit reservoir units at perforated interval；And the vector of the unknown reservoir gesture of adjacent reservoir units executes the solution.

Description

The continuous fully implicit solution well with tridiagonal matrix structure for reservoir simulation Model

Cross-reference to related applications

The application is the pending jointly owned U.S. Patent application No.14/ of the applicant submitted on the 30th of September in 2013 040,930 part, which continues, applies and requires its priority, and requires on 2 9th, 2011 U.S. Patent applications submitted The priority of No.13/023,728 (being the United States Patent (USP) No.9,164,191 of publication on October 20th, 2015 now).

The application further relates to entitled " Sequential Fully Implicit Well Model with Tridiagonal Matrix Structure for Reservoir Simulation, " (act on behalf of Reference Number No.004159.005477) and with inventor same as the present application the jointly owned United States Patent (USP) Shen submitted on the same day It please No.15/061,572.

Technical field

This application involves the computerization of the hydrocarbon reservoir in earth simulations, and in particular it relates to along reservoir In the flow profile along each well simulation.

Background technology

Well model plays an important role in numerical value reservoir simulation.Well model be used for calculating from oil and The oil of well in natural gas reservoirs, water and natural gas productivity.If well productivity is known, they are used for calculating edge The flow profile of the perforated interval of well.With the increase of the ability of the flow velocity for measuring the perforated interval along well, number appropriate Value well model is necessary for calculating correct flow profile with the measurement matched in reservoir simulator.

Commonly known to be, if for any vertical well in reservoir simulator, all reservoirs vertically communicate, then simply Well model (e.g., explicit or half implicit model) may be enough.For these models, well productivity and the layer production index (or Overall flow rate) proportional distribute to perforation.Therefore, it calculates simple and computationally cheap.Reservoir unknown number it is obtained The structure of coefficient matrix remains unchanged.Specifically, coefficient matrix maintains regular sparsity structure.Therefore, any this sparse matrix Solver can be used for each time step for entire reservoir simulation model, solve the linear of gridblock pressure and saturation degree System.

However, for some vertical non-highly non-uniform reservoirs for communicating layer, well model above-mentioned cannot Generate correct physical solution.However, they generate incorrect flow profile, and in some cases, simulator is caused to be restrained Problem.

Increase with the complexity of reservoir model, the quantity of perpendicular layers is to hundreds of orders of magnitude to indicate reservoir unevenness Even property.The fully implicit solution unity couping well model of simultaneous solution with reservoir and well equation cuts open correct simulation along the flow of well Face has been necessary, and is necessary also for the numerical stability of reservoir simulation.It is usually first in order to solve unity couping system First eliminate well equation.This generates the unstructured coefficient matrixes of reservoir unknown number to be solved.The solution of such matrix Need the special solver with special preconditioner.For many wells in well and simulation model with many completions Speech, it is expensive that this method, which became to calculate in terms of the processor time,.

Invention content

Brief, the present invention provides a kind of new and improved the method implemented by computer, utilizes solution reservoir simulation Well equation in model and reservoir equation, reservoir simulation model have the stratum for carrying vertical fluid flowing and without vertical The flow obstacle layer of fluid flowing.Pass through the flowing that will there is flowing only in the reservoir model near well between fluid barrier Layer group becomes single fluidized bed, and the method implemented by computer is formed in the horizontal and multiple vertical well in reservoir simulation model Simplify well model system.For specifying the well of productivity, this method then (to simplify well model system for the grid block that well passes through System) well and the bottom hole pressure of reservoir unknown number (use direct sparse matrix solver) calculated by matrix solve simplification Well model system.This method is recycled and reused for one or more of reservoir simulator well.Then, this method is then by will be every The processing of a well based on layer completion rate (completion rate), strata pressure and is led for determining bottom hole pressure The steady state volumes equilibrium relation of water rate solves full reservoir simulation model, with the layer of one or more wells of the full reservoir model of determination Completion rate, and determine according to the completion rate of the determination of the layer of one or more wells total rate of well.This method and then forming layer Determination completion rate and one or more wells determination total rate record.By using distributing to simulator grid block The perforation rate of reservoir data and calculating, simulator calculate the pressure distribution of the Multiphase Flow fluid saturation distribution in reservoir. Therefore, simulator can calculate the pressure in reservoir in the case of the production of the Jing Chu under given time and given spraying rate And saturation distribution.

The present invention provides a kind of new and improved data processing system, is used to store up by the underground from reservoir model The reservoir simulation of layer forms the well model of horizontal and multiple vertical wells, reservoir model have stratum with fluid flowing and The flow obstacle layer flowed without fluid.Data processing system includes processor, executes following steps：Pass through matrix meter It calculates, solves and simplify well model system model, obtain the storage around the bottom hole pressure and one or more wells of one or more wells Layer unknown number, and by will (one or more) well processing be with determining bottom hole pressure, and be based on layer completion rate, The steady state volumes equilibrium relation of stressor layer and hydraulic conductivity solves full reservoir model, with (one or more of the full reservoir model of determination It is a) the completion rate of each layer of well.Processor determines (one also according to the completion rate of the determination of each layer of (one or more) well Or multiple) total rate of well.Data processing system further includes memory, one for forming the completion rate of the determination of each layer and determining The record of a or multiple total rates.

The present invention also provides a kind of new and improved data storage devices, store the operable instruction of computer In computer-readable medium, described instruction make data processor by the reservoir simulation of the subsurface reservoir from reservoir model come The well model of horizontal and multiple vertical wells is formed, reservoir simulation has the stratum for carrying fluid flowing and flowed without fluid Flow obstacle layer, to execute following steps：Calculated by matrix, solve simplified well model system model, obtain one or The bottom hole pressure and reservoir unknown number of multiple wells, and by being with determining base apertures by the processing of (one or more) well Pressure simultaneously solves full reservoir model based on the steady state volumes equilibrium relation of layer completion rate, strata pressure and hydraulic conductivity, with determination The completion rate of the layer of (one or more) well of full reservoir model.The instruction being stored in data storage device further includes instruction, It makes data processor and is determined according to the completion rate of the determination of the layer of one or more wells total rate of (one or more) well, And form the record of the completion rate of the determination of each layer and total rate of the determination of (one or more) well.

Description of the drawings

Figure 1A and Figure 1B is to be formed as the more of single layer above and below the fluid barrier according to the present invention in reservoir The schematic diagram of well model in the reservoir simulator of a subsurface formations.

Fig. 2A and Fig. 2 B are above and below several fluid barriers being spaced vertically apart from according to the present invention in reservoir Be formed as the schematic diagram of the well model in the reservoir simulator of multiple subsurface formations of single layer.

Fig. 3 is the schematic diagram of the well model for simulating based on explicit model method shown in radial geometry.

Fig. 4 is the schematic diagram of the well model for simulating based on fully implicit solution fully-coupled model method.

Fig. 5 is the signal of the well model of the reservoir simulation of the comparison with the flow profile obtained from the model of Fig. 3 and Fig. 4 Figure.

Fig. 6 is the schematic diagram of the Finite Difference Meshes system of the well model for Fig. 3 and Fig. 4.

Fig. 7 A and Fig. 7 B are the well mould for showing unmodified conventional well layer model and well layer model according to the present invention respectively The schematic diagram of type reservoir.

Fig. 8 A and Fig. 8 B are the schematic diagram for the flow profile for showing the comparison between the model of Fig. 7 A and Fig. 7 B respectively.

Fig. 9 is that there are two the schematic diagrames of the well layer model of zone of fracture for tool in radial coordinate.

Figure 10 is the data of the method and system of the fully implicit solution unity couping well model of the simultaneous numerical solution with reservoir simulation The functional block diagram or flow chart of processing step.

Figure 11 is at the data of the method and system of the continuous fully implicit solution well model according to the present invention for reservoir simulation Manage the functional block diagram or flow chart of step.

Figure 12 is the explicit well model for the one-dimensional reservoir simulator with vertical well with three diagonal coefficient squares The schematic diagram of the linear system of the equation of battle array.

Figure 13 is the fully implicit solution well model for the one-dimensional reservoir simulator with vertical well with three diagonal coefficients The schematic diagram of the linear system of the equation of matrix.

Figure 14 is for the constant bottom hole pressure according to the present invention for the one-dimensional reservoir simulator with vertical well Schematic diagram of the linear system of the equation with three diagonal coefficient matrixes of well model to be handled.

The functional block diagram or flow chart for the step of Figure 15 is the analysis method according to the present invention for showing reservoir simulation.

Figure 16 is the schematic diagram of the computer network of the continuous fully implicit solution well model of reservoir simulation according to the present invention.

Figure 17 is the schematic diagram of three diagonal coefficient matrixes.

Figure 18 is the schematic diagram of the linear system of the equation with three diagonal coefficient matrixes.

Figure 19 is the signal of the Finite Difference Meshes system of the Model of Horizontal Well according to the present invention being orientated in the y-axis direction Figure.

Figure 20 A and Figure 20 B be it is according to the present invention it is before and after being formed as simplified Model of Horizontal Well, in reservoir In flow obstacle region in multiple subsurface formations reservoir simulator in Model of Horizontal Well schematic diagram.

Figure 20 C are the schematic diagrames of the linear system of the equation of the Model of Horizontal Well of the simplification of Figure 20 B.

Figure 21 and Figure 22 is the data processing step of the method and system of the Model of Horizontal Well of simplification according to the present invention Functional block diagram or flow chart.

Figure 23 is the schematic diagram of the Finite Difference Meshes system of the reservoir model according to the present invention with multiple wells.

Figure 24 A and Figure 24 B are the signals for the symbol for showing the grid name in multiple well reservoir models of Figure 19 and Figure 23 Figure.

Figure 25 A, Figure 25 B, Figure 25 C and Figure 25 D are the schematic diagrames for the grid number for showing three dimensional reservoir framework.

Figure 26 is multiple in the reservoir for the fully implicit solution unity couping well model of the simultaneous numerical solution with reservoir simulation The functional block diagram or flow chart of the data processing step of the method and system of the modification of the block diagram of Figure 11 of well.

Figure 27 be simplified two wells, three dimensional reservoir framework equation linear system schematic diagram.

Specific implementation mode

By way of introduction, the present invention provides the continuous fully implicit solution well model of reservoir simulation.Reservoir simulation is reservoir work The mathematical modeling science of journey.The fluid stream in oil or natural gas reservoirs (porous media) is described by one group of partial differential equation It is dynamic.It is any moment of these equations description during the service life for the reservoir for generating oil, natural gas and water, any in reservoir Point place pressure (energy) distribution, oil, water and natural gas rate distribution, oil, water, natural gas volume fraction (be saturated Degree).Fluid flowing in mobile describing reservoir by tracking the ingredients of a mixture.With unit mass or mole expression composition (e.g., methane, ethane, CO₂, nitrogen, H₂S and water) amount.

Due to these equations and relevant thermodynamics and physical laws of description fluid flowing be it is complicated, they It can only be solved on digital computer to obtain pressure distribution, VELOCITY DISTRIBUTION and the stream at any point in the reservoir at any moment The amount of body saturation degree or composition quality or molar distribution.This can only be completed by numerically rather than analytically solving equation. Numerical solution needs reservoir to be subdivided into calculating elements (unit or net on region and vertical direction (x, y, z --- three dimensions) Lattice block), and the time be subdivided all day long or the interval of the moon.For each element, determined by the math equation of solving complexity unknown Number (pressure, speed, volume fraction etc.).

In fact, reservoir simulator model can be considered as the collection of rectangular prism body (as the brick in building walls) It closes.Since the oil, water and natural gas of the Jing Chu being distributed in reservoir produce, it may occur that the variation of pressure and velocity field.With when Between (t) simulated.In general, during the production history of reservoir, the production of each well or injection rate are known.However, by Several reservoirs (element) are passed through in well, therefore contribution of each formation elements (well perforation) to production is calculated by distinct methods. The present invention relates to the calculating how many contribution the perforation of each well produces total well.Since these calculating are expensive, and right Boundary condition is very important for simulator, the method suggestion proposed is used for being computed correctly the flow along well track The practical methods of section.As will be described, it can show which：The some other methods used will cause incorrect flow to cut open The problem of face, this can cause to obtain correct numerical solution, and these methods can make it is very expensive in calculating.

Name

Δ x, Δ y, the size of mesh opening on the directions Δ z=x, y and z

k_x, k_y, k_zPermeability on the direction=x, y and z

P=pressure

ρ=fluid (oil) density

G=gravitational constants

The vertical depth of z=self-reference plane depth

r_oThe Δ x of=Pi Siman (Peaceman) radius=0.2

r_w=drilling well pore radius

T_x,T_y,T_zRock hydraulic conductivity on the direction=x, y and z

It proposes to describe common reservoir simulation model in equation (1) below and indicate related to the present invention interested The equation of well item：

Wherein, Δ_xIt is the difference operator on the directions x of reservoir, T_x,T_y,T_zIt is the x such as limited in following equation, y With the rock hydraulic conductivity on the directions z), j represents the number of fluid phase, n_pIt is the sum of fluid phase, usually 3 (oil, water and natural Gas), Σ is sum term, ρ_i,jIt is the density of the component i in fluid phase j, λ_jIt is the mobility (equation 6) of phase j, Φ_jIt is fluid phase j Fluid potential (pressure of base-level correction), similar, Δ_yIt is the difference operator on the directions y of reservoir, and Δ_zIt is reservoir Difference operator on the directions z, q_{I, w, k}It is the well item (source or remittance (sink)) of the component i of grid block (unit) k, Δ_tIt is in time domain Difference operator, n_iIt is the total moles and n of component i_cBe the component in fluid system sum (methane, ethane, propane, CO₂Deng).

Equation (1) is the partial differential equations of the Non-linear coupling of the fluid flow in describing reservoir.In above-mentioned equation group In, n_iIndicate i-th of component of fluid.n_cIt is the sum of the component of the hydrocarbon and water flowed in reservoir.Herein, Component means such as methane, ethane, propane, CO₂、H₂S, water etc..The quantity of component depends on the carbon that can be used for interested reservoir Hydrogen compound water system.In general, the quantity of component can change to 10 from 3.Equation (1) is combined with continuity equation and momentum side Journey.

In equation (1), q_{I, w, k}It is the position x of component i_k, y_k, z_kThe well perforation rate and k at place are that grid block is (single Member) number.Furthermore the calculating for carrying out this according to the defined productivity at well head is subject of the present invention.

In addition to the differential equation in equation (1), the pore volume constraint at any point (element) in reservoir must satisfy：

Wherein, V_pIt is grid block pore volume, p (x, y, z) is the Fluid pressure at point x, y, z, N_jIt is in fluid phase j Total moles, ρ_jIt is the density of fluid phase j.

There are n in equation (1) and equation (2)_c+ 1 equation and n_c+ 1 is unknown.For each component i, these equations are logical It crosses equation (3) and mutually constrains simultaneous solution with thermodynamics：

Wherein, f_iIt is component fugacity (component fugacity), subscript V represents steam phase, and L represents liquid phase, n_iIt is The total moles of component i, P are that pressure and T indicate temperature.

In fluid system, usually there are three kinds of fluid phases in reservoir：Oil phase, natural gas phase and water phase.Each stream Body mutually may include the different amounts of said components based on reservoir pressure and temperature.Fluid phase is described by symbol j.Symbol j has Maximum value 3 (oil phase, water phase and natural gas phase).Symbol n_pBeing the largest the number of phases (sometimes, can be 1 (oil)；2 (oil With natural gas or oil and water)；Or 3 (oil, water and natural gas)).Number of phases n_pBecome based on reservoir pressure (p) and temperature (T) Change.Symbol n_iIt is the total moles of the component i in fluid system.Symbol n_cIt is the maximum number of the component in fluid system.Root The number of phases and each phase n are determined according to equation (3)_{I, j}Each of point score and phase density ρ_jAnd ρ_{I, j}.In equation (3), V Steam (natural gas) phase is represented, L represents liquid phase (oil or water).

The total moles in fluid phase j are defined by following formula：

Total component molar is defined by equation (5).

Phase mobility in equation (1), the relationship between phase, fluid potential and differential sign are defined on equation (6) to (9) In be defined.

λ_j=k_{R, j}/μ_j (6)

In equation (6), the molecular definition relative permeability of phase, and denominator is the viscosity of phase.

The capillary pressure between the phase relative to phase pressure is defined by equation (7)：

The fluid potential of phase j is defined by following formula：

Φ_j=P_j-gρ_jz (8)

X is defined by following formula, the discrete differential operator on the directions y and z：

Wherein, Δ defines discrete differential symbol and U is aleatory variable.

For each grid block (unit) (grid more in the block may include well) in reservoir simulator, control is used Constraint and definition in equation (1) and equation (3) to equation (9) is written in finite volume difference method.Obtained by simultaneous solution Equation.The given well productivity q of each well is found after the completion_TN_i(x, y, x and t), P (x, y, z and distribution t), according to The present invention gives well productivity using new well model formulation, according to this come the component rate in accounting equation (1).In order to ask (1) and equation (2) are solved equation, by the reservoir boundaries in the space (x, y, z), rock properties are distributed K (x, y, z), and rock is porous Property distribution and fluid properties and saturation degree related data input simulation in.

According to the present invention and as will be described below, simplified well model system is formed, is generated and complicated calculating The determination of the bottom hole pressure of the time-consuming first identical calculating of unity couping well model.

According to the present invention, it has been determined that the grid block for vertically communicating place by several stratum for well track communicates layer It can be combined to which processing is single layer, as schematically indicated for well model in Fig. 1 and Fig. 2.This passes through accomplished below：Know The various flow barrier layer above or belows of perpendicular flow barrier and combination well unit in the reservoir of other well unit.Cause This, full well model system is reduced to the well model system with the reduced size for having lacked many layers, come be incorporated to well model into Row processing.

As shown in Figure 1, well model L indicates that the complicated subsurface reservoir grid block that well passes through is (single with simplified schematic form Member), it is made of seven individual stratum 10, each of stratum 10 is vertically in fluid communication with adjacent layer 10.Model L packets Another group of ten stratum 12 around well 12 are included, each of stratum 12 is vertically in fluid communication with adjacent layer 12.In mould The group on the stratum 10 and 12 that the adjacent layer similar with other in type L is in fluid communication is as shown in Figure 1 by vertical fluid The fluid-tight barrier layer 14 that stream constitutes obstacle separates.

According to the present invention, by the purpose for determining gesture Φ and completion rate, by 14 top of fluid barrier of well model L Layer 10 is combined or is combined into the composite layer 10a in reduced model R, and for processing intent, well model L is reduced to simplified Or the well model R (Figure 1A) simplified.Similar, the layer 12 of 14 lower section of fluid barrier of well model L is combined into simplified well model R In composite layer 12a.

Similar, as Fig. 2 is indicated, reservoir model L-1 is made of the individual stratum in five tops 20, every in stratum 20 It is a to be vertically in fluid communication with adjacent layer 20.Model L-1 includes another group of seven stratum 22, each of stratum 22 with it is adjacent Layer 22 be vertically in fluid communication.The group on the stratum 20 and 22 that adjacent layer similar with other in model L-1 is in fluid communication is such as The fluid-tight barrier layer 24 by constituting obstacle to vertical fluid stream indicated in Fig. 2 separates.Such as indicated in model L-1 , nine stratum 25 of another group be in fluid communication each other below the fluid barrier layer 26 as vertical fluid stream obstacle and layer 22 separate.Nethermost one group of ten stratum 27 being in fluid communication each other are located in model L-1 below fluid stream barrier layer 28.

According to the present invention, by the purpose for determining well layer gesture Φ and completion rate, by the fluid barrier 24 of well model L-1 The layer 20 of top is combined or is combined into the composite layer 20a in reduced model R, and for processing purposes, well model L-1 is simple Turn to well model R-1 (Fig. 2A) that is simplified or simplifying.Similar, by 24,26 and 28 lower section of fluid barrier of well model L-1 Other layers 22,25 and 27 are combined into composite layer 22a, 25a and 27a in simplified well model R-1.

Simplified well model system according to the present invention or well model are solved for reservoir unknown number and bottom pore pressure force. Then, the well in full reservoir simulation model system is considered as defined base apertures well pressure, and it is unknown impliedly to solve reservoir Number.The diagonal element and right end item vector of the coefficient matrix of reservoir model are only points changed in processing according to the present invention Amount, and these are only slightly modified.Then regular sparse solver techniques or methods are used for solving reservoir unknown number.By making With reservoir unknown number (pressure and saturation degree) calculation punch rate.Then it sums these rates to calculate total well rate.Root It will be with simulator for each time step according to the error between total well rate and input well rate of the determination of the present invention Nonlinear Newtonian iteration and reduce.

For the entire reservoir simulation model including many wells, the flow velocity that is calculated according to the present invention is also as passing through unity couping Simultaneous solution calculate flow velocity and restrain.Because the present invention needs to solve The Small Well system model, therefore, it is relatively low to calculate cost.? It was found that if by when combination communicates layer suitably using roughening (upscaling) by be suitably constructed simplified well system, The method convergence of the present invention.

Commonly known to be, if all reservoirs vertically communicate, simple well model (e.g., explicit or half implicit model) is logical Often it has been sufficient.As shown in Figure 3, the quantity that explicit well model E is communicated by vertical current is that the reservoir 30 of Nz forms, every layer With the permeability k for indicating to limit in such as Fig. 3_{X, i}(herein, i expressions level number, rather than component) and thickness deltat z_iAnd perforation Layer rate q_i.Then, in figure 3, as indicated in the equation (3) in same figure, the total output q of explicit model E_TIt is explicit The single productivity q of Nz layer of model_iSum.

For explicit model, well productivity proportionally distributes to perforation with the layer production index (or overall flow rate).Therefore, It calculates simple.The coefficient matrix of obtained unknown number remains unchanged, that is, maintains regular sparsity structure, matrix such as in fig. 12 Shown in format.Therefore, any sparse matrix solver can be used for being directed to each time step, solve gridblock pressure and saturation The linear system of degree.

Well model

The method of several vertical well models of reservoir simulator is also based simply on can slightly compress in reservoir The simple fluid system of the liquid form of single-phase petroleum streams is presented.It will be appreciated, however, that being commonly available to this hair of reservoir It is bright and can be used in any amount of well in regular reservoir simulation model and fluid phase.

Fig. 6 shows the Finite Difference Meshes G for vertical well model used in this specification.As shown, well exists In vertical direction at the center of central location.The model being set forth below is additionally contemplates that the completion on the vertical directions Nz Gesture in well and adjacent cells：Φ_BW,Φ_BE,Φ_BN,Φ_BSFor constant and according to dry run (time step before or repeatedly Generation value) it is known.Here, subscript B refers on " boundary ", and W indicates that the adjacent cells in west side, E represent the adjacent cells in east side, and N Bei Henan is respectively represented with S.Furthermore Φ describes fluid potential (pressure of base-level correction).

As shown in Figure 6, well upwardly penetrates through several reservoirs in the Vertical Square indicated by index I (i=1,2,3 ... Nz), In, Nz is the sum of the perpendicular layers in reservoir model.It is adjacent there are four at same area plane (x, y) for each layer of i Unit.These adjacent cells are located at thing (directions x) He Beinan (directions y).

The gesture Φ of the adjacent cells in east, west, north and southern side is known according to the calculating of simulator time step, and those Gesture is arranged to assume that these gesture do not change with simulator time step.Then consider to be located at the vertical gesture change that well location sets place Change, but assumes that the adjacent gesture that can change in vertical direction is known.

The steady state volumes equilibrium equation of unit (i) in Fig. 6 is as follows：

T_wi(Φ_Bw-Φ_i)+T_Ei(Φ_BE-Φ_i)+T_Ni(Φ_BN-Φ_i)+T_Si(Φ_BS-Φ_i)+ (10)

T_Up,i(Φ_i-1-Φ_i)+T_Down,i(Φ_i+1-Φ_i)-q_i=0.

In equation (10), T indicates the hydraulic conductivity between unit.Subscript W, E, N and S indicate west, east, south and it is northern to, And (i) indicate unit index.

Hydraulic conductivity between the unit in three directions is limited by following equation (11)：

In above-mentioned equation (11), Δ x_iIt is that (unit is big for grid block size on the directions x of grid block (unit) number i It is small).Similar, Δ y_iIt is grid block size (cell size) and the Δ z on the directions y of grid block (unit) number i_iIt is grid Grid block size (cell size) on the directions z of block (unit) number i or grid layer thickness.K_{Z, i-1/2}It is unit i and unit i-1 Intersection vertical permeability.Similar, K_{Z, i+1/2}It is the vertical permeability of the intersection of unit i and unit i+1.Such as Fig. 6 Shown, unit i is at center, intersections of the i+1/2 between unit i and unit i+1.For simplicity in expression thing When stream, increase subscript j.Similar, (i, j-1) is the adjacent cells of the north side central location (i, j).Therefore, symbol (i, j-1/2) It is the boundary between central location in y-direction and the adjacent cells of north side.

The same symbol is used to the adjacent cells in southern side：(i, j+1/2) indicates central location in y-direction and southern side phase Boundary between adjacent unit (i, j+1).For simplicity, in above-mentioned equation, when indicating the direction y (or j), increase subscript i.It will be evident that the hydraulic conductivity in defining equation (11) in a similar manner.

In fig. 14, diagonal item is limited by following equation (17)And the water guide between the unit in three directions It is limited in rate equation 11 as described above.

Item in Figure 14(right-hand vector) is indicated by equation (17b) herein.It is extracted by following formula：

Wherein, i=1,2 ... Nz, wherein the grid sum that Nz is still vertically oriented, the number of perpendicular layers.In side In the extraction of journey (17b), layer productivity index PI_iIt is limited by equation (17) and gesture item Φ_BIt is the west, east, north of central location With the known boundary gesture of the boundary of each adjacent cells in south.In general, item, unit and grid block are used interchangeably.

Conventional well model can be generally classified as three groups：(a) explicit well model；(b) well model as defined in bottom hole pressure； With fully implicit solution well model (Aziz K, Settari A, Petroleum Reservoir Simulation, Applied Science Publishers Ltd,London 1979).In order to better understand the present invention, each well model is briefly introduced.

Explicit well model

For explicit well model, the source item q in equation (10)_iIt is defined according to following equation (12)：

Wherein, q_iIt is the productivity of unit (grid block) i, wherein well passes through wherein and perforated.

Equation (12) equation (10) is substituted into for unit i to obtain

T_UpiΦ_i-1+T_C,iΦ_i+T_Down,iΦ_i+1=b_i (13)

Wherein,

T_c,i=-(T_Up,i+T_{Down, i}+T_Wi+T_Ei+T_Ni+T_Si) (14a)

And

For all unit i=1 around the well of well unit, Nz writes out equation (13) and only obtains with class shown in Figure 12 The linear system of the equation of three diagonal coefficient matrixes of type, the coefficient matrix can be written as in matrix-vector symbol：

In equation (15), A_RRIt is (Nz x Nz) tridiagonal matrix, and Φ_RAnd b_RIt is (Nzx 1) vector.Utilize three Diagonal line linear system solution device (e.g., Thomas algorithms), by well grid block therethrough carry out computer disposal come Equation (15) is solved, the unknown gesture Φ of reservoir is obtained_R。

Tridiagonal matrix and system

Tridiagonal matrix is that only in centre tool, there are three cornerwise matrixes, wherein real number or plural number are used as diagonal line In entry.These diagonal lines are referred to as " Lower diagonal line ", " center diagonal " and " Upper diagonal line ".Tridiagonal matrix In surplus element or entry be zero.For example, Figure 17 shows the tridiagonal matrix A of 8x8 matrix sizes (or rank n=8).Scheming In 17, elements A；I-th of elements A, that is, a₁,a₂,a₃,a_i, i=1,8 expression Lower diagonal lines；b₁,b₂,b₃,b_i, diagonal for center Line；And c₁,c₂,c₃,c_i, it is Upper diagonal line element.

The example trigonal system of equation is shown in Figure 18.It is with tridiagonal matrix as described above and in Figure 18 Shown in the solution of linear system of equation easily solved by Gaussian elimination method.In https:// Include non trivial solution by this method at en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm Example.

Therefore, by solve following equation (16a) executed to the matrix relationship of (16e) it is as shown below and in Figure 18 In x_iSolution：

b′_i=1；I=1,2 ..., n (16b)

Well model as defined in bottom hole pressure

Φ_wIt is use according to document (e.g., the text of Muskat, " Physical Principles of Oil Production ", McGraw-Hill Book Co. (1949) and " The Flow of Homogeneous Fluids Through Porous Media ", McGraw-Hill Book Co. (1937)) in explain technology routine techniques along Well bore is open to the uniform gesture of production.For the purpose modeled in the present context, decline quilt along the score pressure of well It is considered insignificant.Assuming that Φ_wIt is known (or defined), the oil yield (oil rate) from perforation passes through such as lower section Journey (17) calculates：

Wherein, PI_iIt is a layer productivity index, Φ w are defined base apertures gesture (pressure of base-level correction), and Φ i are wells The reservoir gridblock pressure of grid block (unit) i therethrough, r_o,iIt is referred to as the well block radius of the Pi Siman of grid block i, It is defined as 0.2 Δ x, r_wIt is the radius of well.

Variable in equation (17) is explained in " name " part above.Equation (17) is substituted into equation (10) and received There is following result for unit i in the item of collection unit i：

T_UpiΦ_i-1+T_C,iΦ_i+T_Down,iΦ_i+1W=b_i (18)

It enables

T_ci=-(T_Up,i+T_down,i+T_Wi+T_Ei+T_Ni++T_Si+PI_i) (18a)

b_i=-(PI_iΦ_W+T_WiΦ_BW+T_EiΦ_BE+T_NiΦ_BN+T_SiΦ_BS) (18b)

When writing out equation (18) for all grid block i=1, Nz, matrix system shown in Figure 14 is obtained.It can see Going out, the matrix for well model as defined in bottom hole pressure in Figure 14 is similar to the matrix of Figure 12, and in a comparable manner, Equation (18) is similar to equation (13).It can handle by using the matrix computer of three diagonal equation solver methods Easily solve well model as defined in bottom hole pressure.

Fully implicit solution well model

According to equation (19) regulation according to the total output q of the well of fully implicit solution well model_T。

Individual completion rate q is calculated by equation (17)_i.For implicit well model, well bore gesture Φ_WIt is assumed to be entire It is constant in well but it is unknown.

The item that equation (19) is substituted into equation (10) and collector unit i obtains following formula for unit (i)：

T_UpiΦ_i-1+T_C,iΦ_i+T_Down,iΦ_i+1W-PI_iΦ_W=b_i (20)

It enables

T_ci=-(T_Up,i+T_down,i+T_Wi+T_Ei+T_Ni+PI_i) (20A)

b_i=-(T_WiΦ_BW+T_EiΦ_BE+T_NiΦ_BN+T_SiΦ_BS) (20B)

Equation (20), which is write out, for all units generates the linear system with the equation of form shown in Figure 13, top pair Linea angulata solid line indicates the T such as limited by equation (11)_Up,iAnd Lower diagonal line solid line describes described in equation (11) Be referred to as T_Down,iElement.Central item T_C,IIt is limited by equation (20A) and right-hand vector b_iIt is limited by equation (20B).

The linear system of the matrix (equation 20) of Figure 13 can utilize vector matrix symbolic indication as follows：

In equation (21), A_RRIt is (Nz x Nz) tridiagonal matrix, A_RWIt is (Nz x 1) vectorial (reservoir PI), A_WR It is (1x Nz) vectorial (PI) and A_WWIt is (1x1) scalar.For the example：

It is write out with quantic：

Φ is solved according to equation 22_WIt obtains：

It substitutes into equation (21)

The item collected in equation (25) generates：

Coefficient matrix (the A of equation (26)_RR-A_ww ^-1A_wRA_RW) it is (Nz x Nz) complete matrix.

Obtained coefficient matrix can be defined as：

And

Equation (26) can be written as：

Direct solution by using Gaussian elimination method or any other conventional solver appropriate for complete matrix The matrix of equation (28) can be solved.

If in implicit well model, the quantity Nz of layer is big, and well complete completion in all layers, then equation (28) Solution becomes expensive on calculating the time.Thus and in addition, if being related to many wells, equation is usually solved by alternative manner (28).The example is in " A Fully-Implicit Fully-Coupled Well Model for Parallel Mega- Cell Reservoir Simulation”,SPE Technical Symposium of Saudi Arabia Section, It is described in 14-16 May 2005.

Flow chart I (Figure 10) indicates the basic of the fully implicit solution unity couping well model simultaneous solution of matrix-type shown in Figure 13 Computer disposal sequence.During step 100, start to simulate by reading reservoir data and creation data.Reservoir data includes Reservoir geological information --- the degree (length) in its size, x, y and z directionss, and reservoir attribute (e.g., Permeability Distribution, porous point Cloth), the thickness of layer, relative permeability data, capillary pressure force data, fluid properties data (e.g., fluid density table, stratum body Product factor table, viscosity table), the position of the well in reservoir, the position of well perforation.

Limited in step before creation data is included in is well measurements or defined oil, the production of water and natural gas Rate.Creation data further includes the minimum bottom hole pressure of each well.

In many cases, if natural gas creation data is unavailable, only input oil productivity and aquatic yield.Such as Fruit does not generate water in oil field, then only Oil Generation yield is read as input.

During the step 102, time step increases by one, and the nonlinear iteration executed during this time step number The iteration count of amount is arranged to zero.During step 104, the Jacobian matrixes of reservoir data are formed.In step 106, Then sparse preprocessor (Yousef Saad, Iterative Methods for Sparse are used by alternative manner Linear Systems,Society of Industrial&Applied Mathematics(SIAM)Publication, 2003) the result linear system of equation (19) is solved.

During step 108, convergence step is executed to determine whether nonlinear iteration has restrained.What step 106 obtained The individual residual error (residual) of equation is examined for user-defined tolerance (tolerance).If these tolerances are full Foot then exits nonlinear iteration cycle, and for current time step, solution output is written into file, processing back to step 102 with For time step to be increased, continue on for the processing of increased time step, as shown in the figure.If during step 108 User-defined tolerance is confirmed as being unsatisfactory for, then back to step 104 and is continued according to the processing of nonlinear iteration cycle.Such as The number of fruit nonlinear iteration becomes too big, then can determine the size for reducing time step, and return to step 102.

However, being based on pretreated intensity, this method can also be very expensive on calculating the time, because there is no indicate The exact way of complete matrix in equation (19).For difficult equation, in the case of extreme inhomogeneities and small layer thickness, Alternative manner can not restrain.

In addition, for the highly non-uniform reservoir with some vertical non-layers communicated, above-mentioned well model does not generate just True physical solution.On the contrary, they can generate incorrect fluid flow profile, and in some cases, they can result in simulation The problem of device is restrained.

The present invention

In the case of some vertical fluid obstacles, explicit and fully implicit solution model can generate entirely different fluid and cut open Face.This is shown in Figure 5.As shown in figure 5, reservoir model V is made of the upper layer 50 of opposite low-permeability, and in vertical current In the case of, it is located at 52 top of high permeability layer of isolation, is communicated without the vertical current with adjacent layer.Fluid barrier layer 52 is located at Medium permeability in reservoir V and the top of layer 54 that is communicated with vertical current.As that can see in Figure 5, well model V Productivity qT be identical for implicitly and explicitly well modeling method.

However, productivity q1, q2 and q3 of layer 50,52 and 54 are significantly different.The productivity q2 of layer 52 is contributed largely in logical Cross the total output qT shown in the curve 56 of explicit model.On the contrary, for the implicit model shown in by curve 58, the life of layer 52 Yield q2 is obviously less.

In fully implicit solution well model, the productivity of layer 2 is obviously less, this is because following facts：This method considers internal Inhomogeneities around reservoir inhomogeneities (all reservoir attributes in reservoir) and well block other than index of perforating is (or aobvious The property of layer 2 is used only in formula method, as the unique data that rate fraction is distributed to the perforation).For example, fully implicit solution well method Focus on being supplied to layer 2 from layer 1 and layer 3 without fluid due to the fluid-tight obstacle between layer 2 and layer 1 and 3.Therefore, one Denier generates some fluids from well model layer 2, and the pressure of layer 2 should decline and the layer should not be supplied to high-speed in well, Although the layer has very high permeability.

On the other hand, explicit method is based only upon the permeability of this layer without considering to connect with the layer of the upper and lower to distribute The rate of layer 2.Based on this, explicit method will be the very high rate of the Layer assignment and will be kept in next simulated time step-length It.In simulated time step-length later, this will lead to the unstability of productivity, that is, simulator is grown up time step is reduced It is small and will take a very long time completion simulation.For identical input data, reservoir simulation makes model be based on selection The modeling technique that uses and provide different the result is that unsatisfactory.

According to the present invention, more comprehensive well model is provided.Well model according to the present invention is referred to as coupling reservoir well model. Associated numerical solution is referred to as fully implicit solution unity couping and simultaneous solution.Fully implicit solution unity couping reservoir well model is generated along perforation Well interval correct fluid flow profile, as will be described.As shown in Figure 4, reservoir model O is a individually by the i of quantity z Layer 1 to Nz (limit as shown in Figure 4 middle finger, each layer have permeability k_x,iWith thickness deltat z_iAnd gesture Φ_i) and it is opposite The upper and lower 40 of low-permeability are constituted, and in the case of vertical current, are located at 42 top of high permeability layer of isolation It is communicated without the vertical current with adjacent layer with lower section.As that can see in Fig. 4, the productivity q of the layer i of well model O_iIt is logical Cross the expression formula instruction listed in Fig. 4.

Fully implicit solution coupling reservoir well model V propositions in above-mentioned equation (20) of Fig. 5, and also herein further to retouch The matrix format stated proposes：

The present invention is based on following facts：The bottom hole pressure and the system according to the present invention of layered reservoir formation with vertical well It is identical.By identifying fluid barrier and in order to which the reservoir communicated in vertical direction around well is combined by processing intent Or merge, form the system according to the present invention.It care should be taken to when forming simplified system according to the present invention.It is necessary to simplify system It is properly formed, otherwise, the sum of nonlinear Newtonian iteration can be increased by simplifying the mistake during system is constituted.

The system according to the present invention solves bottom hole pressure.By being that the execution of defined bottom hole pressure well should by well processing Solution.Process is fully implicit solution；However, it is not simultaneous solution.On the contrary, solution is continuous.The method of the present invention is convergent, because It is a part for the global Newton iteration of simulator.Therefore, if model according to the present invention is correctly configured, in rate calculations In any possible error will be small and will reduce with simulator Newton iteration.

Flow chart F (Figure 11) shows basic computer processing sequence according to the present invention and in the reservoir using the present invention During the exemplary embodiments of the continuous fully implicit solution well model (sequential fully implicit well model) of simulation The computational methods of generation.

During step 200, started by reading reservoir data and creation data, simulation.It is read during step 200 Reservoir data and creation data are above-mentioned types.During the step 200, also initialize reservoir simulator, will the simulation date and Time step is set as zero.During step 202, time step increases by one, and is executed during this time step non-thread The iteration count of the quantity of property iteration is arranged to zero.

During step 204, the Jacobian matrixes of reservoir data are formed.In step 206, it is formed according to the present invention such as The simplification system limited by above-mentioned model R, and in the manner described above for base apertures gesture Φ_wSolution matrix.In step 208 phase Between, solve the linear system matrix A of modification in the above described manner according to the present invention.

Figure 15 shows to form simplified well model system matrix R and solves base apertures gesture according to the step 204 of Figure 11 and 206 Φ_wMethod.As step 210 indicates, the vertical current obstacle in original well model system is identified.This can be based on log data Or by being completed according to the regulation of the data in initial reservoir model by Reservoir Analysis person.

After step 210, then during step 212, by data processing system D-shaped at simplified well model system System.In well model be located between stream barrier layer and be grouped together for analysis model with those of vertical current layer Purpose.

Then, during step 214, pass through computer disposal, the technology discussed using equation (17), (17a) and (17b) Obtained simplification well model system is solved, base apertures gesture Φ is obtained_wWith reservoir unknown number.Then step 216 utilizes data processing System D solves the full well model system structure matrix of equation (27) using direct solution or other above-mentioned proper technologies.So Completion rate q is determined based on the result of step 216 using data processing system D afterwards_iWith total well stream speed q_T。

Referring again to Figure 11, during step 220, convergence step is executed to determine whether nonlinear iteration restrains.Step The individual residual error of 216 obtained equations is examined for user-defined tolerance.If these tolerances are satisfied, exit non- Linear iteraction recycles, and solution output is written into file in current time step, and handles increased for waiting for back to step 202 Time step, and continue for the processing of increased time step, as shown in the figure.If user-defined during step 208 Tolerance is confirmed as being unsatisfactory for, then back to step 204 and is continued according to the processing of nonlinear iteration cycle.If non-linear change Generation number becomes too big, then can determine adjustment model.

Horizontal well

Figure 19 shows the horizontal well 300 being orientated in the y-direction in three dimensional reservoir framework H.The model H of Figure 19 is similar to Fig. 6, difference are the Finite Difference Meshes of horizontal well 300.Well 300 is located at extending in the y-direction along vertical plane or prolongs in the z-direction That stretches descends in a series in the center of the central location in each of stratomere 304.Figure 24 A and Figure 24 B show the level of Figure 19 The symbol of grid cell in well model H.As shown in Figure 19, well 300 passes through N in the horizontal or directions y_yA unit completion, and Gesture in adjacent cells：ΦB_Up, Φ B_E, Φ B_Down, Φ B_WIt is constant and according to dry run (time step before or iteration Value) it is known.Here, subscript B refers on " boundary ", and the adjacent cells of Up instructions top, E indicates the adjacent cells in east side, down Indicate the adjacent cells of lower section and the adjacent cells in W instructions west side, wherein furthermore, Φ describes fluid potential or base-level correction Pressure.In Figure 19, well 300 extends through a series of each of grid blocks 302 in the directions y along longitudinal axis horizontal.If Know total well rate q_TIt is input in reservoir simulator, then reservoir simulator is calculated in Figure 19 for each time step and shown Each grid block 302 gesture value Φ.

As that will state, known total rate q is added up to being formed_TLayer perforation rate q_iMeasurement in terms of, this hair The bright computer performance for improving reservoir simulator.

Similar to equation (16), the perforation rate q of each grid block_iIt can be represented as：

q_i=PI_i(Φ_i-Φ_w) (29)

In equation (29), productivity index is dominated by following formula：

The above-mentioned expression formula and its physical relation of reservoir parameter are according to Chen and Zhang, " Well Flow Models for Various Numerical Methods”.International Journal of Numerical Analysis and Modeling,Vol 6,No 3,pages 378-388。

Figure 21 and Figure 22 indicates flow chart G jointly, show basic computer processing sequence according to the present invention and The typical case such as the continuous fully implicit solution well model for the Model of Horizontal Well being shown in FIG. 19 using the reservoir simulation of the present invention is real Apply the computational methods occurred during example.During step 400, by the reservoir model H in initialization data processing system D and read Reservoir data and creation data, simulation is taken to start.The reservoir data and creation data read during step 400 is to be related to flow Scheme the above-mentioned type of F and the vertical well model of step 200 (Figure 11).

During step 400, the unit gesture Φ of the k unit for Model of Horizontal Well H is formed_kEstimation.In Figure 19 It schematically indicates, k=nx*ny*nz.During step 402, estimating for wellbore gesture is formed according to the equation (32) being listed below Meter：

Wherein, it is distributed come the Φ in accounting equation (32) according to the initial gesture determined in step 400_i。

During step 404, the perforation rate q of each unit I is determined according to the relationship indicated in equation above (29)_i。 In a step 406, reservoir simulator is initialized, and simulator iteration count υ is arranged to zero.For primary iteration, when simulation Between step-length t be also initialized to zero.During subsequent iteration before executing step 406, iteration count and time step meter Number then increases, as that will be set forth below.Three-dimensional is solved according to the equation (1) of single-phase petroleum streams by using reservoir simulator Potential equation determines the grid block gesture of initial time step-length during step 406.During step 408, for each perforation (i) (i=1,2, Ny) determines each perforation Φ_BW, Φ_{B, Up}, Φ_BE, Φ_{B, Down}The boundary gesture of surrounding, and boundary gesture is stored in data In the memory of processing system D.

Then step 410 is executed with by by the grid block 302 of flow obstacle is not combined (as schemed in-between In 20A 306 at show), form simplified well model H-1 (Figure 20 B).In order to execute step 410, identification such as original level In well model system H 303 at the bottom horizontal flow sheet obstacle that shows.It can be based on log data or by Reservoir Analysis person according to original The regulations of data in beginning reservoir model is completed.The well model system model of the simplification formed during step 410 can lead to well Layer 302 those of in model H, these layers 302 are between flow obstacle grid block or layer 303 and have water between them Advection is dynamic, these layers 302 are grouped together for analysis model purpose.

As in the example of Figure 20 A and Figure 20 B it can be seen that, as step 410 as a result, Figure 20 A flowing hinder The Model of Horizontal Well H with two groups of four grid blocks 302 on the opposite side of grid block 303 is hindered to be transformed to simplified well model H-1, Simplified well model H-1 has on every side of flow obstacle grid block 303 there are two grid block 308.

Therefore the simplification well model equation of the simplification well model system H-1 in Figure 20 B becomes as shown in Figure 20 C：

Wherein,

Then, during step 410 (Figure 22), simplified well model is solved, is obtained And Φ_w, wherein Indicate the grid block gesture of simplified Model of Horizontal Well H-1.

During step 412, Φ is used_wModel of Horizontal Well is converted into fixed flowing base apertures gesture well.Step 414 relates to And the tridiagonal matrix system of Figure 20 C of the Model of Horizontal Well H-1 simplified shown in Figure 20 B is solved, to determine in Figure 19 The gesture of each grid 302 shown.

Referring again to Figure 22, during step 416, convergence step is executed to determine whether nonlinear iteration has restrained. The individual residual error for the equation that step 414 obtains is examined for user-defined tolerance.If these tolerances meet, exit Nonlinear iteration recycles, and solution output is written into file in current time step, and is stored in the memory of data processing system D In, and handle back to step 418 for waiting for increased time step, and the processing for increased time step proceeds to Step 406, as shown in the figure.If user-defined tolerance is confirmed as being unsatisfactory for during step 416, changed according to non-linear The processing of generation cycle back to step 406 and continues.If nonlinear iteration number becomes too big, can determine to adjust model.

Multiple vertical wells

3 dimension reservoir model M are shown in FIG. 23 as with multiple vertical wells 500, it is assumed that each slightly may be used with single-phase The petroleum streams of compression.According to the name tissue model in the symbol being illustrated schematically in Figure 25 A, Figure 25 B, Figure 25 C and Figure 25 D M。

The model M of Figure 23 is similar to Fig. 6, and difference is the Finite Difference Meshes of several vertical wells 500.Well 500 is respectively positioned on edge The central location descended in a series in each of stratomere 504 that is horizontal-extending in the z-direction or extending in the y-direction with the directions z Center in.

As shown in Figure 6, each well 500 in the vertical or directions z by Nz unit completion, and the gesture in adjacent cells： ΦB_N, Φ B_E, Φ B_S, Φ B_WIt is constant and is known according to dry run (time step before or iterative value).Such as Fig. 6 Shown in, here, subscript B refers on " boundary ", and the adjacent cells of N instructions top, E indicates the adjacent cells in east side, under down instructions The adjacent cells in adjacent cells and W the instruction west side of side, wherein furthermore, Φ describes the pressure of fluid potential or base-level correction Power.

Each well 500 extends through a series of each of grid blocks 502 in the directions z along longitudinal axis horizontal, and in various depths At particular formation at degree, each well 500 is divided by completion 506.Figure 25 A to Figure 25 D schematically show the well 50 in Figure 23 The number symbol of well 1 and well 2 is designated as in model M.

If known total well rate qT of each well 500 is provided as the input parameter into reservoir simulator, Reservoir simulator calculates the gesture value Φ of each grid block 502 shown in Figure 23 for each time step.

For the ordinary circumstance of the vertical reservoir model M of 3 dimensions, n_wIndicate the quantity of well, and total well rate of each well 500 By q_T(l) it provides and indicates, l=1, n_wIt is known according to creation data and is provided as the input into reservoir simulator Parameter.

Figure 26 indicates the flow chart of processing according to the present invention, wherein multiple vertical wells is presented, such as in the model M of Figure 23 In situation.Therefore, for the vertical reservoir model M of 3 dimensions, in step 602, reservoir simulator is initialized, and read from memory Take reservoir and production for processing.This by with the step 200 of Figure 11 it is comparable in a manner of complete.In 602 be shown in FIG. 26 Simulation process continues, wherein the estimation Φ w (l), l=1,2 of the wellbore gesture of each well are formed according to equation (32) ... nw, In, productivity index PI is determined according to the measurement expressed in equation (17).During step 604, based on what is obtained from step 602 The estimation of wellbore gesture calculates the perforation rate q of each well_i(l), l=1 ... nw.

Step 606 is related to being formed each well l=1 of the multiple well model Ms of 3 dimensions, and 2 ... the well model of the simplification of nw.Pass through Combination merges the adjacent well unit 502 in stratum to execute reservoir simulation, and the adjacent well unit 502 has between each other It is in fluid communication and is also located between the flow obstacle layer that does not flow therebetween.This is completed in a manner of three-dimension layer, such as, such as to scheme Schematically shown in 1A and Figure 1B with the mode of the single neighbouring layer of vertical low well complete and Figure 20 A and Figure 20 in it is schematic The mode of the layer neighbouring with bottom horizontal flow sheet well shown is completed.

In step 608, the simplification well model based on step 606, with the letter of Figure 20 C of the Model of Horizontal Well H with Figure 19 Change the comparable mode of well matrix, forms the simplification well matrix of each well.

Then, in step 610, the simplification well matrix obtained by solution procedure 608, determines the base apertures gesture of each well Φ w (l), l=1,2 ... nw.In step 612, according to equation 18 above is related to, the relationship of 18A and 18B are each to change The diagonal line of the principal matrix of wellWith right-hand vector b_i?.Figure 27 is to use to solve in Figure 24 A and Figure 24 B and Figure 25 A to Figure 25 D The simplification well of the linear system of the equation of two wells of the simplification for the numbering system released, three dimensional reservoir framework or 3x3x2 reservoir models The example schematic diagram of matrix.

During step 614, for all grid blocks 502 of the unknown number of model M, exemplary total square such as Figure 27 is formed Battle array.Then the processing after the step 614 of Figure 26 proceeds to the convergence test of the mode of the step 416 of Figure 27, and iteration and Time step increases with back to step 602 for being further processed including n_wTotal system in each of a well changes Generation.

Data processing system

As shown in Figure 16, data processing system D according to the present invention includes computer 240, with processor 242 Be coupled to processor 242 memory 244 with by operational order, control information and data library record storage is wherein.It is necessary When, computer 240 can be portable digital processor, e.g., the personal computer of following form：Laptop computer, notebook meter Calculation machine or other programmings appropriate or programmable digital data processing equipment, e.g., desktop computer.It is also understood that meter Calculation machine 240 can be that there is the multi-core processor of node (e.g., to come from Intel Company or United States advanced Micro Devices (AMD)) or the mainframe computer of any general type of proper treatment ability (e.g., is purchased from Armonk, the U.S. world quotient of N.Y. With machines corporation (IBM)) or other sources.

Computer 240 has user interface 246 and Output Display Unit 248, is used for showing the survey executed according to the present invention The output data or record of the processing of well DATA REASONING, to obtain measured value and be formed in (one or more) well of subsurface formations Stratum determination well production model.Output Display Unit 48 includes component can carry such as printer and output display screen Visual output for printout information or following form is recorded as output or image：Figure, tables of data, graph image, number According to figure etc..

The user interface 246 of computer 240 further includes user input apparatus appropriate or input/output control unit 250 It is accessed with providing user to control or access information and data-base recording and operate computer 240.Data processing system D further includes The database 252 of storage in computer storage, computer storage can be internal storage 244 or external networking or non- The memory of networking, as indicated at 254 in associated database server 256.

Data processing system D includes the program code 260 being stored in the non-transitory memory 244 of computer 240.Root Program code 260 according to the present invention is the form of the operable instruction of computer, makes data processor 242 to have already mentioned above Mode form the continuous fully implicit solution well model of reservoir simulation according to the present invention.

It should be noted that program code 260 can be microcode, and program, the shape of routine or the operable language of symbolic computer Formula, the function of providing specific one group of control data processing system D and guide the operation of its sequence operated.Program code 260 instruction can be stored in the memory 244 of computer 240 or in computer disk, tape, conventional hard disk driving, electricity Sub- read-only memory, light storage device or be stored on it computer can use non-transitory medium other numbers appropriate According on storage device.Program code 260, which is further included on data storage device (e.g., server 64), is used as non-transitory meter Calculation machine readable medium, as shown in the figure.

Two illustrated examples model equations are presented below：Seven layers of uniform formations, there are one be broken flow obstacle for tool (Fig. 7)；With two Floor 12 uniform formations, there are two fracture flow obstacles (Fig. 9) for tool.

Seven layers of uniform well model

Fig. 7 A show seven reservoirs and attribute of archetype 70 and simplified model 71.As shown in the figure, it is assumed that reservoir has Seven layers.The layer thickness of each layer 72 with vertical current is 10ft (foot).In the presence of indicating the not fracture of vertical current and thickness Degree is the layer 73 of 1ft.It is further assumed that layer 73 is not communicated with above and below layer 72.Assuming that in centre, there are vertical wells, such as It is indicated by arrow.For the model of Fig. 7 A, initial reservoir gesture (pressure of base-level correction) is 3,000psi.Assuming that each Plane permeability k of the layer 72 with 10md_xAnd k_yAnd the vertical permeability k of 1md_z。

Table 1 summarizes the reservoir and grid property of model 70.Assuming that the sizing grid on in-plane (square net) is 840ft.Oil viscosity is arranged to 1cp, and assumes that petroleum-bearing formation volume factor is 1.The total Oil Generation yield of well is arranged to 1, 000B/D.The layer productivity index PI of each layer completion is calculated by Pi Siman methods, as described, and also in table 1 It shows.

Table 1：Equation 1 --- reservoir attribute

Fully implicit solution unity couping simultaneous solution

Reservoir pressure and bottom are formed with similar fashion illustrating above for equation (18-19) and in figure 13 illustrates The coefficient matrix of the solution of portion's pore pressure force.It can recognize there are only 8 unknown numbers (7 gesture or the pressure and one of base-level correction A base apertures gesture), and coefficient matrix is non-sparse.For unknown number reservoir (layer) gesture Φ_i, i=1,7 and other unknown numbers Φ_W, the linear system of equation can be solved by direct method (e.g., Gaussian elimination method).

As a result

The layer gesture of the model 70 of 2 overview diagram 7A of table being computed, well bore gesture and layer (completion) flow velocity.

The accurate solution of table 2 --- original equation

Seen according to result of calculation, the base apertures gesture Φ of calculating_W=1257.36psi.

The formation of equation according to the present invention

According in table 1 and reservoir data as shown in Figure 7A, there is only communicate a layers not vertical with other layer 73.Therefore, as shown in Figure 7 B, according to the present invention, the layer 72 of 73 top of zone of fracture is combined into single layer to form simplified well Model.Similar, the layer 72 of 73 lower section of layer is combined into single layer.Now it can be seen that the model simplified has only three layers.Phase Than 8 in full model, the sum of unknown number is 4.

Table 3 summarizes the attribute for the simplification well model 71 that processing according to the present invention is formed.

Table 3 --- simplify well model

Layer	Thickness, ft	Kx=Ky,md	Kz,md	PI,b/d/psi
					1	20	10	1	0.24
2	1	10,000	0	12.14
					3	40	10	1	0.40

The linear system for simplifying the equation (equation 20) of system also has non-structural coefficient matrix, but has 50% or less The unknown number of quantity.In actual reservoir, in the case of hundreds of layers and only some flow obstacles, well mould according to the present invention The reduction of type size will be strong, for example, simplified well model system model according to the present invention can be system-wide size 1 percent.Simplified system is solved for layer gesture and base apertures gesture by direct solution.Table 4 lists result.

Table 4 --- simplify the result of system

It can be seen that the base apertures gesture being computed：Φ_WThe Φ that=1257.36psi is calculated with full model_WIt is identical.

Determining well gesture is the unique information needed for next step-length.The well is then processed into defined bottom hole pressure (gesture) model.Computer disposal followed by the process according to matrix and equation (16 and 18) description using Figure 14 is flowed with calculating Measure section (layer rate) and total well rate.In fig. 14, the Upper diagonal solid line of matrix is indicated through equation (11) restriction T_Up,i, and the Lower diagonal solid line description of matrix is also referred to as T by what equation (2) limited_Down,iElement.Central item T_C,iIt is logical Equation (17a) restriction is crossed, and right-hand vector bi is limited by equation (17b).

Summary results in table 5.It is noted that the well rate of total calculating and the input value of 1,000b/d are identical.

Table 5 --- utilize the result of the total system of the present invention

Layer	Gesture, psi	Rate/d
			1	2630.61	166.67
2	2630.61	166.67
			3	1257.37	0.0
4	2630.61	166.67
			5	2630.61	166.67
6	2630.61	166.67
			7	2630.61	166.67
It amounts to		1,000

The result presented in table 5 is identical as the table 1 of fully implicit solution well model.Between the well rate of well calculate and input Poor or error is zero in this situation, and does not need additional iteration.This is because following facts：Reservoir is uniform, and in shape At the error for not generating roughening when simplifying system.Matrix diagonals line element and right-hand vector are identical with Figure 14, that is, Lower diagonal is real Line indicates the T limited by equation (11)_up,i, and the Upper diagonal solid line of matrix describes and above-mentioned is referred to as T_Down,iElement. Central item T_C,iIt is limited by equation (17a), and right-hand vector bi is limited by equation (17b).

The item PI occurred on equation (17) is the perforation productivity index of square net, is defined by following formula：

Wherein, r_wIt is wellbore radius.

With the comparison of explicit well model

In several reservoir simulators, half implicit well model or explicit well model are used.If the formulation of well model is Semi-implicit but it is collapsed under pressure variations to be explicit, then the formulation collapses to explicit well model.Explicit well model is directed to The equation obtained by the computer processing procedure of the matrix and equation (12-14) of deferring to Figure 12.In fig. 12, in diagonal line The item T occurred on element_Down,i,T_Up,iIt is limited by equation (11), and T_c,i,b_iIt is limited by equation (14a) and equation (14b). It is mutually isostructural with the model of Fig. 7 A and Fig. 7 B that Fig. 8 A show that seven reservoirs of implicit well model 80 and attribute and Fig. 8 B are shown Explicit well model 81.As shown in the figure, it is assumed that reservoir has seven layers.Layer 82 all has the gesture Φ of 2630psi.Layer 83 has The gesture Φ of 1257psi, layer 83 indicates fracture and the layer 72 being further assumed not with above and below communicates.Assuming that being deposited in centre In vertical well, as indicated by arrow.Fig. 8 A and Fig. 8 B compare the result of implicitly and explicitly model.Summarize in table 6 and calculates Perforation (layer) rate.

The comparison of perforation (layer) rate of table 6 --- different well models

As shown, the model 81 according to explicit method is inaccurate；Its complete erroneous calculations perforation rate.Explicitly The practical distribution of model method is come all wells production of the thin zone of fracture 83 indicated in the fig. 8b freely, because the layer is with highest Productivity index.

The implicit method simplified model of the present invention (computation-intensive fully implicit solution model according further to) without this distribution, It is opposite to be determined that layer 83 come from above and the support of the fluid of the layer of lower section 82.Show at 83 present in actual reservoir The onty fiuidic that the such zone of fracture gone out can obtain supports the adjacent cells from its plane.However, due to fracture Layer is very thin layer, the hydraulic conductivity on these directions study carefully its person's character for be small.Therefore, zone of fracture cannot be supported by aobvious Fluid under the rate of formula modeling.

In fact, conventional implicit well model was shown during transition time, zone of fracture supports most wells production, such as explicit Method is such.However, stressor layer rapid decrease in layer 83 and assuming the value (constant bottom hole pressure) of uniform well bore gesture. After the pressure reduction, reach stable state, and well productivity actually by from the top of perpendicular flow obstacle 83 and The contribution of the layer 82 of lower section generates.

Two Floor 12 uniform formations models

Model meshes system 90 includes 22 layers as shown in Figure 9.The position of the zone of fracture 6 and 12 of high permeability It is seen as moving downwards through layer and schematically be indicated at 91 and 92.There are 5 layers 93 above layer 91, number be 1 to 5, each layer has vertical current.There are five layers 94 with vertical current also in model 90, in flow obstacle layer 91 and 92 Between, and ten with the vertical current layer 95 positioned at 92 lower section of flow obstacle layer.The reservoir number of model 90 is shown in table 7 According to.

Table 7 --- the reservoir data of 22 layers of equation

Plane adjacent cells permeability=20mD

Total well productivity=2,500B/D

The well of completion in all layers.

As a result：

Fully implicit solution unity couping simultaneous solution

The base apertures gesture Φ of calculating_W=1421.247psi

Table 8 --- gesture is distributed, psi

Layer	ΦW	Φi
			1	1421.25	2901.32
2	1421.25	2900.82
			3	1421.25	2900.38
4	1421.25	2900.38
			5	1421.25	2900.38
6	1421.25	1421.25
			7	1421.25	2864.43
8	1421.25	2864.37
			9	1421.25	2864.10
10	1421.25	2863.88
			11	1421.25	2864.03
12	1421.25	1421.25
			13	1421.25	2841.62
14	1421.25	2841.58
			15	1421.25	2841.48
16	1421.25	2841.33
			17	1421.25	2841.13
18	1421.25	2840.65
			19	1421.25	2840.08
20	1421.25	2839.65
			21	1421.25	2839.37
22	1421.25	2839.24

Then, the flow-data that model 90 is subjected to the explicit model method and technology of the above-mentioned type and determines.It is proposed in table 9 The comparison of the fully implicit solution and the velocity flow profile explicitly handled of the reservoir model 90 of the technology described before use.

The comparison of table 9 --- fully implicit solution unity couping well method and the flow velocity of explicit well method

Layer	Implicitly	Explicitly
			1	35.93	14.56
2	89.78	36.39
			3	53.85	21.83
4	179.48	72.78
			5	89.74	36.39
6	0.00	727.80
			7	105.09	43.67
8	52.54	21.83
			9	157.60	65.50
10	210.10	87.34
			11	87.55	36.39
12	0.00	727.80
			13	129.29	54.59
14	129.28	54.59
			15	129.28	54.59
16	129.26	54.59
			17	129.24	54.59
18	158.49	66.96
			19	158.42	66.96
20	158.37	66.96
			21	158.34	66.96
22	158.33	66.96

Simplified model is constituted

Since there is only two perpendicular flow barrier layers 91 and 92, the layer 93 of 91 top of layer in fig.9 can be combined into One layer；Layer 94 below layer 91 is combined into a layer, and the layer 95 below layer 92 is combined into another single layer.Therefore, The sum of layer according to the present invention is 5.The attribute of simplified model is as follows：

Table 10 --- simplify the attribute of well model

The base apertures gesture of calculating

Φ_W=1421.34psi

Simplified model according to the present invention continues to determine base apertures gesture Φ_w.There are five the results of the simplified model of layer to exist for tool It is listed in table 11.

Table 11 --- gesture is distributed

Layer	Pot wf	Pot
			1	1421.34	2900.53
2	1421.34	1421.35
			3	1421.34	2864.16
4	1421.34	1421.35
			5	1421.34	2740.61

By using the base apertures gesture Φ calculated according to simplified model_wAnd the defined base apertures gesture Φ using full model_w Gesture is calculated, completion layer rate is calculated according to equation (17).As a result it is indicated below in table 12.

The well layer rate of table 12 --- calculating

The q newly calculated_t=2499.84b/d

Error=2,500.-2499.8488

=0.15b/d

The error of total rate and the bottom hole pressure of calculating disappear with the nonlinear Newtonian iteration of simulator.It compares In the fully implicit solution unity couping treatment technology acquisition by the prior art as a result, the present invention obtains the life for having and can receive accuracy The simplified model of yield, but wherein, model complexity and computer processing time are substantially reduced.

The reservoir and the special of well non trivial solution that the present invention as already described above need not be used to couple linearly are asked Solve device.In contrast, the coefficient matrix of the reservoir of the coupling for using before and well equation does not have the sparsity structure of rule. Therefore, the reservoir of the coupling of general type and well equation need special solver that is expensive and can facing convergence equation.

It can be seen that as described above, any special solver of reservoir and well non trivial solution that the present invention need not couple. Utilize the identical solver for reservoir equation.Only to the modification of coefficient matrix in diagonal term.

The present invention, which solves wherein vertical well, has the reservoir simulation equation of many completions (this is commonplace in reservoir).? It is very universal with wells more than 100 perpendicular layers (completion) in nearest analog study.Unity couping with simultaneous solution is entirely hidden Formula well model is very expensive in these cases.The present invention can save a large amount of computer times.

The present invention is highly useful for the well with hundreds of perforation of completion in the reservoir of high uniformity.Pass through identification Be advantageously employed involved physical principle, the present invention will be with the well modeling mould in the reservoir of a large amount of layer (completion) equation Quasi- big time-consuming equation is reduced to small equation.In the present case, it has been found that, the layer vertically communicated can be by Synthesize single layer.The simplified model being thusly-formed retains bottom hole pressure identical with original full model.Simplify mould once solving Type obtains bottom hole pressure, and then the well in big system is just treated as specific bottom hole pressure, and is asked by conventional linear Solution device easily solves.Therefore, the present invention eliminates the needs to being written and obtaining non-structural linear solution device, the non-structural line Property solver be used for many wells with hundreds of completions and may be expensive.

The present invention has been fully described it so that the people of the average knowledge with the theme can reproduce and obtain the present invention In the result mentioned.Nevertheless, any technical staff in the technical field of subject of the present invention can execute and ask herein In the modification that does not describe, by these modifications applied to determining structure, or in its manufacturing process, need in appended right The theme being claimed in it is required that；These structures should be included within the scope of the invention.

It should be noted and appreciated that, in the feelings for not departing from the spirit or scope of the present invention described in appended claims Under condition, the present invention detailed above can be improved and be changed.

Claims (23)

1. a kind of the method implemented by computer forming model, which refers to utilizing coupling during the service life of subsurface reservoir For well reservoir model with time step during the reservoir simulation that well produces, what is determined according to the production of total well of measurement comes from underground The well productivity of composition fluid of multiple wells in reservoir and the layer completion rate of the perforated interval in the well of determination Model, the coupling well reservoir model are organized into reservoir grid, which is subdivided at the perforated interval Multiple reservoir units, the perforated interval in the reservoir is located at multiple stratum, the stratum have with it is described when Between step-length composition fluid unknown well gesture and fluid completion rate, and the stratum includes the vertical of fluid therefrom perpendicular flow Fluid flow layer and flow obstacle layer without fluid therefrom perpendicular flow, the stratum also has permeability, thickness and layer gesture, institute State multiple well units at position of the coupling well reservoir model also with the well in the stratum of the reservoir, the calculating The method that machine is implemented determines the layer completion rate of the composition fluid on the stratum from the well and from the well The well productivity of the composition fluid, described the method implemented by computer include the following steps：

(a) the full calculating matrix reservoir model of the reservoir data of the unit of the model, the full calculating matrix reservoir model are formed It is included in the reservoir data of the reservoir units at the perforated interval, the reservoir data includes the institute on the stratum State permeability, thickness and gesture；

(b) by by it is in the reservoir model, therebetween with vertical fluid flowing and between flow obstacle layer The data group of the vertical fluid fluidized bed of vertical fluid fluidized bed becomes the single perpendicular flow layer in the matrix, shape At simplified well model system matrix；

(c) bottom hole pressure of the well is determined；

(d) coupling reservoir well model is formed, which includes the full calculating matrix reservoir model and the simplified well model system Well processing is the specified well of base apertures, pressure with determining bottom hole pressure, the coupling reservoir well by system matrix Model is the form of following matrix：

Wherein, A_RRIt is the tridiagonal matrix of the reservoir data, A_RWIt is the life on the stratum adjacent with the perforated interval The diagonal matrix of yield index；A_WRIt is the diagonal matrix of the productivity index on the stratum from well to the reservoir；A_WWIt is institute State the matrix of the productivity index of well；It is the matrix of the unknown reservoir gesture of the unit around the well；It is The matrix of unknown well gesture in the well bore of the well；It is the reservoir number of the reservoir units around the well According to the matrix of constant；AndIt is the matrix of the well data constant of the well；

(e) solve the coupling reservoir well model, obtain the fluid in the reservoir units on the stratum flow and With the productivity and gesture of the reservoir units of the time step in each of the stratum；

(f) the coupling reservoir well model is solved, is obtained with the time step at the perforated interval of the reservoir The well unit productivity index；

(g) based on the life of the determination of the reservoir units at the perforated interval of the reservoir of the time step The productivity index of the determination of yield index and well unit, come the vertical fluid fluidized bed for determining the well and the stream The layer completion rate of the composition fluid in dynamic barrier layer；

(h) according to described in the vertical fluid fluidized bed of the well of the time step and the flow obstacle layer The layer completion rate for forming the determination of fluid, determines total well productivity of the well；And

(i) it is formed with described in the vertical fluid fluidized bed of the well of the time step and the flow obstacle layer Form the record of the layer completion rate of the determination of fluid and total well productivity of the determination of the well.

2. the method implemented by computer according to claim 1, wherein the tridiagonal matrix of the reservoir units The diagonal line of the hydraulic conductivity including the unit for indicating the reservoir.

3. the method implemented by computer according to claim 1, wherein the step of solving the coupling well reservoir model is wrapped It includes and applies complete matrix solver.

4. the method implemented by computer according to claim 3, further comprising the steps of：

For defined tolerance, reduce the residual error for the step of carrying out self-application complete matrix solver；And

If the tolerance is met at the time step, and for each of during the service life of the subsurface reservoir Tolerance described in time step is all met, then proceeds to the step of forming record；And

If in the time step strong point, the tolerance is unsatisfactory for, then is directed to another iteration of the processing of the time step strong point, Back to step (a) and step (b) is repeated to step (h).

5. the method implemented by computer according to claim 3, further comprising the steps of：

For defined tolerance, reduce the residual error for the step of carrying out self-application complete matrix solver；And

If tolerance is met as defined in described at the time step, but is not for the lifetime in the reservoir Between each time step described in defined tolerance can be met, then increase the simulator time step；And

Back to step (a) and step (b) is repeated to step (h) for the another of the processing of increased simulator time step One iteration.

6. the method implemented by computer according to claim 1, wherein the step of solving the coupling well reservoir model is wrapped It includes：Based on the permeability on the stratum, thickness and layer gesture come the step of solving fluid flowing.

7. the method implemented by computer according to claim 1, wherein the multiple well includes multiple vertical wells.

8. a kind of number form is at model according to processing system, which refers to being stored up using coupling well during the service life of subsurface reservoir For layer model with time step during the reservoir simulation that well produces, what is determined according to the production of total well of measurement comes from subsurface reservoir In multiple wells composition fluid well productivity and determination the well in perforated interval layer completion rate model, institute It states coupling well reservoir model to be organized into reservoir grid, which is subdivided into multiple storages at the perforated interval Layer unit, the perforated interval in the reservoir are located at multiple stratum, and the stratum has with the time step Form fluid unknown well gesture and fluid completion rate, and the stratum include fluid therefrom perpendicular flow vertical fluid flowing Layer and flow obstacle layer without fluid therefrom perpendicular flow, the stratum also has permeability, thickness and layer gesture, the coupling well Reservoir model also has multiple well units at the position of the well in the stratum of the reservoir, the data processing system The layer completion rate for determining the composition fluid on the stratum from the well and the composition fluid from the well Well productivity, and the data processing system includes：

Processor executes following steps：

(a) the full calculating matrix reservoir model of the reservoir data of the unit of the model, the full calculating matrix reservoir model are formed It is included in the reservoir data of the reservoir units at the perforated interval, the reservoir data includes the institute on the stratum State permeability, thickness and gesture；

(b) by by it is in the reservoir model, therebetween with vertical fluid flowing and between flow obstacle layer The data group of the vertical fluid fluidized bed of vertical fluid fluidized bed becomes the single perpendicular flow in the matrix Layer forms and simplifies well model system matrix；

(c) bottom hole pressure of the well is determined；

(d) coupling reservoir well model is formed comprising the full calculating matrix reservoir model and the simplified well model system square Well processing is the specified well of base apertures, pressure with determining bottom hole pressure, the coupling reservoir well model by battle array With a matrix type：

Wherein, A_RRIt is the tridiagonal matrix of the reservoir data, A_RWIt is the life on the stratum adjacent with the perforated interval The diagonal matrix of yield index；A_WRIt is the diagonal matrix of the productivity index on the stratum from well to the reservoir；A_WWIt is institute State the matrix of the productivity index of well；It is the matrix of the unknown reservoir gesture of the unit around the well；It is The matrix of unknown well gesture in the well bore of the well；It is the reservoir number of the reservoir units around the well According to the matrix of constant；AndIt is the matrix of the well data constant of the well；

(e) solve the coupling reservoir well model, obtain the fluid in the reservoir units on the stratum flow and With the productivity and gesture of the reservoir units of the time step in each of the stratum；

(f) the coupling reservoir well model is solved, is obtained with the time step at the perforated interval of the reservoir The well unit the productivity index；

(g) based on the reservoir units of the time step at the perforated interval of the reservoir and well unit The productivity index of the determination determines the composition of the vertical fluid fluidized bed and the flow obstacle layer of the well The layer completion rate of fluid；

(h) according to described group of the vertical fluid fluidized bed of the well of the time step and the flow obstacle layer At the layer completion rate of the determination of fluid, total well productivity of the well is determined；And

Memory executes following steps：It is formed with the vertical fluid fluidized bed of the well of the time step and institute State the note of the layer completion rate of the determination of the composition fluid of flow obstacle layer and total well productivity of the determination of the well Record.

9. data processing system according to claim 8, wherein the tridiagonal matrix of the reservoir units includes Indicate the diagonal line of the hydraulic conductivity of the unit of the reservoir.

10. data processing system according to claim 8, wherein the processor is solving the coupling well reservoir mould The step of being executed using complete matrix solver when type.

11. data processing system according to claim 10 further includes that the processor executes following steps：

For defined tolerance, reduce the residual error for the step of carrying out self-application complete matrix solver；And

If the tolerance is met at the time step, and when for each during the service life of the subsurface reservoir Between tolerance described in step-length all met, so that the memory is executed the step of forming record；And

If being unsatisfactory for the tolerance in the time step strong point, back to step (a) and repeat step (b) to (h) for With another iteration of the processing of the time step.

12. data processing system according to claim 10 further includes that the processor executes following steps：

For defined tolerance, reduce the residual error for the step of carrying out self-application complete matrix solver；And

If tolerance is met as defined in described at the time step, but during being not the service life for the reservoir Each time step described in defined tolerance met, then increase the simulator time step；And

Back to step (a) and repeat step (b) to (h) for increased simulator time step processing it is another repeatedly Generation.

13. data processing system according to claim 8, wherein the processor is solving the coupling well reservoir mould It is executed when type based on the permeability on the stratum, thickness and layer gesture come the step of solving fluid flowing.

14. data processing system according to claim 8, wherein the multiple well includes multiple vertical wells.

15. the operable instruction of computer is stored in non-transitory computer-readable medium, institute by a kind of data storage device Stating instruction makes data processor form model, which refers to utilizing coupling well reservoir model during the service life of subsurface reservoir With time step during the reservoir simulation of well production, is produced and determined more in subsurface reservoir according to total well of measurement The model of the well productivity of composition fluid of a well and the layer completion rate of the perforated interval in the well of determination, the coupling Well reservoir model is organized into reservoir grid, which is subdivided into multiple reservoir lists at the perforated interval Member, the perforated interval in the reservoir are located at multiple stratum, and the stratum has forming with the time step The unknown well gesture and fluid completion rate of fluid, and the stratum includes the vertical fluid fluidized bed of fluid therefrom perpendicular flow With the flow obstacle layer of no fluid therefrom perpendicular flow, also there is permeability, thickness and layer gesture, the coupling well to store up on the stratum Layer model also has multiple well units at the position of the well in the stratum of the reservoir, and the computer of the storage can Operational order makes the data processor by executing following steps to determine the composition on the stratum from the well The well productivity of the layer completion rate of fluid and the composition fluid from the well：

(a) the full calculating matrix reservoir model of the reservoir data of the unit of the model, the full calculating matrix reservoir model are formed It is included in the reservoir data of the reservoir units at the perforated interval, the reservoir data includes the institute on the stratum State permeability, thickness and gesture；

(b) by by it is in the reservoir model, therebetween with vertical fluid flowing and between flow obstacle layer The data group of the vertical fluid fluidized bed of vertical fluid fluidized bed become single vertical current in the matrix Dynamic layer, to form simplified well model system matrix；

(c) bottom hole pressure of the well is determined；

(d) coupling reservoir well model is formed comprising the full calculating matrix reservoir model and the simplified well model system square Well processing is the specified well of base apertures, pressure with determining bottom hole pressure, the coupling reservoir well model by battle array With a matrix type：

Wherein, A_RRIt is the tridiagonal matrix of the reservoir data, A_RWIt is the life on the stratum adjacent with the perforated interval The diagonal matrix of yield index；A_WRIt is the diagonal matrix of the productivity index on the stratum from well to the reservoir；A_WWIt is institute State the matrix of the productivity index of well；It is the matrix of the unknown reservoir gesture of the unit around the well；It is The matrix of unknown well gesture in the well bore of the well；It is the reservoir number of the reservoir units around the well According to the matrix of constant；AndIt is the matrix of the well data constant of the well；

(e) solve the coupling reservoir well model, obtain the fluid in the reservoir units on the stratum flow and With the productivity and gesture of the reservoir units of the time step in each of the stratum；

(f) the coupling reservoir well model is solved, is obtained with the time step at the perforated interval of the reservoir The well unit the productivity index；

(g) based on the reservoir units of the time step at the perforated interval of the reservoir and well unit The productivity index of the determination, come determine the well the vertical fluid fluidized bed and described group of the flow obstacle layer At the layer completion rate of fluid；

(h) according to described group of the vertical fluid fluidized bed of the well of the time step and the flow obstacle layer At the layer completion rate of the determination of fluid, to determine total well productivity of the well；And

(i) described group with the vertical fluid fluidized bed of the well of the time step and the flow obstacle layer is formed At the record of total well productivity of the determination of the layer completion rate and well of the determination of fluid.

16. data storage device according to claim 15, wherein solve the coupling for making the processor execute The instruction of the step of well reservoir model includes the instruction using complete matrix solver.

17. data storage device according to claim 15, wherein described instruction further include so that the processor is executed with The instruction of lower step：

For defined tolerance, reduce the residual error for the step of carrying out self-application complete matrix solver；And

If the tolerance is met at the time step, and for the subsurface reservoir service life during it is each when Between tolerance described in step-length all met, so that the memory is executed the step of forming record；And

If being unsatisfactory for the tolerance in the time step strong point, back to step (a) and repeat step (b) to (h) with In with another iteration of the processing of the time step.

18. data storage device according to claim 15, wherein described instruction further include so that the processor is executed with The instruction of lower step：

For defined tolerance, reduce the residual error for the step of carrying out self-application complete matrix solver；And

If tolerance is met as defined in described at the time step, but during being not the service life for the reservoir Each time step described in defined tolerance met, then increase the simulator time step；And

Back to step (a) and step (b) to (h) is repeated, is changed for the another of processing of increased simulator time step Generation.

19. data storage device according to claim 15, wherein solve the coupling for making the processor execute The instruction of the step of well reservoir model includes：For solving fluid flowing based on the permeability on the stratum, thickness and layer gesture Instruction.

20. data storage device according to claim 15, wherein the multiple well includes multiple vertical wells.

21. a kind of the method implemented by computer forming model, which refers to utilizing coupling during the service life of subsurface reservoir For well reservoir model with time step during the reservoir simulation that well produces, what is determined according to the production of total well of measurement comes from underground The well productivity of composition fluid of horizontal well in reservoir and the layer completion rate of the perforated interval in the horizontal well of determination Model, the coupling well reservoir model is organized into reservoir grid, which is subdivided into the perforated interval Multiple reservoir units at place, the perforated interval in the reservoir are located at multiple stratum, and the stratum has with described The unknown well gesture and fluid completion rate of the composition fluid of time step, and the stratum includes the water of fluid therefrom bottom horizontal flow sheet Flat fluid flow layer and flow obstacle layer without fluid therefrom bottom horizontal flow sheet, the stratum also have permeability, thickness and layer gesture, The coupling well reservoir model also has multiple well units at the position of the horizontal well in the stratum of the reservoir, institute State the method implemented by computer determine the stratum from the horizontal well the composition fluid layer completion rate and The well productivity of the composition fluid from the horizontal well, described the method implemented by computer include the following steps：

(a) the full calculating matrix reservoir model of the reservoir data of the unit of the model, the full calculating matrix reservoir model are formed It is included in the reservoir data of the reservoir units at the perforated interval, the reservoir data includes the institute on the stratum State permeability, thickness and gesture；

(b) by by it is in the reservoir model, therebetween with horizontal fluid flow and between flow obstacle layer The data group of the horizontal fluid flow layer of horizontal fluid flow layer become single level stream in the matrix Dynamic layer, to form simplified well model system matrix；

(c) bottom hole pressure of the horizontal well is determined；

(d) coupling reservoir well model is formed comprising the full calculating matrix reservoir model and the simplified well model system square Horizontal well processing is the specified well of base apertures, pressure with determining bottom hole pressure, the coupling reservoir well by battle array Model is with a matrix type：

Wherein, A_RRIt is the tridiagonal matrix of the reservoir data, A_RWIt is the life on the stratum adjacent with the perforated interval The vector of yield index；A_WRIt is the vector of the productivity index on the stratum from well to the reservoir；A_WWIt is the horizontal well Productivity index linear scalar；It is the vector of the unknown reservoir gesture of the unit around the horizontal well；It is the vector of the unknown well gesture in the well bore；It is the reservoir of the reservoir units around the horizontal well The vector of data constant；AndIt is the vector of the well data constant of the horizontal well；

(e) solve the coupling reservoir well model, obtain the fluid in the reservoir units on the stratum flow and With the productivity and gesture of the reservoir units of the time step in each of the stratum；

(f) the coupling reservoir well model is solved, is obtained with the time step at the perforated interval of the reservoir The well unit the productivity index；

(g) based on the reservoir units of the time step at the perforated interval of the reservoir and well unit The productivity index of the determination, come determine the horizontal well the horizontal fluid flow layer and the flow obstacle layer institute State the layer completion rate of composition fluid；

(h) according to the institute of the horizontal fluid flow layer and the flow obstacle layer of the horizontal well of the time step The layer completion rate for stating the determination of composition fluid, to determine total well productivity of the horizontal well；And

(i) it is formed with the institute of the horizontal fluid flow layer and the flow obstacle layer of the horizontal well of the time step State the record of the layer completion rate of the determination of composition fluid and total well productivity of the determination of the horizontal well.

22. a kind of data processing system forming model, which refers to being stored up using coupling well during the service life of subsurface reservoir For layer model with time step during the reservoir simulation that well produces, what is determined according to the production of total well of measurement comes from subsurface reservoir In horizontal well composition fluid well productivity and determination the horizontal well in perforated interval layer completion rate mould Type, in the coupling well reservoir model tissue to reservoir grid, which is subdivided into more at the perforated interval A reservoir units, the perforated interval in the reservoir are located at multiple stratum, and the stratum has with the time step The unknown well gesture and fluid completion rate of long composition fluid, and the stratum includes the horizontal liquid of fluid therefrom bottom horizontal flow sheet Fluidized bed and flow obstacle layer without fluid therefrom bottom horizontal flow sheet, the stratum also has permeability, thickness and layer gesture, the coupling Close multiple well units at position of the well reservoir model also with the horizontal well in the stratum of the reservoir, the data Processing system determines the layer completion rate of the composition fluid on the stratum from the horizontal well and comes from the water The well productivity of the composition fluid of horizontal well, and the data processing system includes：

Processor executes following steps：

(a) the full calculating matrix reservoir model of the reservoir data of the unit of the model, the full calculating matrix reservoir model are formed It is included in the reservoir data of the reservoir units at the perforated interval, the reservoir data includes the institute on the stratum State permeability, thickness and gesture；

(b) by by it is in the reservoir model, therebetween with horizontal fluid flow and between flow obstacle layer The data group of the horizontal fluid flow layer of horizontal fluid flow layer become single level stream in the matrix Dynamic layer, to form simplified well model system matrix；

(c) bottom hole pressure of the horizontal well is determined；

(d) coupling reservoir well model is formed comprising the full calculating matrix reservoir model and the simplified well model system square Horizontal well processing is the specified well of base apertures, pressure with determining bottom hole pressure, the coupling reservoir well by battle array Model is with a matrix type：

Wherein, A_RRIt is the tridiagonal matrix of the reservoir data, A_RWIt is the life on the stratum adjacent with the perforated interval The vector of yield index；A_WRIt is the vector of the productivity index on the stratum from well to the reservoir；A_WWIt is the horizontal well Productivity index linear scalar；It is the vector of the unknown reservoir gesture of the unit around the horizontal well；It is the vector of the unknown well gesture in the well bore；It is the reservoir of the reservoir units around the horizontal well The vector of data constant；AndIt is the vector of the well data constant of the horizontal well；

(e) solve the coupling reservoir well model, obtain the fluid in the reservoir units on the stratum flow and With the productivity and gesture of the reservoir units of the time step in each of the stratum；

(f) the coupling reservoir well model is solved, is obtained with the time step at the perforated interval of the reservoir The well unit the productivity index；

(g) based on the reservoir units of the time step at the perforated interval of the reservoir and well unit The productivity index of the determination, come determine the horizontal well the horizontal fluid flow layer and the flow obstacle layer institute State the layer completion rate of composition fluid；

(h) according to the institute of the horizontal fluid flow layer and the flow obstacle layer of the horizontal well of the time step The layer completion rate for stating the determination of composition fluid, to determine total well productivity of the horizontal well；And

Memory is formed with the horizontal fluid flow layer of the horizontal well of the time step and the flow obstacle The record of the layer completion rate of the determination of the composition fluid of layer and total well productivity of the determination of the horizontal well.

23. the operable instruction of computer is stored in non-transitory computer-readable medium, institute by a kind of data storage device Stating instruction makes data processor form model, which refers to utilizing coupling well reservoir model during the service life of subsurface reservoir Well production reservoir simulation during, with time step, according to total well of measurement produce and determine in subsurface reservoir The model of the well productivity of composition fluid of horizontal well and the layer completion rate of the perforated interval in the horizontal well of determination, institute It states coupling well reservoir model to be organized into reservoir grid, which is subdivided into multiple storages at the perforated interval Layer unit, the perforated interval in the reservoir are located at multiple stratum, and the stratum has with the time step The unknown well gesture and fluid completion rate of fluid is formed, and the stratum includes the horizontal fluid flow of fluid therefrom bottom horizontal flow sheet Layer and flow obstacle layer without fluid therefrom bottom horizontal flow sheet, the stratum also has permeability, thickness and layer gesture, the coupling well Reservoir model also has multiple well units at the position of the horizontal well in the stratum of the reservoir, the meter of the storage The operable instruction of calculation machine makes the data processor by executing following steps to determine the stratum from the horizontal well The composition fluid layer completion rate and the composition fluid from the horizontal well well productivity：

(a) the full calculating matrix reservoir model of the reservoir data of the unit of the model, the full calculating matrix reservoir model are formed It is included in the reservoir data of the reservoir units at the perforated interval, the reservoir data includes the institute on the stratum State permeability, thickness and gesture；

(b) by by it is in the reservoir model, therebetween with horizontal fluid flow and between flow obstacle layer The data group of the horizontal fluid flow layer of horizontal fluid flow layer become single level stream in the matrix Dynamic layer, to form simplified well model system matrix；

(c) bottom hole pressure of the horizontal well is determined；

(d) coupling reservoir well model is formed comprising the full calculating matrix reservoir model and the simplified well model system square Horizontal well processing is the specified well of base apertures, pressure with determining bottom hole pressure, the coupling reservoir well by battle array Model is with a matrix type：

Wherein, A_RRIt is the tridiagonal matrix of the reservoir data, A_RWIt is the life on the stratum adjacent with the perforated interval The vector of yield index；A_WRIt is the vector of the productivity index on the stratum from well to the reservoir；A_WWIt is the horizontal well Productivity index linear scalar；It is the vector of the unknown reservoir gesture of the unit around the horizontal well；It is the vector of the unknown well gesture in the well bore；It is the reservoir of the reservoir units around the horizontal well The vector of data constant；AndIt is the vector of the well data constant of the horizontal well；

(e) solve the coupling reservoir well model, obtain the fluid in the reservoir units on the stratum flow and With the productivity and gesture of the reservoir units of the time step in each of the stratum；

(f) the coupling reservoir well model is solved, is obtained with the time step at the perforated interval of the reservoir The well unit the productivity index；

(g) based on the reservoir units of the time step at the perforated interval of the reservoir and well unit The productivity index of the determination, come determine the horizontal well the horizontal fluid flow layer and the flow obstacle layer institute State the layer completion rate of composition fluid；

(h) according to the institute of the horizontal fluid flow layer and the flow obstacle layer of the horizontal well of the time step The layer completion rate for stating the determination of composition fluid, to determine total well productivity of the horizontal well；And

(i) it is formed with the institute of the horizontal fluid flow layer and the flow obstacle layer of the horizontal well of the time step State the record of the layer completion rate of the determination of composition fluid and total well productivity of the determination of the horizontal well.