A Flow Simulation and Transient Well Analysis Method Based on Generalized Tube Flow and Percolation Coupling

Abstract

This invention discloses a multi-phase flow simulation analysis method based on generalized mobility, which comprises the following steps: S1: The generalized mobility describes fluid flow laws in different subset of study area by using the generalized mobility models with the same form; S2: On the basis of generalized mobility, the multi-component multi-phase flow simulation equations are established. Through solving the above mentioned multi-component multi-phase flow simulation equations, the pressure, temperature, saturation, and mole percentage of each component and each phase of multicomponent multiphase flow fluids in study area are obtained; S3: The corresponding application software are formed by using the established multi-component multi-phase flow simulation and analysis equations. The invention plays an important role in solving the single and multi-well flow simulation of complex multicomponent multiphase flow reservoirs, multi-well interference analysis, deliverability analysis, transient pressure analysis, transient rate analysis, transient temperature analysis, well test design, and permanent downhole monitoring data analysis.

Claims

1. A flow simulation and transient well analysis method based on generalized tube flow and percolation coupling, characterized in comprising the following steps: S1: Based on defining the generalized mobility, fluid flow laws in different subset of study area are characterized by using the generalized mobility models with the same form; S2: On the basis of generalized mobility, the multi-component multi-phase flow governing equations are established by considering convection term, diffusion term, accumulation term, adsorption term, and source/sink term. Then a whole set of multi-component multi-phase flow simulation equations are formed through combining energy conservation equation, auxiliary equations with saturation and capillary pressure, initial saturation equation, initial pressure equation, initial temperature equation, phase equilibrium equation, and boundary condition equations; The pressure, temperature and saturation of multi-phase fluid as well as the mole percentage of each component in each phase at any point in the study areas are obtained by solving the above-mentioned multi-component multi-phase flow simulation equations; S3: The corresponding application software are formed by using the established multi-component multi-phase flow simulation and analysis equations.

2. The flow simulation and transient well analysis method based on generalized tube flow and percolation coupling as claimed according to claim 1, characterized in that the details of step S1 is as follows: S11: For any type of fluid motion equation, if it can be written as v=−λ∇p, λ is called generalized mobility. where v denotes the fluid flow velocity, ∇p denotes the pressure gradient, and the generalized mobility λ is a function of space position and time; S12 The general mobility models with the same form are used to characterize the flow laws of fluid in the tube flow area of wellbores, pipes, fractures, vugs, holes, cavities, caves, fracture caves, karst caves and caverns, and the percolation area of porous media.

3. The flow simulation and transient well analysis method based on generalized tube flow and percolation coupling as claimed according to claim 1, characterized in that the details of step S2 is as follows: S21: The generalized mobility of three-dimensional multi-phase flow ${\lambda }_k=\left[\begin{array}{ccc}{\lambda }_{k,\mathrm{xx}}\left(x,t\right)& {\lambda }_{k,\mathrm{xy}}\left(x,t\right)& {\lambda }_{k,\mathrm{xz}}\left(x,t\right)\\ {\lambda }_{k,\mathrm{yx}}\left(x,t\right)& {\lambda }_{k,\mathrm{yy}}\left(x,t\right)& {\lambda }_{k,\mathrm{yz}}\left(x,t\right)\\ {\lambda }_{k,\mathrm{zx}}\left(x,t\right)& {\lambda }_{k,\mathrm{zy}}\left(x,t\right)& {\lambda }_{k,\mathrm{zz}}\left(x,t\right)\end{array}\right],$ where k=1, 2, . . . , n.sub.p, n.sub.p denotes the total phase number, x denotes space position, t denotes time, 9 components of generalized mobility λ.sub.k,xx, λ.sub.k,xy, λ.sub.k,xz, λ.sub.k,yx, λ.sub.k,yy, λ.sub.k,yz, λ.sub.k,zx, λ.sub.k,zy, λ.sub.k,zz are all functions of space position x and time t; S22: Rewriting the three-dimensional multi-phase flow motion equation into a standard form v.sub.k=−λ.sub.k∇p.sub.k by applying three-dimensional multi-phase flow generalized mobility, where k=1, 2, . . . , n.sub.p, n.sub.p denotes the total phase number, $\nabla ={\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y},\frac{\partial }{\partial z}\right)}^T$ is Hamilton operator, v.sub.k is fluid velocity of k-phase, p.sub.k is pressure of k-phase; S23: Multi-component multi-phase flow simulation equations The multi-component multi-phase flow equations considering convection term, diffusion term, accumulation term, adsorption term, and source/sink term are written as: Component governing equation: $\left[\nabla .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.\left({\rho }_k.Math.C_{\mathrm{ik}}.Math.{\lambda }_k.Math.\nabla p_k\right)\right]+\left[\nabla .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.\left({\mathrm{\phi \rho }}_k.Math.S_k.Math.D_{\mathrm{ik}}.Math.\nabla C_{\mathrm{ik}}\right)\right]=\frac{\partial }{\partial t}\left[\phi .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.\left({\rho }_k.Math.S_k.Math.C_{\mathrm{ik}}\right)\right]+\frac{\partial }{\partial t}\left[\left(1-\phi \right).Math.{\rho }_r.Math.V_i\right]+\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.C_{\mathrm{ik}}.Math.q_k,.Math..Math.\left(x,t\right)\in \Omega \times \left(0,t_{\max }\right],i=1,2,...,n_c$ Energy conservation equation: $\left[\nabla .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math.\left({}_k.Math.{\rho }_k.Math.{\lambda }_k.Math.\nabla p_k\right)\right]+\left[\nabla .Math.\kappa .Math.\nabla T\right]=\frac{\partial }{\partial t}\left[\phi .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math.\left({\rho }_k.Math.S_k.Math.U_k\right)\right]+\frac{\partial }{\partial t}\left[\left(1-\phi \right).Math.{\rho }_t.Math.U_r\right]+.Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.{}_k.Math.q_k,.Math..Math.\left(x,t\right)\in \Omega \times \left(0,t_{\max }\right]$ Auxiliary equations with saturation and capillary pressure: Saturation equation $\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.S_k=1,\left(x,t\right)\in \Omega \times \left(0,t_{\max }\right]$ Capillary pressure equation at α-β phase interface
p.sub.cαβ(S.sub.w)=p.sub.α-p.sub.β,(x,t)∈Ω×(0,T],α=1, . . . ,n.sub.p;β=1, . . . ,n.sub.p Phase equilibrium equation: The phase equilibrium constant of component i in α and β phase:
K.sub.iαβ=C.sub.iα/C.sub.iβ,α=1, . . . ,n.sub.p;β=1, . . . ,n.sub.p;i=1, . . . ,n.sub.c Mole percentage normalization condition of components in each phase: $\underset{i=1}{\overset{n_c}{.Math.}}.Math.C_{\mathrm{ik}}=1,k=1,...,.Math.n_p$ The total mole percent of component i: $Z_i=\stackrel{n_p}{\underset{k=1}{.Math.}}.Math.C_{\mathrm{ik}}.Math.,i=1,...,n_c$ Normalization conditions for ratio of moles of each phase to the whole system: $\underset{i=1}{\overset{n_c}{.Math.}}.Math.=1,k=1,...,n_p$ Boundary condition equations: Boundary condition equation of pressure for each phase $\left(c_{k,1}.Math.p_k+c_{k,2}.Math.{\lambda }_k.Math.\frac{\partial p_k}{\partial n^{\partial \Omega }}\right)=g_k\left(x,t\right),\left(x,t\right)\in \partial \Omega \times \left(0,t_{\max }\right],.Math..Math.k=1,...,n_p$ Boundary condition equation of temperature $\left(d_1.Math.T+d_2.Math.\kappa .Math.\frac{\partial T}{\partial n^{\partial \Omega }}\right)=w\left(x,t\right),\left(x,t\right)\in \partial \Omega \times \left(0,t_{\max }\right]$ Initial saturation equation: Initial saturation equation of each phase
S.sub.k(x,0)=l.sub.k(x),x∈Ω,k=1, . . . ,n.sub.p Initial pressure equation: Initial pressure equation of each phase
p.sub.k(x,0)=ƒ.sub.k(x),x∈Ω,k=1, . . . ,n.sub.p Initial temperature equation:
T(x,0)=τ(x),x∈Ω Variable symbols description in multi-component multi-phase flow simulation analysis equations: ϕ is porosity, which is the function of average pressure p, %; n.sub.c is the total component number, dimensionless; n.sub.p is the total phase number, dimensionless; λ.sub.k is the generalized mobility of k-phase, m.sup.2/(Pa.Math.s); ρk, ρ.sub.r are densities of k-phase and rock, kg/m.sup.3; π.sub.k is enthalpy of k phase, J/kg; U.sub.k, U.sub.r are internal energy of k-phase and rock, J/kg; V.sub.i is adsorption amount of component i, dimensionless; S.sub.k is the k-phase saturation; t is time, s; t.sub.max is the maximum of time, s; p.sub.k is the k-phase pressure, Pa; l.sub.k is the k-phase initial saturation distribution function; q.sub.k is the source/sink term of k-phase, kg/(m.sup.3.Math.s); ƒ.sub.k is the k-phase initial pressure distribution function, Pa; K.sub.iαβ is the phase equilibrium constant of component i in α and δ phase, dimensionless; C.sub.ik is the mole percent of component i in k-phase, dimensionless; D.sub.ik is the diffusion coefficient of component i in k-phase, m.sup.2/s; Z.sub.i is the total mole percent of component i, dimensionless; is the ratio of k-phase to the whole system, dimensionless; g.sub.k is the boundary functions of k-phase on reservoir boundary, dimensionless; w is boundary function of temperature on reservoir boundary, dimensionless; τ is the initial temperature distribution function of reservoirs, K; pcαβ is the capillary pressure at α-β phase interface, Pa; Ω is reservoir space; ∂Ω is reservoir boundary including internal boundary and outer boundary; c.sub.k,1 is the pressure term coefficients of k-phase on reservoir boundary, 1/Pa; c.sub.k,2 is the k-phase coefficients of derivative terms along the outer normal direction on reservoir boundary condition, s/m; d.sub.1 is the temperature term coefficient on reservoir boundary, 1/K; d.sub.2 is temperature coefficients of derivative terms along the outer normal direction on reservoir boundary condition, s.Math.m.sup.2/J; κ is coefficient of the thermal conductivity, J/(m.Math.s.Math.K); n.sup.∂Ω is the outer normal direction on reservoir boundary, m; ∇ is the Hamilton operator; ∂ is the partial derivative sign; S24: The analytical solution algorithms of above multi-component multi-phase flow simulation equations include direct solving method, Laplace transformation, Fourier transformation, and orthogonal transformation method Their numerical methods include finite difference method, finite volume method, boundary element method, and finite element method; After solving, the pressure, temperature, and saturation in multi-phase flow are obtained, including the pressure, temperature, saturation, and mole percent of each component in each phase at any time and at any position in the study area.

4. The flow simulation and transient well analysis method based on generalized tube flow and percolation coupling as claimed according to claim 1, characterized in that the details of step S3 is as follows: S31: The application software described above includes 5 main parts: data pre-processing system, numerical simulation system, analytical analysis system, and analysis results output system, data input and output management system. Its analysis process includes reservoir definition, setting of initial and boundary conditions, wellbore and fracture setting, numerical simulator selection or fluid type and composition setting, generalized mobility model definition, grid design of numerical simulation, wellbore storage model setting, flow period and regime definition, coordinate transformation, setting of models and their type curve analysis, parameter adjustment and history matching, dynamic prediction; S32: The application software described above can be used for single and multi-well flow simulation of complex multi-component multi-phase flow reservoirs, multi-well interference analysis, deliverability analysis, transient pressure analysis, transient rate analysis, transient temperature analysis, well test design, and permanent downhole monitoring data analysis. The complex multi-component multi-phase flow reservoirs mentioned above include all types of fluid reservoirs, oil and gas reservoirs with fluid injection, underground gas storage, ground water reservoirs, and geothermal reservoirs.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0055] FIG. 1 shows 18 types of basic structural units of a reservoir, which is provided by this invention. A complex reservoir can be formed through nesting, stacking and splicing of a plurality of basic reservoir structural units.

[0056] FIG. 2 shows 32 types of basic reservoir bodies, which is provided by this invention. One basic reservoir structural unit can be selected from one of 32 types of basic reservoirs. A complex reservoir can be formed through nesting, stacking and splicing of a plurality of basic reservoir structural units. Every basic reservoir structural unit can be selected from one of 32 types of basic reservoirs.

[0057] FIG. 3 is a schematic diagram of a complex reservoir model for tube flow and percolation coupling, which is provided in embodiment 1 of the invention. The whole reservoir space is consisted of four reservoir bodies and three different wells (wells may also be considered as a part of the reservoirs). Four reservoir bodies include two low-permeability reservoirs, one karst cave reservoir, and one high permeability reservoir; three wells include one inclined well, one vertical well, and one multi-segment fractured horizontal well. This is merely a kind of special case of the example.

[0058] FIG. 4 is an auxiliary plot of the product of the oil phase generalized mobility and the normal vector at the internal discrete control point, which is provided in embodiment 1 of the invention.

[0059] FIG. 5 is an auxiliary plot of the product of the water phase generalized mobility and the normal vector at the internal discrete control point, which is provided in embodiment 1 of the invention.

[0060] FIG. 6 is an auxiliary plot of the product of the gas phase generalized mobility and the normal vector at the internal discrete control point, which is provided in embodiment 1 of the invention.

[0061] FIG. 7 is an auxiliary plot of the product of the oil phase generalized mobility and the normal vector at the boundary discrete control point, which is provided in embodiment 1 of the invention.

[0062] FIG. 8 is an auxiliary plot of the product of the water phase generalized mobility and the normal vector at the boundary discrete control point, which is provided in embodiment 1 of the invention.

[0063] FIG. 9 is an auxiliary plot of the product of the gas phase generalized mobility and the normal vector at the boundary discrete control point, which is provided in embodiment 1 of the invention.

[0064] FIG. 10 is a schematic diagram of black oil model of oil/gas/water three-phase flow through treating the wellbore as inner boundary for a partially open vertical well located in a homogeneous reservoir, which is provided in embodiment 2 of the invention.

[0065] FIG. 11 is a schematic diagram of black oil model of oil/gas/water three-phase flow through treating the wellbore as reservoir for a partially open vertical well located in a homogeneous reservoir, which is provided in embodiment 3 of the invention.

[0066] FIG. 12 is a schematic diagram of oil phase flow model for a partially open vertical well located in a homogeneous reservoir through treating the wellbore as inner boundary, which is provided in embodiment 4 of the invention.

[0067] FIG. 13 is a schematic diagram of oil phase flow model for a partially open vertical well located in a homogeneous reservoir through treating the wellbore as reservoir, which is provided in embodiment 5 of the invention.

[0068] FIG. 14 is a schematic diagram of oil phase flow model for a fully open vertical well located in a homogeneous radial 2-zone composite reservoir through treating the wellbore as inner boundary, which is provided in embodiment 6 of the invention.

[0069] FIG. 15 is a schematic diagram of oil phase flow model for a fully open vertical well located in a homogeneous radial 2-zone composite reservoir through treating the wellbore as reservoir, which is provided in embodiment 7 of the invention.

[0070] FIG. 16 is a schematic diagram of a tube-shaped reservoir model, which is provided in embodiment 8 of the invention.

[0071] FIG. 17 is the typical pressure draw-down type curves of a tube-shaped reservoir model, which is provided in embodiment 8 of the invention. The typical characteristics include: {circle around (1)} early-stage of pressure derivative curve is parallel to the pressure difference curve with a ½ slope (linear flow); {circle around (2)} late-stage pressure derivative curve has a 1 slope (pseudo-steady flow).

[0072] FIG. 18 is the pressure build-up type curve corresponding to FIG. 17.

[0073] FIG. 19 is the typical pressure draw-down type curves corresponding to FIG. 17 with consideration of skin and wellbore storage effect. The early-stage pressure derivative curve overlaps with the pressure difference curves with a slope of 1. The pressure derivative curve then goes upward due to (positive) skin action.

[0074] FIG. 20 is a schematic diagram of a cylindrical reservoir model, which is provided in embodiment 9 of the invention.

[0075] FIG. 21 is the typical pressure draw-down type curves of a cylindrical reservoir model, which is provided in embodiment 9 of the invention. The typical characteristics include: {circle around (1)} early-stage of pressure derivative curve has a −½ slope (spherical flow); {circle around (2)} middle-stage pressure derivative curve has a 0 slope (radial flow); {circle around (3)} late-stage pressure derivative curve has a 1 slope (pseudo-steady flow).

[0076] FIG. 22 is the pressure build-up type curve corresponding to FIG. 21.

[0077] FIG. 23 is the typical pressure draw-down type curves corresponding to FIG. 21 with consideration of skin and wellbore storage effect. The early-stage pressure derivative curve overlaps with the pressure difference curves with a slope of 1. The pressure derivative curve then goes upward due to (positive) skin action.

[0078] FIG. 24 is the typical pressure draw-down type curves of a cylindrical reservoir model, which is provided in embodiment 9 of the invention. The typical characteristics include: {circle around (1)} early-stage of pressure derivative curve has a −½ slope (spherical flow); {circle around (2)} middle-stage pressure derivative curve has a ½ slope (linear flow); {circle around (3)} late-stage pressure derivative curve has a 1 slope (pseudo-steady flow).

[0079] FIG. 25 is the pressure build-up type curve corresponding to FIG. 24.

[0080] FIG. 26 is the typical pressure draw-down type curves corresponding to FIG. 24 with consideration of skin and wellbore storage effect. The early-stage pressure derivative curve overlaps with the pressure difference curves with a slope of 1. The pressure derivative curve then goes upward due to (positive) skin action.

[0081] FIG. 27 is a schematic diagram of a spherical reservoir model, which is provided in embodiment 10 of the invention.

[0082] FIG. 28 is the typical pressure draw-down type curves of a spherical reservoir model, which is provided in embodiment 10 of the invention. The typical characteristics include: (Dearly-stage of pressure derivative curve has a −½ slope (spherical flow); late-stage pressure derivative curve has a 1 slope (pseudo-steady flow).

[0083] FIG. 29 is the pressure build-up type curve corresponding to FIG. 28.

[0084] FIG. 30 is the typical pressure draw-down type curves corresponding to FIG. 28 with consideration of skin and wellbore storage effect. The early-stage pressure derivative curve overlaps with the pressure difference curves with a slope of 1. The pressure derivative curve then goes upward due to (positive) skin action.

[0085] FIG. 31 is a schematic diagram of a composite model of hollow cylinder and cylindrical reservoir, which is provided in embodiment 11 of this invention.

[0086] FIG. 32 is the typical pressure draw-down type curves of a composite model of hollow cylinder and cylindrical reservoir, which is provided in embodiment 11 of the invention. The typical characteristics is that pressure derivative curve successively shows a −½ slope (spherical flow), 0 slope (radial flow of cylindrical reservoir), the junction of cylindrical reservoir and hollow cylindrical reservoir goes upward (transition regime), 0 slope (radial flow of hollow cylindrical reservoir), 1 slope (pseudo-steady flow).

[0087] FIG. 33 is the pressure build-up type curve corresponding to FIG. 32.

[0088] FIG. 34 is the typical pressure draw-down type curves corresponding to FIG. 32 with consideration of skin and wellbore storage effect. The early-stage pressure derivative curve overlaps with the pressure difference curves with a slope of 1. The pressure derivative curve then goes upward due to (positive) skin action.

[0089] FIG. 35 is the typical pressure draw-down type curves of a hollow cylinder and cylindrical reservoir composite model, which is provided in embodiment 11 of the invention. The typical characteristics is that pressure derivative curve successively shows a −½ slope (spherical flow), 0 slope (radial flow of cylindrical reservoir), the junction of cylindrical reservoir and hollow cylindrical reservoir goes downward with a −½ slope (transient period), 0 slope (radial flow of hollow cylindrical reservoir), 1 slope (pseudo-steady flow).

[0090] FIG. 36 is the pressure build-up type curve corresponding to FIG. 35.

[0091] FIG. 37 is the typical pressure draw-down type curves corresponding to FIG. 35 with consideration of skin and wellbore storage effect. The early-stage pressure derivative curve overlaps with the pressure difference curves with a slope of 1. The pressure derivative curve then goes upward due to (positive) skin action.

[0092] FIG. 38 is a schematic diagram of a composite model of tube-shaped and spherical reservoir, which is provided in embodiment 12 of the invention.

[0093] FIG. 39 is the typical pressure draw-down type curves of tube-shaped and spherical reservoir composite model, which is provided in embodiment 12 of the invention. The typical characteristics include: {circle around (1)} early-stage of pressure derivative curve is parallel to the pressure difference curve with a ½ slope (linear flow); {circle around (2)} middle-stage pressure derivative curve goes downward; {circle around (3)} late-stage pressure derivative curve has a 1 slope (pseudo-steady flow).

[0094] FIG. 40 is the pressure build-up type curve corresponding to FIG. 39.

[0095] FIG. 41 is the typical pressure draw-down type curves corresponding to FIG. 39 with consideration of skin and wellbore storage effect. The early-stage pressure derivative curve overlaps with the pressure difference curves with a slope of 1. The pressure derivative curve then goes upward due to (positive) skin action.

[0096] FIG. 42 is the vertical well Blasingame type curves based on a cylindrical reservoir model, which is provided in embodiment 9 of the invention.

[0097] FIG. 43 is the vertical well Agarwal-Gardner type curves based on a cylindrical reservoir model, which is provided in embodiment 9 of the invention.

[0098] FIG. 44 is the vertical well NP1 type curves based on a cylindrical reservoir model, which is provided in embodiment 9 of the invention.

[0099] FIG. 45 is the vertical well Transient type curves based on a cylindrical reservoir model, which is provided in embodiment 9 of the invention.

[0100] FIG. 46 is the schematic diagram of a one-dimensional tube-shaped reservoir model, which is provided in embodiment 13 [high-speed non-Darcy flow regime I (0≤x≤X.sub.1)+Darcy flow regime II (X.sub.1≤x≤X.sub.2)].

[0101] FIG. 47 is the generalized mobility type curves of Darcy flow regime II of a one-dimensional tube-shaped reservoir model, which is provided in embodiment 13 [high-speed non-Darcy flow regime I (0≤x≤X.sub.1)+Darcy flow regime II (X.sub.1≤x≤X.sub.2)].

[0102] FIG. 48 is the generalized mobility type curves of Darcy flow regime II of a one-dimensional tube-shaped reservoir model, which is provided in embodiment 13 [β=0 in high-speed non-Darcy flow regime I becomes Darcy flow regime I (0≤x≤X.sub.1)+Darcy flow regime II(X.sub.1≤x≤X.sub.2)].

[0103] FIG. 49 is the comparison plot between high-speed non-Darcy flow and Darcy flow of a one-dimensional tube-shaped reservoir model, which is provided in embodiment 13 (it becomes Darcy flow when β=0 in high-speed non-Darcy flow).

[0104] FIG. 50 is a generalized mobility variation type curve of transient temperature analysis model with oil single-phase flow, which is provided by embodiment 18: The larger the generalized mobility, the lower the position of temperature difference and temperature difference derivative curve.

DETAILED DESCRIPTION OF EMBODIMENTS

[0105] The name of the invention is a flow simulation and transient well analysis method based on generalized tube flow and percolation coupling.

[0106] Through comparing the Hagen-Poiseuille (1839,1840) formula of laminar flow and tube flow with the Darcy formula of conventional percolation, it is found that there is no essential difference in mathematical expression form between the two formulas. They both are linear functions with positive proportion i.e. The average velocity of fluid is proportional to the pressure gradient. Fluid flow laws described by both laminar and tube flow formula and Darcy's conventional percolation formula is linear flow. However, Fluid flow laws, which is described by turbulent tube flow equation, Forchheimer (1901) binomial high-speed non-Darcy formula, power-law non-Newton formula, low-speed non-Darcy percolation equation, free flow Navier-Stokes equation, etc, is a nonlinear flow. In order to unify the linear flow and nonlinear flow and consider the anisotropy of the reservoir or others, this invention proposed that the fluid flow velocity of underground fluid flow motion equations are uniformly written as the product of the generalized mobility and the pressure gradient of the negative direction, which means

v.sub.k=−λ.sub.k∇p.sub.k

[0107] where

[00011]${\lambda }_k=\left[\begin{array}{cccc}{\lambda }_{k,\mathrm{xx}}\left(x,t\right)& & {\lambda }_{k,\mathrm{xy}}\left(x,t\right)& {\lambda }_{k,\mathrm{xz}}\left(x,t\right)\\ {\lambda }_{k,\mathrm{yx}}\left(x,t\right)& & {\lambda }_{k,\mathrm{yy}}\left(x,t\right)& {\lambda }_{k,\mathrm{yz}}\left(x,t\right)\\ {\lambda }_{k,\mathrm{zx}}\left(x,t\right)& & {\lambda }_{k,\mathrm{zy}}\left(x,t\right)& {\lambda }_{k,\mathrm{zz}}\left(x,t\right)\end{array}\right]$

is the generalized mobility of three-dimensional multi-phase flow [unit: m2/(Pa.Math.s)]. where k=1, 2, . . . , n.sub.p, n.sub.p denotes the total phase number, x denotes space position, t denotes time, 9 components of generalized mobility λ.sub.k,xx, λ.sub.k,xy, λ.sub.k,xz, λ.sub.k,yx, λ.sub.k,yy, λ.sub.k,yz, λ.sub.k,zx, λ.sub.k,zy, λ.sub.k,zz are all functions of space position x and time t.

[00012]$\nabla ={\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y},\frac{\partial }{\partial z}\right)}^T$

is Hamilton operator (unit: 1/m), v.sub.k is fluid velocity (m/s), p.sub.k is pressure (Pa).

[0108] Because the generalized mobility defined in this invention itself covers the concepts of generalized permeability or equivalent permeability, equivalent percolation coefficient or converted percolation coefficient, when solving some relatively simple practical problems, the use of these definitions can also achieve the similar effect of generalized mobility. If the parameters (such as fluid flow coefficient or transmissibility, pressure transmitting coefficient, etc.) related to fluid flow capacity and reservoir conductivity are defined by using the definition method of this invention, the application scope of the relevant parameters can be extended.

[0109] The invention can be applied to fluid types including multi-component oil/gas/water single-phase flow and multi-component multi-phase flow. The fluid types can be Newton and non-Newtonian. The flowing laws can be linear, non-linear, or locally linear and locally nonlinear.

[0110] The further detailed descriptions of this invention are given in following by combining with drawings. The following embodiments are described only for the purpose of illustration and not for any limitation to the scope of the present disclosure.

[0111] A flow simulation and transient well analysis method based on generalized tube flow and percolation coupling includes:

[0112] S1: Based on defining the generalized mobility, fluid flow laws in different subset of study area are characterized by using the generalized mobility models with the same form;

[0113] S2: On the basis of generalized mobility, the multi-component multi-phase flow governing equations are established by considering convection term, diffusion term, accumulation term, adsorption term, and source/sink term. Then a whole set of multi-component multi-phase flow simulation equations are created through combining energy conservation equation, auxiliary equations with saturation and capillary pressure, initial saturation equation, initial pressure equation, initial temperature equation, phase equilibrium equation, and boundary condition equations; The pressure, temperature and saturation of multi-phase fluid as well as the mole percentage of each component in each phase at any point in the study areas are obtained by solving the above multi-component multi-phase flow simulation equation;

[0114] S3: The corresponding application software are formed by using the established multi-component multi-phase flow simulation and analysis equations.

[0115] wherein S1 (Based on defining the generalized mobility, fluid flow laws in different subset of study area are characterized by using the same formed generalized mobility models) includes:

(1) For any type of fluid motion equation, if it can be written as v=−λ∇p, λ is called generalized mobility, where v denotes the fluid flow velocity, ∇p denotes the pressure gradient, and the generalized mobility λ is a function of space position and time.
(2) The general mobility models with the same form are used to characterize the flow laws of fluid in the tube flow area of wellbores, pipes, fractures, vugs, holes, cavities, caves, fracture caves, karst caves and caverns, and the percolation area of porous media.

[0116] wherein S2 (On the basis of generalized mobility, the multi-component multi-phase flow governing equations are established by considering convection term, diffusion term, accumulation term, adsorption term, and source/sink term. Then a whole set of multi-component multi-phase flow simulation equations are created through combining energy conservation equation, auxiliary equations with saturation and capillary pressure, initial saturation equation, initial pressure equation, initial temperature equation, phase equilibrium equation, and boundary condition equations; The pressure, temperature and saturation of multi-phase fluid as well as the mole percentage of each component in each phase at any point in the study areas are obtained by solving the above multi-component multi-phase flow simulation equation) includes:

(1) The generalized mobility of three-dimensional multi-phase flow

[00013]${\lambda }_k=\left[\begin{array}{cccc}{\lambda }_{k,\mathrm{xx}}\left(x,t\right)& & {\lambda }_{k,\mathrm{xy}}\left(x,t\right)& {\lambda }_{k,\mathrm{xz}}\left(x,t\right)\\ {\lambda }_{k,\mathrm{yx}}\left(x,t\right)& & {\lambda }_{k,\mathrm{yy}}\left(x,t\right)& {\lambda }_{k,\mathrm{yz}}\left(x,t\right)\\ {\lambda }_{k,\mathrm{zx}}\left(x,t\right)& & {\lambda }_{k,\mathrm{zy}}\left(x,t\right)& {\lambda }_{k,\mathrm{zz}}\left(x,t\right)\end{array}\right],$

where k=1, 2, . . . , n.sub.p, n.sub.p denotes the total phase number, x denotes space position, t denotes time, 9 components of generalized mobility λ.sub.k,xx, λ.sub.k,xy, λ.sub.k,xz, λ.sub.k,yx, λ.sub.k,yy, λ.sub.k,yz, λ.sub.k,zx, λ.sub.k,zy, λ.sub.k,zz are all functions of space position x and time t.

[0117] From the generalized mobility of three-dimensional multi-phase flow, it is easy to obtain the generalized mobility of two-dimensional multi-phase flow

[00014]${\lambda }_k=\left[\begin{array}{cc}{\lambda }_{k,\mathrm{xx}}\left(x,t\right)& {\lambda }_{k,\mathrm{xy}}\left(x,t\right)\\ {\lambda }_{k,\mathrm{yx}}\left(x,t\right)& {\lambda }_{k,\mathrm{yy}}\left(x,t\right)\end{array}\right],$

where k=1, 2, . . . , n.sub.p, n.sub.p denotes the total phase number, x denotes space position, t denotes time, 4 components of generalized mobility λ.sub.k,xx, λ.sub.k,xy, λ.sub.k,yx, λ.sub.k,yy, are all functions of space position x and time t.

[0118] From the generalized mobility of three-dimensional multi-phase flow, it is easy to obtain the generalized mobility of one-dimensional multi-phase flow λ.sub.k=λ.sub.k,xx(x,t), where k=1, 2, . . . , n.sub.p, n.sub.p denotes the total phase number, x denotes space position, t denotes time.

[0119] Rewriting the three-dimensional multi-phase flow motion equation into a standard form v.sub.k=−λ.sub.k∇p.sub.k by applying three-dimensional multi-phase flow generalized mobility, where k=1, 2, . . . , n.sub.p, n.sub.p denotes the total phase number,

[00015]$\nabla ={\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y},\frac{\partial }{\partial z}\right)}^T$

is Hamilton operator, v.sub.k is fluid velocity of k-phase, p.sub.k is pressure of k-phase.

[0120] Rewriting the two-dimensional multi-phase flow motion equation into a standard form v.sub.k=−λ.sub.k∇p.sub.k by applying two-dimensional multi-phase flow generalized mobility, where k=1, 2, . . . , n.sub.p, n.sub.p denotes the total phase number.

[00016]$\nabla ={\left(\frac{\partial }{\partial x},\frac{\partial }{\partial y}\right)}^T$

is Hamilton operator, v.sub.k is fluid velocity of k-phase, p.sub.k is pressure of k-phase.

[0121] Rewriting the one-dimensional multi-phase flow motion equation into a standard form v.sub.k=λ.sub.k∇p.sub.k by applying one-dimensional multi-phase flow generalized mobility, where k=1, 2, . . . , n.sub.p, n.sub.p denotes the total phase number,

[00017]$\nabla =\frac{\partial }{\partial x}$

is Hamilton operator, v.sub.k is fluid velocity of k-phase, p.sub.k is pressure of k-phase.
(3) Multi-component multi-phase flow simulation equation

[0122] The multi-component multi-phase flow equations considering convection term, diffusion term, accumulation term, adsorption term, and source/sink term are written as:

[0123] Component governing equation:

[00018]$\left[\nabla .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.\left({\rho }_k.Math.C_{\mathrm{ik}}.Math.{\lambda }_k.Math.\nabla p_k\right)\right]+\left[\nabla .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.\left({\mathrm{\phi \rho }}_k.Math.S_k.Math.D_{\mathrm{ik}}.Math.\nabla C_{\mathrm{ik}}\right)\right]=.Math.\frac{\partial }{\partial t}\left[\phi .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.\left({\rho }_k.Math.S_k.Math.C_{\mathrm{ik}}\right)\right]+\frac{\partial }{\partial t}\left[\left(1-\phi \right).Math.{\rho }_r.Math.V_i\right]+.Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.C_{\mathrm{ik}}.Math.q_k,\left(x,t\right)\in \Omega \times \left(0,t_{\max }\right],i=1,2,...,.Math.n_c$

[0124] Energy conservation equation:

[00019]$\left[\nabla .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.\left({}_k.Math.{\rho }_k.Math.{\lambda }_k.Math.\nabla p_k\right)\right]+\left[\nabla .Math.\kappa .Math.\nabla T\right]=.Math.\frac{\partial }{\partial t}\left[\phi .Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.\left({\rho }_k.Math.S_k.Math.U_k\right)\right]+\frac{\partial }{\partial t}\left[\left(1-\phi \right).Math.{\rho }_r.Math.U_r\right]+.Math.\underset{k=1}{\overset{n_p}{.Math.}}.Math.{}_k.Math.q_k,\left(x,t\right)\in \Omega \times \left(0,t_{\max }\right]$

[0125] Auxiliary equations with saturation and capillary pressure:

[0126] Saturation equation

[00020]$\underset{k=1}{\overset{n_p}{.Math.}}.Math.S_k=1,\left(x,t\right)\in \Omega \times \left(0,t_{\max }\right]$

[0127] Capillary pressure equation at α-β phase interface

p.sub.cαβ(S.sub.w)=p.sub.α-p.sub.β,(x,t)∈Ω×(0,T],α=1, . . . ,n.sub.p;β.sub.1, . . . ,n.sub.p

[0128] Equilibrium equation:

[0129] The equilibrium constant of component i in α and β phase:

K.sub.iαβ=C.sub.iα/C.sub.iβ,α=1, . . . ,n.sub.p;β=1, . . . ,n.sub.p;i=1, . . . ,n.sub.c

[0130] Mole percentage normalization condition of components in each phase:

[00021]$\underset{i=1}{\overset{n_c}{.Math.}}.Math..Math.C_{\mathrm{ik}}=1,k=1,...,n_p$

[0131] The total mole percent of component i:

[00022]$Z_i=\underset{k=1}{\overset{n_p}{.Math.}}.Math..Math.C_{\mathrm{ik}}.Math.,i=1,...,n_c$

[0132] Normalizing conditions for ratio of moles of each phase to the whole system:

[00023]$\underset{i=1}{\overset{n_c}{.Math.}}.Math..Math.=1,k=1,...,n_p$

[0133] Boundary condition equation:

[0134] Boundary condition equation of pressure for each phase

[00024]$\left(c_{k,1}.Math.p_k+c_{k,2}.Math.{\lambda }_k.Math.\frac{\partial p_k}{\partial n^{\partial \Omega }}\right)=g_k\left(x,t\right),\left(x,t\right)\in \partial \Omega \times \left(0,t_{\max }\right],k=1,...,n_p.Math.\left(d_1.Math.T+d_2.Math.\kappa .Math.\frac{\partial T}{\partial n^{\partial \Omega }}\right)=w\left(x,t\right),\left(x,t\right)\in \partial \Omega \times \left(0,t_{\max }\right]$

[0135] Initial saturation equation:

[0136] Initial saturation equation of each phase

S.sub.k(x,0)=l.sub.k(x),x∈Ω,k=1, . . . ,n.sub.p

[0137] Initial pressure equation:

[0138] Initial pressure equation of each phase:

p.sub.k(x,0)=ƒ.sub.k(x),x∈Ω,k=1, . . . ,n.sub.p

[0139] Initial temperature equation:

T(x,0)=τ(x),x∈Ω

[0140] Variable symbols description in multi-component multi-phase flow simulation analysis equations: ϕ is porosity, which is the function of average pressure p, %; n.sub.c is the total component number, dimensionless; n.sub.p is the total phase number, dimensionless; λ.sub.k is the generalized mobility of k-phase, m.sup.2/(Pa.Math.s); ρ.sub.k, ρ.sub.r are densities of k-phase and rock, kg/m.sup.3; π.sub.k is enthalpy of k phase, J/kg; U.sub.k, U.sub.r are internal energy of k-phase and rock, J/kg; Vi is adsorption amount of component i, dimensionless; Sk is the k-phase saturation; t is time, s; t.sub.max is the maximum of time, s; p.sub.k is the k-phase pressure, Pa; l.sub.k is the k-phase initial saturation distribution function q.sub.k is the source/sink term of k-phase, kg/(m.sup.3.Math.s); ƒ.sub.k is the k-phase initial pressure distribution function, Pa; K.sub.iαβ is the phase equilibrium constant of component i in α and β phase, dimensionless; C.sub.ik is the mole percent of component i in k-phase, dimensionless; D.sub.ik is the diffusion coefficient of component i in k-phase, m.sup.2/s; Z.sub.i is the total mole percent of component i, dimensionless; is the ratio of k-phase to the whole system, dimensionless; g.sub.k is the boundary functions of k-phase on reservoir boundary, dimensionless; w is boundary function of temperature on reservoir boundary, dimensionless; τ is the initial temperature distribution function of reservoirs, K; pee is the capillary pressure at α-β phase interface, Pa; Ω is reservoir space; ∂Ω is reservoir boundary including internal boundary and outer boundary; c.sub.k,1 is the pressure term coefficients of k-phase on reservoir boundary, 1/Pa c.sub.k,2 is the k-phase coefficients of derivative terms along the outer normal direction on reservoir boundary condition, s/m; d.sub.1 is the temperature term coefficient on reservoir boundary, 1/K; d.sub.2 is temperature coefficients of derivative terms along the outer normal direction on reservoir boundary condition, s.Math.m.sup.2/J; κ is coefficient of the thermal conductivity, J/(m.Math.s.Math.K); n.sup.∂Ω is the outer normal direction on reservoir boundary, m; ∇ is the Hamilton operator; ∂ is the partial derivative sign;

(4) The analytical solution of above multi-component multi-phase flow simulation equations include direct solving method, Laplace transformation, Fourier transformation, and orthogonal transformation method; The numerical methods include finite difference method, finite volume method, boundary element method, and finite element method; After solving, the pressure, temperature, and saturation in multi-phase flow are obtained, including the pressure, temperature, saturation, and mole percent of each component in each phase at any time and at any position in the study area.

[0141] Wherein S3 (The corresponding application software are developed on the basis of multi-component multi-phase flow simulation analysis equations) includes:

(1) The application software described above includes 5 main parts: data pre-processing system, numerical simulation system, analytical analysis system, and analysis results output system, data input and output management system, etc. The key technologies include reservoir definition, setting of initial and boundary conditions, wellbore and fracture setting, numerical simulator selection or fluid type and composition setting, generalized mobility model definition (including defining generalized permeability or equivalent permeability, etc., in simplified cases), grid design of numerical simulation, wellbore storage model setting, flow period and regime definition, coordinate transformation, setting of models and their type curve analysis, parameter adjustment and history fitting, dynamic prediction. The details are as follows:

(1.1) Reservoir Definition

[0142] The basic reservoir structural unit is shown in FIG. 1. The reservoir is composed of a single or several basic reservoir structural units, which can be nested, stacked and spliced among each other.

[0143] Reservoir boundaries include constant pressure boundary, variable pressure boundary, sealed boundary, semi-permeable boundary, infinite acting boundary and so on.

[0144] The reservoir types are as shown in FIG. 2. The reservoir types of basic reservoir structural unit include homogeneous, dual porous media, triple porous media, fractal media, cavities caves, caverns, fractures and so on.

(1.2) Setting of Initial and Boundary Conditions

[0145] Initialization of pressure, temperature, saturation and internal and external boundary conditions. The Inner boundary models include effective wellbore diameter model, source function model and so on.

(1.3) Wellbore and Fracture Setting

[0146] Well types include vertical well, horizontal well, inclined well, fractured well, multibranch well and so on.

[0147] The wellbore and fracture types include uniform flux, infinite conductivity, and finite conductivity and so on.

[0148] The wellbore and fracture opening degree includes fully open and partially open.

(1.4) Numerical Simulator Selection or Fluid Type and Composition Setting

[0149] The numerical simulator selection includes black oil model, compositional model and thermal recovery model. The fluid phases (including composition) include oil single-phase, water single-phase, gas single-phase, oil-water two-phase, oil-gas two-phase, gas-water two-phase, oil-gas-water three phase, condensate gas, carbon dioxide, foam fluid, polymer and so on.

(1.5) Generalized Mobility Model Definition (Including Defining Generalized Permeability or Equivalent Permeability, Etc., in Simplified Cases)

[0150] The generalized mobility models include: {circle around (1)} Newton fluid+Darcy flow, {circle around (2)} Newton fluid+Laminar tube flow, {circle around (3)} Newton fluid+High speed Non-Darcy flow, {circle around (4)} Newton fluid+turbulent tube flow, {circle around (5)} Newton fluid+pressure sensitive flow, {circle around (6)} Non-Newton fluid+Darcy flow, {circle around (7)} Non-Newton fluid+Laminar tube flow, {circle around (8)} Non-Newton fluid+pressure sensitive flow, {circle around (9)} Non-Newton fluid+turbulent tube flow, {circle around (10)} Newton fluid+low speed non-Darcy flow, {circle around (11)} Newton fluid+slippage effect, etc.

EXAMPLES

[0151] {circle around (1)} Newton fluid+Darcy flow

[0152] On the basis of Darcy's percolation formula (1856), λ is written as:

[00025]$\lambda =\frac{K}{\mu }$

[0153] Where μ is the viscosity of the fluid, K is the permeability,

[0154] For the form of considering the influence of gravity, λ is written as:

[00026]$\lambda =\frac{K}{\mu }.Math.\left(1-\frac{\rho .Math.g.Math.\cos .Math.\alpha }{.Math.\frac{\partial p}{\partial l}.Math.}\right)$

[0155] Where μ is the viscosity of the fluid, K is the permeability, ρ is the density of the fluid, g is the acceleration of gravity, p is the pressure, l is the distance, α is the angle between pressure gradient direction and gravity direction.

{circle around (2)} Newton fluid+Laminar tube flow

[0156] On the basis of Hagen-poiseuille's formula (1839,1840), λ is written as:

[00027]$\lambda =\frac{d^2}{3.Math.2.Math.\mu }$

[0157] Where μ is the viscosity of the fluid, d is the hydraulic diameter of tubes.

{circle around (3)} Newton fluid+High speed Non-Darcy flow

[0158] On the basis of Forchheimer binomial high-speed Non-darcy formula (1901), λ is written as:

[00028]$\lambda =\frac{1}{\frac{\mu }{K}+\mathrm{\beta \rho }.Math..Math.v.Math.}$

[0159] Where μ is the viscosity of the fluid, K is the permeability, β is the high speed Non-Darcy factor, ρ is the density of the fluid, v is the velocity of the fluid.

{circle around (4)} Newton fluid+turbulent tube flow

[0160] On the basis of Darcy-Weisbach formula (1845) and Colebrook on-way resistance formula (1938), λ is written as:

[00029]$\lambda =0.0232.Math.\frac{d}{\rho .Math..Math.v}.Math.{\left(\frac{d.Math..Math.\rho .Math..Math.v.Math.}{d}.Math.\frac{\mathrm{.Math.}}{d}-18.574.Math.W\left(0.199.Math.\frac{d.Math..Math.\rho .Math..Math.v}{\mu }.Math.e^{0.0538.Math.\frac{d.Math..Math.\rho .Math..Math.v}{\mu }.Math.\frac{\mathrm{.Math.}}{d}}\right)\right)}^2$

[0161] Where μ is the viscosity of the fluid, d is the hydraulic diameter of tubelines, ρ is the density of the fluid, ε/d is the relative roughness, v is the velocity, W(g) is the Lambert W function or product logarithm function.

{circle around (5)} Newton fluid+pressure sensitive flow

[0162] On the basis of permeability modulus (pressure sensitivity) formula, λ is written as

[00030]$\lambda =\frac{K_i.Math.e^{\gamma \left(p-p_i\right)}}{\mu }$

[0163] Where μ is the viscosity of the fluid, K.sub.i is the initial permeability of reservoir, γ is the permeability modulus, P is the pressure, p.sub.i is the initial reservoir pressure.

{circle around (6)} Non-Newton fluid+Darcy flow

[0164] On the basis of the Darcy's percolation formula, Qstwald-DeWaele power-law fluid viscosity formula (1923,1925) and Hirasaki-Pope shear rate formula (1974), λ is written as:

[00031]$\lambda =\left(k/{\mu }_{\mathrm{eff}}\right).Math.{.Math.v.Math.}^{1-n}$${\mu }_{\mathrm{eff}}=\frac{H}{1.Math.2}.Math.{\left(9+\frac{3}{n}\right)}^n.Math.{\left(1.Math.5.Math.0.Math.K.Math.\phi \right)}^{\frac{1-n}{2}}$

[0165] Where μ.sub.eff is the effective viscosity, K is the permeability, n is the power-law index, v is the fluid flow velocity, H is the consistency coefficient, ϕ is the porosity.

{circle around (7)} Non-Newton fluid+Laminar tube flow

[0166] On the basis of Hagen-poiseuille laminar tube flow formula (1839,1840), Qstwald-DeWaele power-law fluid viscosity formula (1923,1925) and Hirasaki-Pope shear rate formula (1974), λ is written as:

[00032]$\lambda =\frac{d^2}{32.Math.{\mu }_{\mathrm{eff}}.Math.v^{n-1}}$

[00033]${\mu }_{e.Math.f.Math.f}=\frac{H}{12}.Math.{\left(9+\frac{3}{n}\right)}^n.Math.{\left(\frac{75.Math.d^2}{1.Math.6}\right)}^{\frac{1-n}{2}}$

[0167] Where μ.sub.eff is the effective viscosity, d is the hydraulic diameter of tubes, n is the power-law index, v is the fluid flow velocity, H is the consistency coefficient.

{circle around (8)} Non-Newton fluid+pressure sensitive flow

[0168] On the basis of the permeability modulus(stress sensitivity) formula, Qstwald-DeWaele power-law fluid viscosity formula (1923,1925) and Hirasaki-Pope shear rate formula (1974), λ is written as:

[00034]$\lambda =\frac{K_i.Math.e^{\gamma \left(p-p_i\right)}}{{\mu }_{\mathrm{eff}}.Math.v^{n-1}}$${\mu }_{\mathrm{eff}}=\frac{H}{1.Math.2}.Math.{\left(9+\frac{3}{n}\right)}^n.Math.{\left(1.Math.5.Math.0.Math.K.Math.\phi \right)}^{\frac{1-n}{2}}$

[0169] Where μ.sub.eff is the effective viscosity, K.sub.i is the initial permeability of reservoir, γ is the permeability modulus, p is the pressure, p.sub.i is the initial reservoir pressure, v is the fluid flow velocity, n is the power-law index, H is the consistency coefficient, K is the permeability, ϕ is the porosity.

{circle around (9)} Non-Newton fluid+turbulent tube flow

[0170] On the basis of Darcy-Weisbach formula (1845), Colebrook resistance formula along the way (1938), Qstwald-DeWaele power-law fluid viscosity formula (1923,1925), and Hirasak-Pope shear rate formula (1974), λ is written as:

[00035]$\lambda =0.0232.Math.\frac{d}{\rho .Math..Math.v}.Math.{\left(\frac{d.Math.\rho }{{\mu }_{\mathrm{eff}}.Math.v^{n-2}}.Math.\frac{\mathrm{.Math.}}{d}-18.574.Math.W\left(0.199.Math..Math.\frac{d.Math.p}{{\mu }_{\mathrm{eff}}.Math.v^{n-2}}.Math.e^{0.0538.Math.\frac{d.Math..Math.\rho }{{\mu }_{\mathrm{eff}}.Math.v^{n-2}}.Math.\frac{\mathrm{.Math.}}{d}}\right)\right)}^2$$.Math.{\mu }_{\mathrm{eff}}=\frac{H}{12}.Math.{\left(9+\frac{3}{n}\right)}^n.Math.{\left(\frac{7.Math.5.Math.d^2}{1.Math.6}\right)}^{\frac{1-n}{2}}$

[0171] Where μ.sub.eff a is the effective viscosity, d is the hydraulic diameter of tubes, ρ is the density of the fluid,

[00036]$\frac{\mathrm{.Math.}}{d}$

is the relative roughness, v is the velocity of the fluid, n is the power-law index, H is the consistency coefficient, W(g) is the Lambert W function or product logarithm function.
{circle around (10)} Newton fluid+low speed non-Darcy flow or Bingham (1919) Non-Newton fluid+Darcy percolation flow

[0172] λ is written as:

[00037]$\lambda =\frac{K}{\mu }.Math.\left(1-\frac{G}{.Math.\frac{\partial p}{\partial l}.Math.}\right)$

[0173] Where μ is the viscosity of the fluid, K is the permeability, G is the low-speed non-darcy factor or pseudo-threshold pressure gradient, p is the pressure, l is the distance.

{circle around (11)} Newton fluid+slippage effect

[0174] λ is written as:

[00038]$\lambda =\frac{k\left(1+b/p\right)}{\mu }$

[0175] Where μ is the viscosity of the fluid, K is the permeability, b is the Klinkenberg factor, p is the pressure.

[0176] For non-Darcy percolation in the above formula, besides Forchheimer (1901) binomial high-speed non-Darcy formula and permeability modulus formula (stress sensitivity), there are also formulas created as a replacement by Irmay (1968), Izbash (1971), Swartzendruber (1962). Huang Tingzhang (1997), Halex (1979) and so on. The Colebrook (1939)'s along-way resistance formula above can be replaced by other formulas such as Blasius (1913). Nikuradse (1933), Churchill (1977). The Qstwald-DeWaele (1923,1925)'s power-law fluid viscosity formula above can be replaced by other formulas such as Bingham (1919), Cross (1979). Carreau (1979). Meter (1964), Sisko (1958); The Hirasaki-Pope (1974) 's shear rate formula can be replaced by other formulas such as Gogarty (1967), W. Littmann (1988), Camillen (1987), Rabinowitsch (1929), Jennings (1971). In addition, for low-density gas flow cases, the generalized mobility expression can be obtained by correcting the non-Darcy effect with the Klinkenberg equation (1941).

(1.6) Grid Design of Numerical Simulation

[0177] The design of grid types (orthogonal, corner, hybrid, etc.), grid directions (considering sediment source direction, flow direction, well location, etc.), grid boundaries, grid layers, grid sizes, grid density, etc.

(1.7) Wellbore Storage Model Setting

[0178] Wellbore storage models include constant wellbore storage model, changing wellbore storage model [Fair model (1981), Hegeman model (1993), leaky packer model and near wellbore fracture storage model (Spivey, 1999.

(1.8) Flow Period and Regime Definition

[0179] The flow period definition includes pressure build-up, pressure falloff, pressure drawdown, injection test, and slug flow. The flow regime definition includes linear flow, bilinear flow, radial flow and spherical flow.

(1.9) Coordinate Transformation (Flow Rate History Influence Treatment)

[0180] The coordinate transformation methods include equivalent constant flow rate history Horner method, equivalent constant flow rate history Agarwal method, variable flow rate history Horner method, variable flow rate history Agarwal method, deconvolution method and so on.

(1.10) Setting of Models and their Type Curve Analysis

[0181] The type curves of model analysis include various model type curves for pressure, flow and temperature analysis. The pressure analysis p type curves include full view pressure history plot, linear plot, semi-logarithmic plot, Gringarten-Bourdet log-log plot, linear flow plot, bilinear flow plot, (hemi)spherical flow plot, PPD plot, SLPD plot, early time plot, ζ function plot, curve comparison of previous single well or well group and so on.

[0182] The type curves of production analysis include blasingame plate, Agarwal Gardner plate, NPI plate and transient plate. The Type curves of temperature analysis include double logarithm diagram of temperature change and its derivative.

[0183] See Table 1 for the combination of tube flow (free flow)—percolation coupling model methods available in the numerical simulation and type curve analysis.

(1.11) Parameter Adjustment and History Matching

[0184] The adjustable parameters include permeability, saturation and porosity. History matching indexes include pressure index, production index, temperature index, fluid property index and so on.

(1.12) Dynamic Prediction

[0185] Dynamic prediction includes various development indexes of the prediction scheme for economic evaluation.

TABLE-US-00001 TABLE 1 The tube flow percolation coupling models combined in the software (Scope cases of Claims) Different combinations of models Application field Well and reservoir coupling, Combination Transient well Numerical Numerical simulation coupling between reservoirs, of fluid types analysis (pressure, simulation of oil of underground boundary coupling, etc and motion laws rate and temperature) and gas reservoir water reservoir Coupling of single/dual Remark 1 This invention Ones published Ones published by medium homogeneous by others limited others limited to use reservoir and vertical or to use of of equivalent horizontal well equivalent percolation coefficient permeability Remark 2 This invention This invention This invention Coupling of triple media Remark 1 This invention This invention Ones published by hreservoir and vertical others limited to or horizontal well use of equivalent percolation coefficient Remark 2 This invention This invention This invention Coupling of different Remark 3 This invention This invention This invention reservoirs and different well types (the above cases are not included) Coupling of different Remark 3 This invention This invention This invention reservoirs Coupling of different Remark 3 This invention This invention This invention reservoirs and different types of boundaries Remark 1 Three cases of fluid type and motion law combinations: {circle around (1)}Newtonian fluid + Darcy flow; {circle around (2)}Newtonian fluid + laminar tube flow; {circle around (4)}Newtonian fluid + turbulent tube flow; Remark 2 Seven cases of fluid type and motion law combinations: {circle around (3)}Newtonian fluid + high speed non-Darcy flow; {circle around (5)}Newtonian fluid + stress-sensitive percolation; {circle around (6)}non-Newionian fluid + Darcy flow; {circle around (7)}non-Newtonian fluid + laminar tube flow; {circle around (8)}non-Newtonian fluid + stress-sensitive percolation; {circle around (9)}non-Newtonian fluid + turbulent tube flow; {circle around (10)}Newtonian fluid + low speed non-Darcy flow (or Binghamnon-Newtonian fluid + Darcy flow). Remark 3 Ten cases of fluid type and motion law combinations: {circle around (1)}Newtonian fluid + Darcy flow; {circle around (2)}Newtonian fluid + laminar tube flow; {circle around (3)}Newtonian fluid + high speed non-Darcy flow; {circle around (4)} Newtonian fluid + turbulent tube flow; {circle around (5)} Newtonian fluid + stress-sensitive percolation; {circle around (6)}Non-Newtonian fluid + Darcy flow; {circle around (7)}Non-Newtonian fluid + laminar tube flow; {circle around (8)}Non-Newtonian fluid + stress-sensitive percolation; {circle around (9)}Non-Newtonian fluid + turbulent tube flow; {circle around (10)}Newtonian fluid + low speed non-Darcy flow (or Bingham Non-Newtonian fluid + Darcy flow). Remark 4 (1)“Published by others” are limited to the use of equivalent permeability and equivalent percolation coefficient (or converted percolation coefficient) in the above numerical simulation models, where the use of generalized mobility (or equivalent mobility) belongs to this invention; (2) This invention includes characteristic methods for tube flow (free flow) or percolation cases realized by the method of this invention; (3) The permanent downhole monitoring data analysis is mainly based on models of “transient well analysis” (pressure or temperature) or/and “transient rate production data analysis”. (4) This invention includes above models that can be applied for geothermal reservoirs. Remark 5 The reservoirs refer to a reservoir of oil, gas or water. Reservoirs contains tubes, porous media, fracture media, caves, or their combinations, connected combinations or overlapping.

[0186] (2) The software developed can be used for single and multi-well flow simulation of complex multi-component multi-component multi-phase flow reservoirs, multi-well interference analysis, deliverability analysis (i.e. steady well test analysis), transient pressure analysis (i.e. unsteady pressure well test analysis), transient rate analysis (i.e. production data analysis or rate decline analysis), transient temperature analysis (i.e. unsteady temperature analysis), well test design, and permanent downhole monitoring data analysis (including pressure, temperature, flow rate, and so on). The above complex multi-component multi-phase flow reservoirs refer to all types of reservoirs, reservoirs with fluid injection, underground gas storage, ground water reservoirs, geothermal reservoirs, etc.

[0187] Embodiment 1 is a general example that can be used to solve tube flow-percolation coupling problems and multi-phase flow problems for various complex reservoirs, embodiments 2 to 7 are simplified on the basis of embodiment 1 to reflect that the wellbore can be treated either as an inner boundary or as a reservoir, the specific solution methods of which can be solved according to the numerical methods given in embodiment 1. However, some embodiments are simple and can be solved by analytical methods. For example, embodiments 4 and 6 can be solved by Laplace integral transformation method and give analytical solutions in Laplace space. embodiments 8 to 12 are ones given where the generalized mobility is a constant and they have analytical solutions in Laplace space. embodiment 13 gives a situation where the generalized mobility is a non-constant. embodiment 14 is an example of multi-component condensate gas flow simulation analysis for complex reservoirs. Embodiment 15 is an example of multi-component carbon dioxide drive flow simulation analysis in complex reservoirs. Embodiment 16 is an example of multi-component polymer drive flow simulation analysis in complex reservoirs. Embodiment 17 is an example of multi-component three-phase flow simulation analysis considering non-isothermal process in complex reservoirs. Embodiment 18 is an example of transient temperature analysis of oil single-phase flow.

[0188] It should be noted that the generalized mobility is a constant (matrix) when only the fluid flow is Darcy flow or laminar tube flow. In general, the generalized mobility is a function of space position and time. In particular, the generalized mobility may be a function of variables such as fluid velocity, pressure gradient, and so on. If the generalized mobility can be expressed as the product of a constant matrix and a scalar function, the non-linear flow problem can be transformed into linear flow through defining pseudo-pressure and pseudo-time functions, which means to transform generalized mobility into a constant matrix; Otherwise, the problem cannot be transformed into a linear flow problem.

[0189] Let a function ƒ(t) be defined on [0,∞), then ƒ(t) is a real-valued function or a complex-valued function of a real variable t. The function

[00039]$\stackrel{\_}{f}\left(z\right)=\underset{0}{\overset{\infty }{\int }}.Math.f\left(t\right).Math.e^{-\mathrm{zt}}.Math.\mathrm{dt}$

determined by the Laplace integral

[00040]$\underset{0}{\overset{\infty }{\int }}.Math.f\left(t\right).Math.e^{-\mathrm{zt}}.Math.\mathrm{dt}$

is called the Laplace transformation of the function ƒ(t). By applying the integral transformation above, the well testing models in embodiments 4, 6 and 8 to 12 can be transformed into homogeneous equation systems in Laplace space, which are used to obtain bottom hole pressure function solutions in Laplace space.

[0190] Let p.sub.wD(z) be dimensionless bottom hole pressure solution function in Laplace space, then the real space dimensionless bottom hole pressure solution p.sub.wD(t.sub.D) can be obtained by Stehfest numerical inversion technique

[00041]$p_{\mathrm{wD}}\left(t_D\right)=\frac{\log .Math..Math.2}{t_D}.Math.\underset{k=1}{\overset{N}{.Math.}}.Math.\left({\left(-1\right)}^{\frac{N}{2}+k}.Math.\underset{n=\left[\left(k+1\right)/2\right]}{\overset{\min \left(k,N/2\right)}{.Math.}}.Math.\frac{n^{N/2}\left(2.Math.n\right)!}{\left(N/2-n\right)!.Math.n!.Math.\left(n-1\right)!.Math.\left(k-n\right)!.Math.\left(2.Math.n-k\right)!}\right).Math.{\stackrel{\_}{p}}_{\mathrm{wD}}\left(k.Math.\frac{\log .Math..Math.2}{t_D}\right)$

where N is an even number, generally the value is between 8 and 16.

[0191] Other numerical inversion techniques can also be used here in addition to Stehfest numerical inversion technology.

[0192] Embodiment of transient rate analysis (i.e., production data analysis) is shown in FIGS. 41 to 45. These figures are realized on the basis of a cylindrical reservoir model provided in embodiment 9 of this invention, which are vertical well Blasingame type curve, vertical well Agarwal-Gardner type curve, vertical well NPI type curve, and vertical well Transient type curve. Similarly, various kinds of analysis type curves of other analysis models can also be created. The basic principles of drawing type curves are detailed in the book “Modern production decline analysis methods for oil and gas wells and their application” written by Sun Hedong (2013).

[0193] For some layered reservoirs, the problem can be further simplified by defining the generalized flow coefficient or the equivalent transmissibility (Kh/μ) to facilitate the analysis and solution of some problems.

[0194] In addition, based on the above principle method, It can be similar to follow the traditional analytical model methods to establish the corresponding deliverability well test analysis, multi-well interference well test analysis, well test analysis method for wells in geothermal reservoirs and permanent downhole monitoring data analysis.

[0195] For example, with regard to the establishment of single-phase flow deliverability test analysis model method for steady-flow vertical wells in unsaturated layered reservoirs, the formula of production index is written as: follows:

[00042]$J=\frac{q}{P_e-P_{\mathrm{wf}}}-\frac{2.Math.\mathrm{\pi \lambda h}}{B\left[\ln \left(r_e/r_w\right)+S\right]}$

[0196] B—Volume factor of the fluid, m.sup.3/m.sup.3; h—Reservoir thickness, m; p.sub.e-reservoir outer boundary pressure (average formation pressure of reservoir can be used in practical application), Pa; pwf-bottom hole flowing pressure, Pa; q—flow rate, m.sup.3/s; r.sub.e-outer radius of reservoir, m; r.sub.w-wellbor radius, m S-skin factor; λ-generalized mobility of the reservoir fluid, m.sup.2/Pa.Math.s.

Embodiment 1

[0197] This embodiment provides a generalized oil/gas/water three-phase flow simulation analysis method for complex reservoirs. FIG. 3 is the physical model, which only denotes a special case of this embodiment.

[0198] Based on the mass conservation principle, the general equation of oil/gas/water three-phase flow can be established (Note: the volume factor is equal to the fluid density divided by the surface reference density and the surface reference density is a constant. Therefore, the two sides of the governing equation in the following examples can be converted into the fluid density by multiplying the surface reference density of each phase, respectively):

[0199] Oil phase governing equation considering source/sink term:

[00043]$\begin{array}{cc}\nabla .Math.\left[\frac{1}{B_o}.Math.{\lambda }_o.Math.\nabla p_o\right]=\frac{\partial }{\partial t}.Math.\left(\frac{1}{B_o}.Math.{\phi .Math.S}_o\right)+\frac{q_o}{{\rho }_{\mathrm{osc}}},\left(x,t\right)\in \Omega \times \left(0,t_{\max }\right]& \left(1\right)\end{array}$

[0200] Water phase governing equation considering source/sink term:

[00044]$\begin{array}{cc}\nabla .Math.\left[\frac{1}{B_w}.Math.{\lambda }_w.Math.\nabla p_w\right]=\frac{\partial }{\partial t}.Math.\left(\frac{1}{B_w}.Math.{\phi .Math.S}_w\right)+\frac{q_w}{{\rho }_{\mathrm{osc}}},\left(x,t\right)\in \Omega \times \left(0,t_{\max }\right]& \left(2\right)\end{array}$

[0201] Gas phase governing equation considering source/sink term:

[00045]$\begin{array}{cc}\nabla .Math.\left[\frac{1}{B_g}.Math.{\lambda }_g.Math.{\nabla }_g\right]+\nabla .Math.\left[\frac{R_g}{B_o}.Math.{\lambda }_o.Math.\nabla p_o\right]=\frac{\partial }{\partial t}.Math.\left(\frac{1}{B_g}.Math.\phi ,S_g\right)+\frac{q_g}{{\rho }_{\mathrm{gsc}}}+\frac{\partial }{\partial t}.Math.\left(\frac{R_s}{B_o}.Math.{\phi .Math.S}_o\right),\left(x,t\right)\in \Omega \times \left(0,t_{\max }\right]& \left(3\right)\end{array}$

[0202] Auxiliary equations with saturation and capillary pressure:

S.sub.o+S.sub.w+S.sub.g=1,(x,t)∈Ω×(0,t.sub.max]  (4)

p.sub.cow(S.sub.w)=p.sub.o−p.sub.w,(x,t)∈Ω×(0,t.sub.max]  (5)

p.sub.cgo(S.sub.w)=p.sub.g−p.sub.o,(x,t)∈Ω×(0,t.sub.max]  (6)

[0203] Boundary condition equations:

[00046]$\begin{array}{cc}\left(c_{o,1}.Math.p_o+c_{o,2}.Math.{\lambda }_o.Math.\frac{\partial p_o}{\partial n^{\partial \Omega }}\right)=g_o\left(x,t\right),\left(x,t\right)\in \partial \Omega \times \left(0,t_{\max }\right]& \left(7\right)\\ \left(c_{w,1}.Math.p_w+c_{w,2}.Math.{\lambda }_w.Math.\frac{\partial p_w}{\partial n^{\partial \Omega }}\right)=g_w\left(x,t\right),\left(x,t\right)\in \partial \Omega \times \left(0,t_{\max }\right]& \left(8\right)\\ \left(c_{g,1}.Math.p_g+c_{g,2}.Math.{\lambda }_g.Math.\frac{\partial p_g}{\partial n^{\partial \Omega }}\right)=g_g\left(x,t\right),\left(x,t\right)\in \partial \Omega \times \left(0,t_{\max }\right]& \left(9\right)\end{array}$

[0204] Initial condition equations:

S.sub.k(x,0)=l.sub.k(x),x∈Ω  (10)

p.sub.k(x,0)=ƒ.sub.k(x),x∈Ω  (11)

Where

S.sub.k=S.sub.k(x,t),k=o,w,g  (12)

p.sub.k=p.sub.k(x,t),k=o,w,g  (13)

q.sub.k=q.sub.k(x,t),k=o,w,g  (14)

λ.sub.k=λ.sub.k(x,t),k=o,w,g  (15)

[0205] ϕ is porosity, which is the function of average pressure p, %; B.sub.w, B.sub.o, B.sub.g are volume factors of water/oil/gas, water volume factor is a function of water-phase pressure p.sub.w, oil volume factor is a function of oil-phase pressure p.sub.o and bubble point pressure p.sub.b, gas volume factor is a function of gas-phase pressure p.sub.g, dimensionless; λ.sub.w, λ.sub.o, λ.sub.g are generalized mobility for water/oil/gas, separately, m.sup.2/(Pa.Math.s); ρ.sub.w, ρ.sub.o, ρ.sub.g are densities of water/oil/gas, kg/m.sup.3; ρ.sub.o,o is oil component density in oil-phase, kg/m.sup.3; ρ.sub.g,o is gas component density in oil-phase, kg/m.sup.3; ρ.sub.wsc, ρ.sub.osc, ρ.sub.gsc are surface reference densities of water/oil/gas, kg/m.sup.3; R.sub.s is the dissolved gas-oil ratio, dimensionless; S.sub.w, S.sub.o, S.sub.g are saturation of water/oil/gas; t is time, s; t.sub.max is the total time, s; p.sub.w, p.sub.o, p.sub.g are pressures of water/oil/gas, Pa; l.sub.w, l.sub.o, l.sub.g are distribution functions of water/oil/gas initial saturation; q.sub.w, q.sub.o, q.sub.g are source/sink terms of water/oil/gas, kg/(m.sup.3.Math.s); ƒ.sub.w, ƒ.sub.o, ƒ.sub.g are distribution functions of water/oil/gas initial pressure, Pa; g.sub.w, g.sub.o, g.sub.g are boundary functions of water/oil/gas on reservoir boundary, dimensionless; p.sub.cow, p.sub.cgo are capillary pressure of oil-water/gas-oil contact, Pa; Ω is reservoir space; ∂Ω is reservoir boundary including internal boundary and outer boundary; c.sub.w1, c.sub.o,1, c.sub.g,1 are pressure term coefficients of water/oil/gas on reservoir boundary, 1/Pa; c.sub.w,2, c.sub.o,2, c.sub.g,2 are respectively directional derivative coefficient of water, oil and gas along the normal direction outside the boundary on reservoir boundary condition, s/m; n.sup.∂Ω is The direction of outer normal of reservoir boundary, m; ∇ is the Hamilton operator; ∂ is the partial derivative sign.

[0206] Problems can be optimized and discretely solved by applying the finite volume method. The following example of tetrahedral element is given to illustrate the process.

[0207] Firstly the region Ω is divided into a series of tetrahedrons Ω.sub.i, which do not overlap each other. The control point of tetrahedrons Ω.sub.i is marked as x.sub.i and the boundary of tetrahedrons Ω.sub.i is marked as ∂Ω.sub.i. If the number of tetrahedrons that have a common surface with tetrahedrons Ω.sub.i is 4, then x.sub.i is called the internal control point, and the set of all internal control points is Ω.sub.I. If the number of tetrahedrons that have a common surface with tetrahedrons is less than Ω.sub.i then x.sub.i is called the boundary control point, and the set of all boundary control points is called Ω.sub.B.

[0208] Oil phase governing equation considering source/sink term is discretized as

[00047]$\begin{array}{cc}\underset{\sigma .Math.\partial {\Omega }_{\eta }}{.Math.}.Math.u_{o,i_2}.Math.a_{o,i_2}.Math.p_o\left(x_{i_2},t_{n+1}\right)-u_{o,i_1}.Math.a_{o,i_1}.Math.p_o\left(x_{i_1},i_{n+1}\right)=e_{o,2}.Math.S_o\left(x_{i_1},i_{n+1}\right)-e_{o,1}.Math.S_o\left(x_{i_1},t_n\right)+V_{i_1}.Math.\frac{q_o\left(x_{i_1},t_{n+1}\right)}{{\rho }_{\mathrm{osc}}}& \left(16\right)\end{array}$

[0209] Where x.sub.i.sub.1∈Ω.sub.I, x.sub.i.sub.2∈Ω.sub.I UΩ.sub.B, and

[00048]$\begin{array}{cc}.Math.e_{o,1}=\frac{V_{i_1}.Math.\phi \left[\stackrel{\_}{p}\left(x_{i_1},t_n\right)\right]}{B_o\left[p_o\left(x_{i_1},t_n\right),p_b\right].Math.\left(t_{n+1}-t_n\right)}& \left(17\right)\\ .Math.e_{o,2}=\frac{V_{i_1}.Math.\phi \left[\stackrel{\_}{p}\left(x_{i_1},t_{n+1}\right)\right]}{B_o\left[p_o\left(x_{i_1},t_{n+1}\right),p_b\right].Math.\left(t_{n+1}-t_n\right)}& \left(18\right)\\ .Math.a_{o,i_1}=S_{i_1,\sigma }\left({\alpha }_{o,i_1}+{\beta }_{o,i_1}+{\gamma }_{o,i_1}\right)& \left(19\right)\\ .Math.{\alpha }_{o,i_2}=S_{i_2,\sigma }\left({\alpha }_{o,i_2}+{\beta }_{o,i_2}+{\gamma }_{o,i_2}\right)& \left(20\right)\\ b_{o,i_1}=S_{i_1,\sigma }\left[{\alpha }_{o,i_1}.Math.p_o\left(d_{j_1},t_{n+1}\right)+{\beta }_{o,i_1}.Math.p_o\left(d_{j_2},t_{n+1}\right)+{\gamma }_{o,i_1}.Math.p_o\left(d_{j_4},t_{n+1}\right)\right]& \left(21\right)\\ b_{o,i_2}=S_{i_2,\sigma }\left[{\alpha }_{o,i_2}.Math.p_o\left(d_{j_2},t_{n+1}\right)+{\beta }_{o,i_2}.Math.p_o\left(d_{j_3},t_{n+1}\right)+{\gamma }_{o,i_2}.Math.p_o\left(d_{j_5},t_{n+1}\right)\right]& \left(22\right)\\ .Math.u_{o,i_1}=\{\begin{array}{cc}1/2,& b_{o,i_1}+b_{o,i_2}=0\\ b_{o,i_2}/\left(b_{o,i_1}+b_{o,i_2}\right),& b_{o,i_1}+b_{o,i_2}\ne 0\end{array}& \left(23\right)\\ .Math.u_{o,i_2}=\{\begin{array}{cc}1/2,& b_{o,i_1}+b_{o,i_2}=0\\ b_{o,i_1}/\left(b_{o,i_1}+b_{o,i_2}\right),& b_{o,i_1}+b_{o,i_2}\ne 0\end{array}& \left(24\right)\end{array}$

[0210] The vertex pressures in the formulas [(21) and (22)] above are treated by the inverse distance weighting method, i. e,

[00049]$\begin{array}{cc}{p}_o\left(d_j,t_{n+1}\right)=\underset{i\in D_j}{.Math.}.Math.\frac{||d_j-x_i.Math.{||}_2^l}{\underset{i\in D_j}{.Math.}.Math.||d_j-x_i.Math.{||}_2^l}.Math.p_o\left(x_i,t_{n+1}\right),j=j_1,j_2,j_3,j_4,j_5& \left(25\right)\end{array}$

[0211] D.sub.j denotes the number set of all discrete control body central points from the vertex d.sub.j, l is a real number less than 0.

[0212] The coefficients α.sub.o,i.sub.1, β.sub.o,i.sub.1, γ.sub.o,i.sub.1 in the formulas [(19) and (21) above are determined by the following formula

[00050]$\begin{array}{cc}\frac{{{\lambda }_o}^T\left(x_{i_1},t_{n+1}\right).Math.n_{i_1,\sigma }}{B_o\left[p_o\left(x_{i_1},t_{n+1}\right),p_b\right]}={\alpha }_{o,i_1}\left(d_{j_1}-x_{i_1}\right)+{\beta }_{o,i_1}\left(d_{j_2}-x_{i_1}\right)+{\gamma }_{o,i_1}\left(d_{j_4}-x_{i_1}\right),x_{i_1}\in {\Omega }_i& \left(26\right)\end{array}$

[0213] As shown in FIG. 4, the subscripts of vertices coordinates j.sub.1, j.sub.2, j.sub.4, are intersected by a ray in the direction of

[00051]$\frac{{{\lambda }_o}^T\left(x_{i_1},t_{n+1}\right).Math.n_{i_1,\sigma }}{B_o\left[p_o\left(x_{i_1},t_{n+1}\right),p_b\right]}$

through the starting point x.sub.i.sub.1 with a sub-surface of the discrete control body Ω.sub.i.sub.1. The three vertices coordinates of the sub-surface are d.sub.j.sub.1, d.sub.j.sub.2, d.sub.j.sub.4, as shown in FIG. 4.

[0214] The coefficients α.sub.o,i.sub.2, β.sub.o,i.sub.2, γ.sub.o,i.sub.2 in the formulas Eq,20 and Eq.22 above are determined by the following formula

[00052]$\begin{array}{cc}\frac{{{\lambda }_o}^T\left(x_{i_2},t_{n+1}\right).Math.n_{i_2,\sigma }}{B_o\left[p_o\left(x_{i_2},t_{n+1}\right),p_b\right]}={\alpha }_{o,i_2}\left(d_{j_2}-x_{i_2}\right)+{\beta }_{o,i_2}\left(d_{j_3}-x_{i_2}\right)+{\gamma }_{o,i_2}\left(d_{j_5}-x_{i_2}\right)& \left(27\right)\end{array}$

[0215] The subscripts of vertices coordinates j.sub.2, j.sub.3, j.sub.5 are intersected by a ray in the direction of

[00053]$\frac{{{\lambda }_o}^T\left(x_{i_2},t_{n+1}\right).Math.n_{i_2,\sigma }}{B_o\left[p_o\left(x_{i_2},t_{n+1}\right),p_b\right]}$

through the starting point x.sub.i.sub.2 with a sub-surface of the discrete control body χ.sub.i.sub.2. The three vertices coordinates of the sub-surface are d.sub.j.sub.2, d.sub.j.sub.3, d.sub.j.sub.5, as shown in FIG. 4.

[0216] Water phase governing equation considering source/sink term is discretized as

[00054]$\begin{array}{cc}\underset{\sigma .Math.\partial {\Omega }_{\eta }}{.Math.}.Math.u_{w,i_2}.Math.a_{w,i_2}.Math.p_w\left(x_{i_2},t_{n+1}\right)-u_{w,i_1}.Math.a_{w,i_1}.Math.p_w\left(x_{i_1},t_{n+1}\right)=e_{w,2}.Math.S_w\left(x_{i_1},t_{n+1}\right)-e_{w,1}.Math.S_w\left(x_{i_1},t_n\right)+V_{i_1}.Math.\frac{q_w\left(x_{i_1},t_{n+1}\right)}{{\rho }_{\mathrm{wsc}}}& \left(28\right)\end{array}$

[0217] Where x.sub.i.sub.1∈Ω.sub.I,

[00055]$\begin{array}{cc}.Math.e_{w,1}=\frac{V_{i_1}.Math.\phi \left[\stackrel{\_}{p}\left(x_{i_1},t_n\right)\right]}{B_w\left[p_w\left(x_{i_1},t_n\right)\right].Math.\left(t_{n+1}-t_n\right)}& \left(29\right)\\ .Math.e_{w,2}=\frac{V_{i_1}.Math.\phi \left[\stackrel{\_}{p}\left(x_{i_1},t_{n+1}\right)\right]}{B_w\left[p_w\left(x_{i_1},t_{n+1}\right)\right].Math.\left(t_{n+1}-t_n\right)}& \left(30\right)\\ .Math.a_{w,i_1}=S_{i_1,\sigma }\left({\alpha }_{w,i_1}+{\beta }_{w,i_1}+{\gamma }_{w,i_1}\right)& \left(31\right)\\ .Math.{\alpha }_{w,i_2}=S_{i_2,\sigma }\left({\alpha }_{w,i_2}+{\beta }_{w,i_2}+{\gamma }_{w,i_2}\right)& \left(32\right)\\ b_{w,i_1}=S_{i_1,\sigma }\left[{\alpha }_{w,i_1}.Math.p_w\left(d_{j_1},t_{n+1}\right)+{\beta }_{w,i_1}.Math.p_w\left(d_{j_2},t_{n+1}\right)+{\gamma }_{w,i_1}.Math.p_w\left(d_{j_4},t_{n+1}\right)\right]& \left(33\right)\\ b_{w,i_2}=S_{i_2,\sigma }\left[{\alpha }_{w,i_2}.Math.p_w\left(d_{j_2},t_{n+1}\right)+{\beta }_{w,i_2}.Math.p_w\left(d_{j_3},t_{n+1}\right)+{\gamma }_{w,i_2}.Math.p_w\left(d_{j_5},t_{n+1}\right)\right]& \left(34\right)\\ .Math.u_{w,i_1}=\{\begin{array}{cc}1/2,& b_{w,i_1}+b_{w,i_2}=0\\ b_{w,i_2}/\left(b_{w,i_1}+b_{w,i_2}\right),& b_{w,i_1}+b_{w,i_2}\ne 0\end{array}& \left(35\right)\\ .Math.u_{w,i_2}=\{\begin{array}{cc}1/2,& b_{w,i_1}+b_{w,i_2}=0\\ b_{o,i_1}/\left(b_{w,i_1}+b_{w,i_2}\right),& b_{w,i_1}+b_{w,i_2}\ne 0\end{array}& \left(36\right)\end{array}$

[0218] The vertex pressures in the formulas Eq.37 and Eq.34 above are treated by the inverse distance weighting method, i. e,

[00056]$\begin{array}{cc}{p}_w\left(d_j,t_{n+1}\right)=\underset{i\in D_j}{.Math.}.Math.\frac{||d_j-x_i.Math.{||}_2^l}{\underset{i\in D_j}{.Math.}.Math.||d_j-x_i.Math.{||}_2^l}.Math.p_w\left(x_i,t_{n+1}\right),j=j_1,j_2,j_3,j_4,j_5& \left(37\right)\end{array}$

[0219] D.sub.j denotes the number set of all discrete control body central points from the vertex d.sub.j, l is a real number less than 0.

[0220] The coefficients α.sub.w,i.sub.1, β.sub.w,i.sub.1, γ.sub.w,i.sub.1 in the formulas [(31) and (33)] above are determined by the following formula

[00057]$\begin{array}{cc}\frac{{\lambda }_w^T\left({\varkappa }_{i_1},t_{n+1}\right).Math.n_{i_1,\sigma }}{B_w\left[p_w\left({\varkappa }_{i_1},t_{n+1}\right)\right]}={\alpha }_{w,i_1}\left(d_{j_1}-{\varkappa }_{i_1}\right)+{\beta }_{w,j_1}\left(d_{j_2}-{\varkappa }_{i_1}\right)+{\gamma }_{w,j_1}\left(d_{j_4}-{\varkappa }_{i_1}\right),{\varkappa }_{i_1}\in {\Omega }_I& \left(38\right)\end{array}$

[0221] The subscripts of vertices coordinates j.sub.1, j.sub.2, j.sub.4 are intersected by a ray in the direction of

[00058]$\frac{{\lambda }_w^T\left({\varkappa }_{i_1},t_{n+1}\right).Math.n_{i_1,\sigma }}{B_w\left[p_w\left({\varkappa }_{i_1},t_{n+1}\right)\right]}$

through the starting point x.sub.i.sub.1 with a sub-surface of the discrete control body Ω.sub.i.sub.1. The three vertices coordinates of the sub-surface are d.sub.j.sub.1, d.sub.j.sub.2, d.sub.j.sub.4 as shown in FIG. 5.

[0222] The coefficients α.sub.w,i.sub.2, β.sub.w,i.sub.2, γ.sub.w,i.sub.2, in the formulas [(32) and (34)] above are determined by the following formula

[00059]$\begin{array}{cc}\frac{{\lambda }_w^T\left({\varkappa }_{i_2},t_{n+1}\right).Math.n_{i_2,\sigma }}{B_w\left[p_w\left({\varkappa }_{i_2},t_{n+1}\right)\right]}={\alpha }_{w,i_2}\left(d_{j_2}-{\varkappa }_{i_2}\right)+{\beta }_{w,j_2}\left(d_{j_3}-{\varkappa }_{i_2}\right)+{\gamma }_{w,i_2}\left(d_{j_5}-{\varkappa }_{i_2}\right)& \left(39\right)\end{array}$

[0223] The subscripts of vertices coordinates j.sub.2, j.sub.3, j.sub.5 are intersected by a ray in the direction of

[00060]$\frac{{\lambda }_w^T\left({\varkappa }_{i_2},t_{n+1}\right).Math.n_{i_2,\sigma }}{B_w\left[p_w\left({\varkappa }_{i_2},t_{n+1}\right)\right]}$

through the starting point x.sub.i.sub.2 with a sub-surface of the discrete control body Ω.sub.i.sub.2. The three vertices coordinates of the sub-surface are d.sub.j.sub.2, d.sub.j.sub.3, d.sub.j.sub.5 as shown in FIG. 5.

[0224] Gas phase governing equation considering source/sink term (Eq.3) is discretized as

[00061]$\begin{array}{cc}\underset{\sigma .Math.\partial {\Omega }_{i_1}}{.Math.}.Math..Math.\left\{\begin{array}{c}{u}_{g,i_2}\left[a_{g,i_2}.Math.p_g\left({\varkappa }_{i_2},t_{n+1}\right)+a_{o,i_2}^{\prime }.Math.p_o\left({\varkappa }_{i_2},t_{n+1}\right)\right]-\\ u_{g,i_1}\left[a_{g,i_1}.Math.p_g\left({\varkappa }_{i_1},t_{n+1}\right)+a_{o,i_1}^{\prime }.Math.p_o\left({\varkappa }_{i_1},t_{n+1}\right)\right]\end{array}\right\}=\hspace{1em}[.Math.\begin{array}{c}{e}_{g,2}.Math.S_g\left({\varkappa }_{t_1,t_{n+1}}\right)-e_{g,1}.Math.S_g\left({\varkappa }_{t_1},t_n\right)+\\ e_{0,2}^{\prime }.Math.S_0\left({\varkappa }_{t_1},t_{n+1}\right)-e_{0,1}^{\prime }.Math.S_0\left({\varkappa }_{t_1},t_n\right)\end{array}]+V_{i_1}.Math.\frac{q_g\left({\varkappa }_{i_1},t_{n+1}\right)}{{\rho }_{\mathrm{gsc}}}& \left(40\right)\end{array}$

[0225] Where x.sub.i.sub.1∈Ω.sub.I,

[00062]$\begin{array}{cc}{e}_{g,1}=\frac{V_{i_1}.Math.\phi \left[\stackrel{\_}{p}\left({\varkappa }_{i_1},t_n\right)\right]}{B_g\left[p_g\left({\varkappa }_{i_1},t_n\right)\right].Math.\left(t_{n+1}-t_n\right)}& \left(41\right)\\ e_{g,2}=\frac{V_{i_1}.Math.\phi \left[\stackrel{\_}{p}\left({\varkappa }_{i_1},t_{n+1}\right)\right]}{B_g\left[p_g\left({\varkappa }_{i_1},t_{n+1}\right)\right].Math.\left(t_{n+1}-t_n\right)}& \left(42\right)\\ e_{o,1}^{\prime }=\frac{V_{i_1}.Math.\phi \left[\stackrel{\_}{p}\left({\varkappa }_{i_1},t_n\right)\right]}{B_o\left[p_o\left({\varkappa }_{i_1},t_n\right),p_b\right].Math.\left(t_{n+1}-t_n\right)}& \left(43\right)\\ e_{o,2}^{\prime }=\frac{V_{i_1}.Math.\phi \left[\stackrel{\_}{p}\left({\varkappa }_{i_1},t_n\right)\right]}{B_o\left[p_o\left({\varkappa }_{i_1},t_{n+1}\right),p_b\right].Math.\left(t_{n+1}-t_n\right)}& \left(44\right)\\ a_{g,i_1}=S_{i_1,\sigma }\left({\alpha }_{g,i_1}+{\beta }_{g,i_1}+{\gamma }_{g,i_1}\right)& \left(45\right)\\ a_{g,i_2}=S_{i_2,\sigma }\left({\alpha }_{g,i_2}+{\beta }_{g,i_2}+{\gamma }_{g,i_2}\right)& \left(46\right)\\ a_{o,i_1}^{\prime }=S_{i_1,\sigma }\left({\alpha }_{o,i_1}^{\prime }+{\beta }_{o,i_1}^{\prime }+{\gamma }_{o,i_1}^{\prime }\right)& \left(47\right)\\ a_{o,i_2}^{\prime }=S_{i_2,\sigma }\left({\alpha }_{o,i_2}^{\prime }+{\beta }_{o,i_2}^{\prime }+{\gamma }_{o,i_2}^{\prime }\right)& \left(48\right)\\ b_{g,i_1}=S_{i_1,\sigma }\left[{\alpha }_{g,i_1}.Math.p_g\left(d_{j_1},t_{n+1}\right)+{\beta }_{g,i_1}.Math.p_g\left(d_{j_3},t_{n+1}\right)+{\gamma }_{g,i_1}.Math.p_g\left(d_{j_1},t_{n+1}\right)\right]& \left(49\right)\\ b_{g,i_2}=S_{i_2,\sigma }\left[{\alpha }_{g,i_2}.Math.p_g\left(d_{j_2},t_{n+1}\right)+{\beta }_{g,i_2}.Math.p_g\left(d_{j_3},t_{n+1}\right)+{\gamma }_{g,i_2}.Math.p_g\left(d_{j_5},t_{n+1}\right)\right]& \left(50\right)\\ b_{o,i_1}^{\prime }=S_{i_1,\sigma }\left[{\alpha }_{o,i_1}^{\prime }.Math.p_o\left(d_{j_1},t_{n+1}\right)+{\beta }_{o,i_1}^{\prime }.Math.p_o\left(d_{j_2},t_{n+1}\right)+{\gamma }_{o,i_1}^{\prime }.Math.p_o\left(d_{j_4},t_{n+1}\right)\right]& \left(51\right)\\ b_{o,i_2}^{\prime }=S_{i_2,\sigma }\left[{\alpha }_{o,i_2}^{\prime }.Math.p_o\left(d_{j_2},t_{n+1}\right)+{\beta }_{o,i_2}^{\prime }.Math.p_o\left(d_{j_3},t_{n+1}\right)+{\gamma }_{o,i_2}^{\prime }.Math.p_o\left(d_{j_5},t_{n+1}\right)\right]& \left(52\right)\\ u_{g,i_1}=\{\begin{array}{c}1/2,\left(b_{g,i_1}+b_{o,i_1}^{\prime }\right)+\left(b_{g,i_2}+b_{o,i_2}^{\prime }\right)=0\\ \left(b_{g,i_2}+b_{o,i_2}^{\prime }\right)/\left(b_{g,i_1}+b_{o,i_1}^{\prime }+b_{g,i_2}+b_{o,i_2}^{\prime }\right),\left(b_{g,i_1}+b_{o,i_1}^{\prime }\right)+\left(b_{g,i_2}+b_{o,i_2}^{\prime }\right)\ne 0\end{array}& \left(53\right)\\ u_{g,i_2}=\{\begin{array}{c}1/2,\left(b_{g,i_1}+b_{o,i_1}^{\prime }\right)+\left(b_{g,i_2}+b_{o,i_2}^{\prime }\right)=0\\ \left(b_{g,i_1}+b_{o,i_1}^{\prime }\right)/\left(b_{g,i_1}+b_{o,i_1}^{\prime }+b_{g,i_2}+b_{o,i_2}^{\prime }\right),\left(b_{g,i_1}+b_{o,i_1}^{\prime }\right)+\left(b_{g,i_2}+b_{o,i_2}^{\prime }\right)\ne 0\end{array}& \left(54\right)\end{array}$