The present invention concerns a method for determining at least one characteristic value of a multiphase flow at each instant.

The wording "multiphase flow" refers to any fluid flow consisting of more than one phase. The phases are distinguished by different characteristics on the molecular level. Hence, it is common to distinguish gas/solids flows, liquid/solids flows, gas/particle flows or gas/liquid flows.

An example of a gas/liquid flow is a hydrocarbon fluids flow in a pipeline system.

Such flow comprises notably gas bubbles distributed over a liquid phase or on the contrary, liquid bubbles distributed over a gas phase. The concentration of the bubbles depends mainly on flow conditions which can change for example in function of exterior environmental conditions.

The study of a multiphase pipeline flow presents an important issue for the oil and gas industry. Indeed, a pressure drop in the pipeline system due to a pipeline wall friction or hold-up regions formed by the bubbles may be determined using characteristic values of the multiphase flow.

Two types of characteristic values are usually defined in the multiphase flow.

A first type of characteristic values, called local characteristic values, allows the characterization of each phase of the multiphase flow independently of the others. An example of such local characteristic values is a local density or a local viscosity of the corresponding phase.

Thus, for example, the local density of a liquid phase is higher than the local density of a gas phase and can be determined independently of this last.

A second type of characteristic values, called global characteristic values, allows the characterization of the multiphase flow independently of its phases. An example of such global characteristic values is a multiphase flow velocity vector or a flow pressure. In particular, the velocity vector may comprise two or three components depending on the pipeline system configuration.

In the art, different methods allowing the determination of at least some of characteristic values of a multiphase flow have been proposed.

The document WO 201 1 159372 A1 discloses a method for modeling properties of a reservoir containing hydrocarbon fluids at least at two different phases. In order to determine characteristic values of such flow (density or viscosity for example), the method uses the Darcy's Law for modelling a multiphase flow. The document US 6,038,389 discloses a method of modeling a hydrocarbon flow in a heterogeneous underground reservoir containing hydrocarbon fluids. The modelling is carried out using the Darcy's Law and Galerkin numerical solutions of this law.

The document WO 02/03103 A2 discloses a numerical method for modeling a multiphase flow using finite element models.

However, the proposed methods of determination of the characteristic values may lead to inaccurate results, their implementation is relatively complex and needs a relatively big quantity of computational resources.

One aim of the invention is to provide a method for determining at least one characteristic value of a multiphase flow which is particularly accurate, allows the determination of this value in a simple and rapid way even with a limited quantity of computational resources.

To this end, the subject-matter of the invention is a method for determining at least one characteristic value of a multiphase flow at each instant, the multiphase flow comprising a flow of hydrocarbon fluids in a flow area, each phase of the flow being characterized by a local density and a local viscosity of this phase;

the method comprising the following steps:

- providing reference parameters corresponding to at least one characteristic value of the multiphase flow defined at each point of the flow area for different sampling instants;

- for at least said characteristic value, determining a temporal correlation matrix for all of the corresponding references parameters;

- for the or each determined correlation matrix:

+ determining the eigenvalues and eigenvectors of the determined correlation matrix;

+ determining space normal modes corresponding to a linear combination of eigenvectors of said correlation matrix and the corresponding references parameters;

+ determining a number of significant normal modes;

- determining coefficients of representative tensors associated to a time differential equation system, said equation system comprising time variables and resulting from the Galerkin projection of time discretized Navier-Stokes type equations on the significant normal modes;

the coefficients of each representative tensor being determined using the significant normal modes of at least one correlation matrix, the coefficients of at least one representative tensor being further determined using reduced constants, each reduced constant being function of a combination of local densities and/or local viscosities of the phases of the multiphase flow;

- determining time variables of said equation system at a next instant using time variables at a current instant and the coefficients of the representatives tensors;

- determining at least one characteristic value at each instant using the temporal variables at this instant and the corresponding significant normal modes.

The method according to the invention may comprise one or more of the following feature(s), taken in isolation, or according to any one of any technically feasible combination:

- the time variables of the time differential equation system correspond to the eigenvectors of least one correlation matrix;

- each characteristic value is a value of a space and time function determined at a given point of the flow area for a given instant, and is chosen in the group comprising:

- flow pressure at the given point of the flow area for the given instant; and - flow velocity according to one or several predetermined directions at the given point of the flow area for the given instant;

- at least one characteristic value comprises a global density, the global density being a function of space and time equal at a given point of the flow area for a given instant to the local density of the flow phase at the given point for the given instant;

- at least several coefficients of the representatives tensors of the time differential equation system are determined using a sum of products of at least two significant normal modes, at least one significant normal mode being a significant normal mode of a correlation matrix corresponding to the global density;

- at least one characteristic value comprises a global viscosity defined as a linear function of the global density;

- the Navier-Stokes type equations are defined by the global density, the velocity, the pressure and the global viscosity of the multiphase flow.

- the first reduced constant and the second reduced constant are defined using respectively the following expressions:

where

i is the local density of a first phase;

p ₂ is the local density of a second phase;

μ _χ is the local viscosity of the first phase; and

μ ₂ is the local viscosity of the second phase; - at least several coefficients of the representative tensors depend on the first reduced constant and at least several coefficients of the representative tensors depend on the second reduced constant;

- the flow area is a portion of a pipeline system comprising an inflow area and an outflow area;

- determining of a pressure drop in the portion of the pipeline system and / or a step of determining hold-up regions, the determining step comprising using at least one characteristic value of the multiphase flow, in particular a global density;

- the references parameters corresponding to at least one characteristic value are provided from a reference database established for the portion of the pipeline system;

- the total number of the significant normal modes is less than twenty and advantageously equals to ten.

The invention also concerns a computer program product having computer- executable instructions which, when carried out on a computer system, perform the method comprising at least some of features mentioned above.

The present invention concerns also a system for determining at least one characteristic value of a multiphase flow at each instant, the multiphase flow comprising a flow of hydrocarbon fluids in a flow area, each phase of the flow being characterized by a local density and a local viscosity of this phase;

the system comprising one or several calculation units able to:

- provide reference parameters corresponding to at least one characteristic value of the multiphase flow defined at each point of the flow area for different sampling instants;

- for at least said characteristic value, determine a temporal correlation matrix for all of the corresponding references parameters;

- for the or each determined correlation matrix:

+ determine the eigenvalues and eigenvectors of the determined correlation matrix;

+ determine space normal modes corresponding to a linear combination of eigenvectors of said correlation matrix and the corresponding references parameters;

+ determine a number of significant normal modes;

- determine coefficients of representative tensors associated to a time differential equation system, said equation system comprising time variables and resulting from the Galerkin projection of time discretized Navier-Stokes type equations on the significant normal modes; the coefficients of each representative tensor being determined using the significant normal modes of at least one correlation matrix, the coefficients of at least one representative tensor being further determined using reduced constants, each reduced constant being function of a combination of local densities and/or local viscosities of the phases of the multiphase flow;

- determine time variables of said equation system at a next instant using time variables at a current instant and the coefficients of the representatives tensors;

- determine at least one characteristic value at each instant using the temporal variables at this instant and the corresponding significant normal modes.

The invention will be better understood, upon reading of the following description, taken solely as an example, and made in reference to the following drawings, in which:

- Figure 1 is a schematic view of a characterization system implementing a method for determining at least one characteristic value of a multiphase flow, according to the invention;

- Figure 2 is a schematic view of a pipeline system in which at least one characterization value of the multiphase flos is able to be determined by the system of Figure 1 ;

- Figure 3 is a schematic view of a portion of the pipeline system of Figure 2; and

- Figure 4 is a general flow chart of steps for carrying out the method implemented by the system of Figure 1 .

Hereinafter, by "significantly equivalent", an equivalence relation with a maximum relative margin of 30% is understood.

The characterization system 10 illustrated on Figure 1 is configured for implementing a method for determining at least one characteristic value of a multiphase flow at each instant t.

Hereinafter, the wording "multiphase flow" refers to any fluid flow consisting of more than one phase flow, in particular a liquid flow and a gaseous flow.

According to a preferred embodiment of the invention, the multiphase flow is a hydrocarbon fluids flow in a pipeline system 1 1 illustrated on Figure 2.

The pipeline system 1 1 is used to transport hydrocarbon fluids such as different types of oils and gases in the oil and gas industry. Thus, the pipeline system 1 1 is configured to connect for example a wellhead to a production vessel or platform, or from an offshore platform to an onshore process plant.

The pipeline system 1 1 advantageously comprises a plurality of pipelines, valves, pumps, connection collectors, control units and other equipment known in the art. Particularly, the pipeline system 1 1 comprises a pipeline 1 3 extending along a central axis Y of a hydrocarbon fluids flow. The axis Y forms for example an angle γ comprised between 5° and 15° with a horizontal axis perpendicular to the gravitational field g.

An axial section of a portion 15 of the pipeline 13 is illustrated in details on Figure 3.

As illustrated on this figure, the portion 15 has a constant circular cross section of diameter D, and comprises an inflow area 1 7 and an outflow area 19. The diameter D is for example comprised between 2 cm and 160 cm.

The portion 1 5 further comprises a cylindrical wall 21 delimiting an interior part of the portion 1 5. Each point of the interior part of the portion 15 is defined by a position vector x which comprises for example three scalar components x, y, z.

The portion 1 5 further comprises a plurality of sampling points Pi to P _n distributed along the axis Y. Each sampling point P to P _n comprises a measurement apparatus (not illustrated) able to measure at least some characteristic values in relation with the multiphase flow as it will be explained later.

Particularly, each measurement apparatus is able to measure said characteristic values at each point x of a cross section of the portion 1 5 associated to the corresponding sampling point Pi to P _n at each predefined sampling instant t _a , a being a natural number comprised between 1 and M, M being a natural number comprised for example between 400 et 500.

The sampling instants t _a are comprised in a predetermined time interval [t ₀ , ... , T] where, for simplicity reasons, t ₀ is considered to be equal to zero.

Each sampling point P to P _n is connected to the characterization system 1 0 by data transmission means and able to transmit the measured characteristic values to the characterization system 1 0.

The hydrocarbon fluids flow in the portion 15 comprises notably a two-phase flow of first phase flow P, and a second phase flow P _M .

The first phase flow P, is for example a gas flow of hydrocarbon fluids in the pipeline system 1 1 . The second phase flow P _N is for example a liquid flow of hydrocarbon fluids in the pipeline system 1 1 .

In particular, in the example of Figure 3, the hydrocarbon fluids flow is a bubble flow wherein the first phase flow P, forms a plurality of gas bubbles distributed over the second liquid phase flow P _M .

Each phase Pi, P _N of the multiphase flow is characterized by local characteristic values such as a local density p ₁ , p ₂ and a local viscosity μ ₁ , μ ₂ - The local density p _t , p ₂ and the local viscosity μ ₁ , μ ₂ of each phase Pi, P _N correspond respectively to the density and to the dynamic viscosity of the corresponding phase considered independently of the other phase like in a monophasic flow.

The multiphase flow is characterized by global characteristic values defined at each point x of the portion 1 5 for each instant of time t.

The global characteristic values comprises notably a velocity vector U and a pressure P.

The velocity vector U is formed by three scalar components u, v, w. The component v corresponds to the multiphase flow velocity according to axis Y. The components u and w correspond respectively to the multiphase flow velocity according two perpendicular axes X and Z which are perpendicular to the axis Y.

The hydrocarbon fluids are able to cross the portion 1 5 from the inflow area 17 at an inflow pressure P ₁ to the outflow area 1 9 at an outflow pressure P ₂ . The difference AP = P ₂ — i between the outflow and inflow pressures is called hereinafter "pressure drop".

The components u, v, w of the velocity vector U and the pressure P of the multiphase flow are able to be measured at the sampling points Pi to P _n at each sampling instant t _a .

According to the invention, the multiphase flow is further characterized by a global density p and a global viscosity μ defined at each point x of the portion 1 5 for each instant t.

For a given instant t, the global density p at a given point x is equal to the local density p ₁ , p ₂ of one of the phases P,, P _N situated at the given point x and instant t.

The global viscosity μ is a function of the global density p. In particular, it is defined as a linear function of the global density as follows:

μ = c _t + pc ₂ ,

where c _x and c ₂ are respectively a first and a second reduced constants determined using the following expressions:

With reference to Figure 1 , the characterisation system 10 comprises a first calculation unit 31 , a second calculation unit 32 and a reference database 34 connected to both calculation units 31 and 32.

The reference database 34 is a memory device like a hard disc able to store a computer data issued from the first and the second calculation units 31 and 32. The first calculation unit 31 is a computer connected to each sampling point Pi to P _n and able to receive the measured characteristic values issued from the sampling points Pi to P _n at each sampling instant t _a .

The first calculation unit 31 is further able to analyse the received characteristic values in order to generate reference parameters of the multiphase flow for each sampling instant t _a .

The reference parameters of the multiphase flow at a given sampling instant t _a comprises a set of global characteristic values of the multiphase flow determined at each point x of the portion 15 for the given sampling instant t _a .

Such set of global characteristic values for a sampling instant t _a is called hereinafter a "snapshot" S _a of the multiphase flow at the sampling instant t _a . The number of snapshots S _a is thus equal to the total number M of the sampling instants t _a .

Finally, the first calculation unit 31 is able to store the generated reference parameters in the reference database 34.

The second calculation unit 32 is able to determine at least one characteristic value of the multiphase flow at each instant t using the references parameters storeed in the reference database 34 as it will be explained hereafter.

The first and the second calculation units 31 and 32 are for example situated in distinct geographical places and connected to the reference database 34 via a computer network such as Internet.

In a variant, the first and the second calculation units 31 and 32 form a unique calculation unit.

The method according to the invention will now be explained with reference to Figure 4 corresponding to a general flow chart of its steps.

During the initial step 1 10, the first calculation unit 31 receives from each sampling point Pi to P _n measured characteristic values of the multiphase flow.

Notably, the measured characteristic values comprise the components u, v, w of the velocity vector U and the pressure P measured at each point x of the cross section associated to the corresponding sampling point Pi to P _n at each sampling instant t _a comprised in the time interval [0, ... , T].

Then, the first calculation unit 31 determines using the received characteristic values, a snapshot S _a of the multiphase flow for each sampling instants t _a .

Each snapshot S _a comprises notably values of the components u, v, w of the velocity vector U, the pressure P and the global density p at each point x of the portion 15 for the corresponding sampling instant t _a . In order to determine the snapshots S _a , the first calculation unit 31 uses a CFD model known in the art. Advantageously, said CFD model is a one-dimensional model which allows the modelling of the global characteristic values of the multiphase flow at each point x of the portion 15 comprised between two neighbouring sampling points P to P _n .

Finally, the first calculation unit 31 stores the snapshots S _a in the reference database

34.

After at least one execution of the initial step 1 10, the following step 120 may be launched independently of the initial step 1 10. All of the following steps are performed by the second calculation unit 32.

During the step 120, the second calculation unit 32 acquires the snapshots S _a storeed in the reference database 34.

During the following step 130, the second calculation unit 32 constructs a temporal correlation matrix 0^β, θ"β respectively for the global density p and the velocity vector U using the snapshots S _a , α, β being indices of the correlations matrices varying between 1 and M.

The matrices C _a ^p _B , C are defined using the following expressions:

where:

- ίβ is a sampling instant distinct from the instant t _a or coinciding with it;

- C ^p (t _a , t ), C ^u (t _a , t ) are operators corresponding to the matrices C^ C^ ;

- Ω is an integration area equal to the portion 15 of the pipeline 13 and comprising a plurality of points x.

The coefficients of the correlation matrices C^, c"p are obtained by a numerical integration according to a numerical method known in the art.

During the following step 140, the second calculation unit 32 determines the eigenvalues /if, /i and the eigenvectors b α _έ respectively for the matrices C^, c"p by resolving respectively the following eigenvalue problems:

where i is a natural number varying between 1 and M, and the eigenvalues λ^,λ" are ordered according to a decreasing order. This means that:

^Λ 1 ^{≥ Λ} 2 ^{≥ Λ} Μ·

— ^2 — — ^M.

These problems are resolved using numerical methods of eigenvalue problem resolution known in the art.

It should be noted that the eigenvectors b α _έ are functions of the time t and the resolution of the eigenvalue problems mentioned above allows the determination of these functions at each sampling instant t _a .

During the following step 143, the second calculation unit 32 determines space normal modes of, <I> respectively of the global density p and the velocity vector U functions.

The space normal modes Of, are functions of the point x and correspond to a linear combination of the corresponding eigenvectors b α _έ and the corresponding global characteristic values at the sampling instants t _a . Thus, the space normal modes Φ?, Φ" are defined as follows:

M

M

*i = _{j i} (t _a )U{ x, t _a ).

a=l

The space normal modes of , Φ" present spatial bases respectively for the functions p(x, t) and U(x, t).

During the following step 148, the second calculation unit 32 determines a number Npof significant normal modes.

The significant normal modes correspond to space normal modes Φ?, Φ" taken in consideration for representing the functions p(x, t) and U(x, t) in the corresponding bases.

Hereinafter, the significant normal modes will also be denoted by of, Φ" . But in this last case, the index i varies between 1 and the number N _P . The number N _P is less than the total number M of the normal modes of, Φ" .

In order to determine the number N _P , according to an embodiment of the invention, the second calculation unit 32 determines a relative energy E _rel (i) associated to each normal mode Φ" of the velocity function U(x, t). The relative energy E _rel (i) is defined as follows: k

Ls=l ^A s

The number N _P is defined using the following expression :

N _P = max{i, 1≤ i≤M \ E _rel (i)≥ E ₀ ],

where E ₀ is a predetermined threshold. Thus, the number N _P corresponds to the maximum index i for which the value of the associated relative energy E _rel (i) is higher or equal than a predetermined threshold E ₀ .

According to another embodiment of the invention, the number N _P is determined in the same way but using a relative energy E _rel (i) associated to each normal mode of of the density function p(x, t).

According to another embodiment of the invention, the second calculation unit 32 makes the number N _P equal to a predetermined number which is equal for example to ten.

During the following step 1 50, the second calculation unit 32 determines coefficients of representative tensors M _ijk , Q _ijk i, Ρ^, Ό^, Ε^, β^, Ν^ and R _ijk .

The coefficients of the representative tensors M _ijk , Q _ijkl , P _i7 , D _i7 , £ _i7fc , G _i7 , N _i; and R _ijk correspond to time-independent coefficients of a time differential equation system DS as defined below.

The time differential equation system DS results from a Galerkin projection of

Navier-Stokes time discretized equations on the significant normal modes of, Φ" of the density p(x, t) and the velocity U(x, t) functions.

In particular, the Navier-Stokes equations for the multiphase flow in the portion 15 of the pipeline 1 1 are defined as follows:

dU

p ~dt ^{+ ρ (} ~ ^υ - ^{ν) ί/ = _ V + μΑ υ +}

dp where:

p and μ are respectively the global density and the global viscosity as defined above;

g is the gravity field;

V is the gradient operator known in the art; and

Δ is the Laplace operator known in the art.

A time discretized version of the Navier-Stokes equations corresponding to a time explicit numerical scheme known in the art, has the following form:

ρη _υ η+1 _{= ρ} η _υ η _{+ Δί} (_ _ρ η( _ί/ η _{ν) ί/} η _ _vp n ₊ )

pn+ ^l _{= p} n _{+ Δί} (_(ί/". V)p ⁿ ) = 0, (2) where At is a time discretization step.

Hereinafter, the values with the index n are determined at a current instant t ⁿ , and the values with the index n + 1 are determined at a next instant t ⁿ⁺¹ corresponding to the current instant t ⁿ increased by the value of time discretization step At.

Of course, other time numerical schemes can be applied to the Navier-Stokes equations so as other discretized versions of these equations may be obtained.

The Galerkin projection of the time discretized Navier-Stokes equations is obtained by multiplying the equation (1 ) by each significant normal mode Φ" and by integrating each member of this equation over the portion 15, and by multiplying the equation (2) by each significant normal mode of and by integrating each member of this equation over the portion 15.

Finally, the time differential equation system DS corresponds to the Galerkin projection of the time discretized Navier-Stokes equations in which the functions p, U of the global density and of the velocity vector are replaced by the corresponding expressions in the bases of the significant normal modes of, Φ" using the following expressions:

P

i=i

and the function P of the pressure is replaced using its representation in a basis of significant normal modes Of defined via the eigenvectors which means:

Thus, time differential equation system DS is represented via the representative tensors M _ijk , Q _ijkl , P _t} , D _tj , E , G _tj , N _tj and R as follows:

+ At{-Q _ijkl b?altf - P _tJ af + D _tJ af + E _i]k a ⁿ _k bf + G _tj b ),

Nt _j b^ ¹ = Nt _j bJ ¹ + At{-R _ijk a ⁿ _k b ),

where the coefficients of the representative tensors are defined as follows:

Mij _k = J ^" _n ^φ _¾ Φ/ ⁷ Φ"άίΙ : mass tensor for U ;

Qi _jk i = f _a (Φ ■ V)O Φ"άΩ. : advection tensor for U;

P _tj = J ^" _n V( ) Φ"άΩ. : pressure tensor for U;

Di _j = - J _n ο-^Φ" V dn + J _r c _± d o /dn Φ?άΓ : first diffusion tensor for U; _¾ νΦ^Φ _; ^ρ V dn + / _Γ c ₂ d Φ%/δη Φ^Φ^ dT

: second diffusion tensor for U; : gravitational tensor for U;

: mass tensor for p;

J? _{£ fc} = /^- V)o o dn : advection tensor for p; and

Γ denotes the border of the portion 15.

The coefficients of the representative tensors are determined using numerical methods of integration over the portion 15, known in the art.

During the following step 160, the second calculation unit 32 launches a time loop for determine the eigenvectors α _έ (ί), ί> _έ (ί) at each instant t different from the sampling instants t _a .

Thus, each next value of the eigenvectors α _έ (ί), ί> _έ (ί) at the next instant t ⁿ⁺¹ is determined using the current values of the eigenvectors α _έ (ί), Ζ? _έ (ί) at the current instant t ⁿ according to the time differential equation system DS.

During the following step 170, the second calculation unit 32 determines at a given point x of the portion 15 of the pipeline 13 at least some global characteristic values of the multiphase flow as the density p, the velocity vector U or the pressure P for a given instant t.

To this end, the second calculation unit 32 uses the eigenvectors α _έ (ί), Ζ? _έ (ί) determined on the previous step 160 and the following expressions:

During an optional step 180, the second calculation unit 32 determines the pressure drop Δ in the portion 15 at a given instant t using the inflow pressure value P ₁ and the outflow pressure value P ₂ respectively in the inflow area 17 and the outflow area 19. The values P _t and P ₂ are determined using the values of the pressure function P(x, t) determined on step 170.

During an optional step 185, the second calculation unit 32 determines gas hold-up regions in the multiphase flow over the portion 15 by figuring out regions with relatively low values of the global density p. Preferably, a region of low density has a global density comprised between 600kg/m3 and 830kg/m3 and advantageously smaller than 800kg/m3.

The method according to the invention thus allows the determination of global characteristic values of a multiphase flow such as pressure or velocity vector. The determination of these values allows the determination of the local characteristic values in relation with each phase forming the multiphase flow such as local density or viscosity.

The characteristic values of the multiphase flow further allows the determination of other parameters of a hydrocarbon fluids flow in a pipeline system such as pressure drops or hold-up regions.

The method according to the invention is a three-dimensional method allowing the determination of all characteristic values of the multiphase flow in a very accurate way. The method is particularly fast and able to determine the characteristic values of a multiphase flow significantly as fast as a classical one-dimensional method event with a limited quantity of the computational resources.

The use of a global viscosity and a global density allows a quite simple mathematical format to be used, while still maintaining the possibility of dissimilarity between the different phases.

According to another embodiment of the invention, the initial step 1 1 0 of determining the snapshots S _k is performed in any other suitable way.

Thus, for example, the snapshots S _k are determined using two- or three-dimensional

CFD methods and/or using a modelling of a turbulent flow.

In addition or in a variant, the snapshots S _k for the portion 1 5 of the pipeline system 1 1 are determined only once during a test phase of exploitation of the pipeline system 1 1 .

In addition or in a variant, the snapshots S _k for the portion 1 5 of the pipeline system 1 1 are determined using other suitable ways of modelling of the pipeline system 1 1 .

In a variant, in order to construct the snapshots S _k , at least some characteristic values are simulated in the sampling points to P _n using modelling methods known in the art.

Hence, this allows stocking the snapshots S _k in the reference database 34 and its using when a determination of at least one characteristic value of the multiphase flow is needed. So, this value can be determined in a very simple and fast way.