In the present disclosure, reference shall be made, in particular, to systems that carry the hydrocarbons extracted from wells to the inlet of the distribution network.

However, the simulation method of the present invention can also be applied to the distribution network.

Currently, simulating the thermo-fluid dynamic behavior of multiphase hydrocarbons through thermo-fluid dynamic simulators that solve the Navier-Stokes equations by finite volume discretization methods is known.

These finite volume thermo-fluid dynamic simulators have good reliability, but they also present a sizable computational cost which determines high simulation times, which grow as complexity and the dimensions of the analyzed system grow.

Simulation times are critical during the studies for the development of a hydrocarbon production and distribution system, design steps in which it can be extremely useful to very rapidly identify potentially critical situations for the productive life of the system. An object of the present invention is to overcome the aforementioned drawbacks and in particular to devise a method for simulating the thermo-fluid dynamic behavior of multiphase fluids in a hydrocarbons production and transport system that is able to obtain reliable results, while entailing requiring shorter simulation times than prior art finite volume thermo-fluid dynamic simulators .

This and other objects according to the present invention are achieved by a method for simulating the thermo-fluid dynamic behavior of multiphase fluids in a hydrocarbons production and transport system according to claim 1.

Additional characteristics of the method for simulating the thermo-fluid dynamic behavior of multiphase fluids in a hydrocarbons production and transport system are the subject of the dependent claims.

The characteristics and advantages of a method for simulating the thermo-fluid dynamic behavior of multiphase fluids in a hydrocarbons production and transport system according to the present invention shall become more readily apparent from the following exemplifying and non-limiting description, referred to the accompanying schematic drawings in which:

- figure la is a schematic view representing a component block of a hydrocarbons production and transport system modeled according to a simplified conduit model;

figure lb is a schematic view representing a component block of a hydrocarbons production and transport system modeled according to a simplified valve model; figure lc is a schematic view representing a component block of a hydrocarbons production and transport system modeled according to a simplified reservoir model;

- figure Id is a schematic view representing a component block of a hydrocarbons production and transport system modeled according to a simplified separator model;

figure 2 is a schematic representation of a hydrocarbons production and transport system to be simulated;

- figures 3a, 3b, 3c and 3d are four elements of an oriented graph that represent respectively a conduit, a valve, a reservoir and a separator;

- figure 4 is an oriented graph that represents the series connection of a reservoir, a well and a valve; figure 5 is a schematic view representative of a conduit in slug flow regime;

- figure 6 is a flow chart that represents a method for the simulation of the thermo-fluid dynamic behavior of multiphase fluids in a hydrocarbons production and transport system according to the present invention; figure 7 is a schematic block diagram that illustrates a learning step of simplified analytical mathematical models of the corresponding component blocks of the system of figure 2;

- figure 8 is a schematic block diagram representing a resolution step in the method for the simulation of the thermo-fluid dynamic behavior of multiphase fluids in a hydrocarbons production and transport system of figure 4;

- figure 9 shows a plurality of charts that show the trend over time of pressure (PT) , temperature (TM) , mass flow rate of liquid (GLT) and gas (GG) , calculated at a node of the system of figure 2 by the simulation method according to the present invention (solid line) and by a prior art finite volume thermo-fluid dynamic simulator (dashed line) ; in particular, the charts on the left are obtained before the learning step of the simplified analytical mathematical models, while the charts on the right are obtained after the learning step of the simplified models.

With reference to the figures, a simulation method is shown for simulating the thermo-fluid dynamic behavior of multiphase fluids in a hydrocarbons production and transport system, indicated in its entirety with the number 100.

Said simulation method 100 comprises an initial step 110 in which the transport system is outlined as a plurality of interconnected component blocks thus obtaining a schematic representation such as the one shown in figure 2. Each component block is then modeled 120 according to a simplified analytical mathematical model selected from the group of models comprising at least one conduit model, a valve model, a reservoir model, a separator model.

The system shown in figure 2, for example, comprises two reservoirs Gl and G2 with related wells PI and P2 and wellhead valves VI and V2 ; from said valves originate two conduits P3 and P4 that join a third conduit P5, which ends in a constant pressure node (e.g. the inlet of a separator) .

It is stressed that each of the two wells PI and P2 of the two reservoirs was modeled according to the simplified conduit model.

Each of these simplified analytical mathematical models is represented by a plurality of constitutive equations adapted to describe the thermo-fluid dynamic behavior of the corresponding component block. In particular, the constitutive equations describe the evolution of a plurality of variables between inlet and outlet of the individual component block, as well as the evolution of these variables over time; the simulation method of the present patent is able to describe the characteristic transients in the operational scenarios described above .

The aforesaid variables are, for example, the flow rates of the different phases, the pressures, the temperatures and so on.

Preferably, the modeling step 120 further comprises the step in which for each simplified analytical mathematical model a plurality of corrective coefficients is applied to the constitutive equations, where the corrective coefficients are estimated so as to adapt the results obtained from the simplified analytical mathematical model to reference data. These reference data can be, in particular, derived from actual field measurements or from known finite volume thermo-fluid dynamic simulators.

The corrective coefficients are defined in such a way as to assume a value equal to 1 in case of perfect agreement between the simplified analytical mathematical model and the reference date. Therefore, greater discrepancies or inadequacies of the simplified model relative to the reference date are expressed by the progressive departure of the corrective coefficients from the unitary value.

In modelling the discrete elements of the system, it must be considered that the conduit element and the separator element are dynamic elements, i.e. provided with memory, whereas the reservoir element and the valve element are static, i.e. without memory.

In the present description, the subscript L refers to a liquid, the subscript G refers to a gas, the subscript M refers to a liquid-gas mixture, the subscript IN refers to an inlet, the subscript OUT refers to an outlet, the subscript UP refers to the upstream side of a valve, the subscript DO refers to the downstream side of a valve, the subscript RES refers to a reservoir, the subscript W refers to a well, the subscript V refers to a valve, the subscript P refers to a conduit ("pipeline"), the subscript SEP refers to a separator. Preferably, the simplified analytical mathematical conduit model can comprise the following constitutive equations :

- two mass conservation equations, one for the liquid phase and one for the gaseous phase: m _L =- [m _L v _L0UT -m _LO v _LJN )-y/ _G

where m _L indicates the mass of liquid per unit of volume, m _G indicates the mass of gas per unit of volume, m _L o indicates the mass of liquid per unit of volume in the previous conduit, m _G o indicates the mass of gas per unit of volume in the previous conduit, AL indicates the length of the conduit, v _L , i _N and V _L , ₀ UT indicate the velocity of the liquid respectively at the inlet and at the outlet, v _G , IN and V _G , ₀ UT indicate the velocity of the gas respectively at the inlet and at the outlet, i|r _G indicates the mass flow rate per unit of volume that changes from the liquid phase to the gaseous phase.

a total momentum conservation equation (describes both phases) :

A

\ ^M L ^V L,OUT ^M L ^V L,IN + ^M G0 ^V G,OUT ^M G0 ^V G,IN )-—{POUT - PIN ) ^~ A T - AR where G _ou t indicates the mass flow rate out, A indicates the cross section area of the conduit, AL indicates the length of the conduit, v _L , i _N and V _L , ₀ UT indicate the velocity of the liquid respectively in and out, v _G , IN and V _G ,OUT indicate the velocity of the gas respectively in and out, m _L indicates the mass of liquid per unit of volume, m _G indicates the mass of gas per unit of volume, m _L o indicates the mass of liquid per unit of volume in the previous conduit, m _G o indicates the mass of gas per unit of volume in the previous conduit, Γ indicates the head loss per unit of length (which depends on the flow regime) , R indicates the friction loss per unit of length (which depends on the flow regime) , pi _n and p _ou t indicate the pressure respective at the inlet and at the outlet.

a total energy E conservation equation (describes both hases :

where m _L indicates the mass of liquid per unit of volume, m _G indicates the mass of gas per unit of volume, m _L o indicates the mass of liquid per unit of volume in the previous conduit, m _G o indicates the mass of gas per unit of volume in the previous conduit, AL indicates the length of the conduit, v _L , i _N and V _L , _0UT indicate the velocity of the liquid respectively in and out, v _G , _IN and V _G , _OUT indicate the velocity of the gas respectively in and out, A indicates the cross section area of the conduit, AL indicates the length of the conduit, S indicates the perimeter of the section of the conduit, U indicates the relative heat transfer coefficient of the walls of the conduit (including the insulation layers) , T _out and T _ex t indicate the temperature respectively at the outlet and external, g indicates the acceleration due to gravity, h _L and h _G indicate the specific enthalpy, respectively, of the liquid and of the gas, h _L o and h _G o indicate the specific enthalpy in the previous conduit respectively of the liquid of the gas, zin and z _out indicate, respectively, the elevation of the inlet and of the outlet of the conduit relative to a reference level;

- an equation that describes the evolution of pressure (describes both phases) obtained from the combination of the two mass conservation e uations:

where m _L and m _G indicate the mass per unit of volume respectively of liquid and of gas, m _L o and m _G o indicate the mass per unit of volume in the previous conduit respectively of liquid and of gas, V _L , _IN and V _L , _0UT indicate the velocity of the liquid respectively in and out, v _G , _IN and V _G , _0UT indicate the velocity of the gas respectively in and out, p indicates pressure, AL indicates the length of the conduit, p _L and p _G indicate the density respectively of the liquid and of the gas, x _G indicates the quality of the gas, ψ indicates a flow rate per unit of volume obtained from the mass flow rate per unit of volume ψ defined above, _L indicates volume percent of liquid.

Using a single momentum conservation equation instead of two distinct ones for each phase, it is necessary to introduce an algebraic relationship to take into account the difference in velocity ("slip") between the two phases. Therefore, in addition to the aforementioned equations for each flow regime (stratified, dispersed bubble, slug, etc.), the simplified conduit model can also comprise the following equations:

- an equation that defines the terms R that take into account, in the momentum conservation equation, the frictions of the fluid with the walls of the conduit and between the different phases (R) ;

- an equation that defines the terms T that take into account, in the momentum conservation equation, hydraulic head losses;

- a slip equation.

For example, for the stratified regime the following equations are considered:

equation defining the friction loss per unit of length R:

where R _L and R _G indicate the friction losses respectively of the liquid phase and of the gaseous phase, A indicates the cross section area of the conduit, f _L and f _G indicate the friction coefficients respectively of liquid and of gas, p _L and p _G indicate respectively the densities of liquid and of gas, v _L and v _G indicate respectively the velocities of liquid and of gas, S _L and S _G indicate the parts of circumference of the section of the conduit that are "wet" respectively by liquid and by gas;

- equation defining the head loss per unit of length Γ r=(m _L +m _G )gsinQ

where m _L and m _G indicate the mass per unit of volume respectively of liquid and of gas, g indicates the acceleration due to gravity, Θ indicates the inclination of the conduit relative to the horizontal;

- slip equation: a _L {l - _L ) l -a _L a _L

where R _L and R _G indicate the friction losses respectively of the liquid and gaseous phase, Ri indicates the friction loss between the two phases, g indicates the acceleration due to gravity, Θ indicates the inclination of the conduit relative to the horizontal, p _L and p _G indicate the density respectively of the liquid and of the gas, _L indicates the volume percent of liquid.

For the bubbly regime, the following equations are considered :

equation defining the friction loss per unit of length R: ^R = —/MPMVM

where A indicates the cross section area of the conduit, f _M indicates the friction coefficient of the liquid-gas mixed phase, p _M indicates the density of the liquid-gas mixture, v _M indicates the velocity of the fluid, D indicates the diameter of the conduit;

- equation defining the head loss per unit of length Γ r=(m _L +m _G )gsir

where m _L and m _G indicate the mass per unit of volume respectively of liquid and of gas, g indicates the acceleration due to gravity, Θ indicates the inclination of the conduit relative to the horizontal;

- slip equation:

v _L and v _G indicate the velocities respectively of liquid and of gas, g indicates the acceleration due to gravity, Θ indicates the inclination of the conduit relative to the horizontal, p _L and p _G indicate the density respectively of the liquid and of the gas, σ indicates the surface tension of the liquid.

For the slug regime, the following equations are considered :

equation defining the friction loss per unit of length R:

1 £ 1 \ £

LA " s+T " s+T where R _s and R _T indicate the friction losses in the "slug" section and in the "Taylor bubble" section of the slug unit, A indicates the section of the conduit, f _s and f _GT indicates the friction coefficients of the fluid in the "slug" and "Taylor bubble" sections, p _s and p _G indicate the densities of the fluid respectively in the "slug" and in the gas section, v _M and V _GT indicate the velocities of the fluid and respectively in the "slug" section and of the gas contained in the "Taylor bubble", S and S _GT indicate respectively the circumference of the section of the conduit and the part of circumference "wet" by the gas of the "Taylor bubble", l _s the length of the "slug" section and 1 _T the length of the "Taylor bubble" section, 1 _S+T indicates the total length of the slug unit;

- equation defining the head loss er unit of length Γ

where g indicates the acceleration due to gravity, Θ indicates, p _s and p _T indicate the density respectively in "slug" regime and in "Taylor bubbly" regime, l _s the length of the "slug" section, 1 _T indicates the length of the "Taylor bubbly" section, 1 _s+T indicates the total length of the slug unit;

- slip equation:

^V _G ^~ C _QT v _M —v _QT =0

where C _OT is a flow distribution coefficient and V _OT is the velocity of a gas bubble that rises along a conduit of stagnating liquid (v _L = 0) .

As corrective coefficients can be selected, for example, multiplicative coefficients that correct:

- the slip equation;

- the friction component;

- the derivative of gas quality with respect to pressure .

In particular, for the stratified regime, the λ ₃ coefficient is inserted in the slip equation as follows :

where R _L and R _G indicate the friction losses respectively of the liquid and gaseous phase, Ri indicates the friction loss between the two phases, _L and _G indicate the volume percentages respectively of liquid and of gas, g indicates the acceleration due to gravity, Θ indicates the inclination of the conduit relative to the horizontal, p _L and p _G indicate the density respectively of the liquid and of the gas.

For the bubbly regime, the λ ₃ coefficient is inserted in the slip equation as follows:

v _L and v _G indicate the velocities respectively of liquid and of gas, g indicates the acceleration due to gravity, Θ indicates the inclination of the conduit relative to the horizontal, p _L and p _G indicate the density respectively of the liquid and of the gas.

For the slug regime, the λ ₃ coefficient is inserted in the slip equation as follows:

^V _G ^~ C _0T v _M —v _0T =0

where v _G is the velocity of the gas, v _M the velocity of the fluid, C _OT is a flow distribution coefficient and V _OT is the velocity of a gas bubble that rises along a conduit of stagnating liquid (vL = 0) . The corrective coefficient _f is multiplied times the friction coefficient of the liquid phase f _L in case of stratified regime, times the friction coefficient of the mixed liquid-gas phase f _M in case of "bubbly", times the friction coefficient of the "slug" phase fs regime in case of "slug" regime.

The corrective coefficient X _dx , _P is multiplied times the derivative of gas quality with respect to pressure dx.

dp

The simplified conduit model can thus be written in state-space representation: where u _P is the vector of the "forcing" variables, e.g.:

- outlet pressure P _OUT

- inlet temperature ΤΊ _Ν ;

- mass flow rate of the liquid at the inlet G _L , _I N;

- mass flow rate of the gas at the inlet G _G , _I N,- x _P is the vector of the state variables, e.g.:

- input pressure pi _N;

- internal energy of the fluid in the conduit element per unit of volume E;

- total mass flow rate at the outlet G _0UT ;

- mass of liquid in the conduit element per unit of volume m _L ;

- mass of gas in the conduit element per unit of volume m _G ;

y _P is the variable output vector, e.g.:

- inlet pressure P _I N/

- outlet temperature Τ ₀ υτ ; - mass flow rate of the liquid at the outlet G _L , OUT ;

- mass flow rate of the gas at the outlet G _G , OUT ;

- average hold-up of the conduit or volume percent of liquid _L , i.e. the ratio between the volume of liquid contained and the total volume of the conduit.

λ _Ρ is the vector containing the aforementioned corrective multiplicative coefficients.

Preferably, the simplified analytical mathematical valve model can comprise the following constitutive equations :

- an equation that ties the total mass flow rate to the pressure drop across the valve; one possible model is the Perkins model, known to persons skilled in the art;

G - fv(Pm )- A _m gPupP

where A _TH is the section area of the throat of the valve, A _UP is the section of the valve inlet, p _T H and p _up indicate the pressure respectively at the throat and at the inlet of the valve, x _G indicates gas quality, n indicates the polytropic expansion exponent, Φ and n are parameters related tothe mass percentages of liquid and gas, G indicates the total mass flow rate that traverses the valve and p indicates the density of the fluid;

an equation that ties the mass flow rate of gas downstream and the mass flow rate of gas upstream; G

dx

I P , ^Gg,UP (™L,UP + m _{G D0} )

ΰΡ I DO ^A UP ^M L,UP

where G _G , UP and G _G , DO indicate the mass flow rates of gas respectively at the inlet and at the outlet of the valve, A _UP is the inlet section area of the valve, <^Q indicates the mass flow rate of gas that becomes liquid per unit of surface area, m _L , _UP and m _L , _D0 indicates the liquid mass per unit of volume respectively at the inlet and at the outlet of the valve, p indicates pressure, Δρ indicates the pressure difference between the inlet and the outlet of the valve;

- an equation that describes the temperature interval between upstream and downstream T _D0 .

TDO ^X G,UP^G,UP ⁺ (l ^{_ X} G,Up)^L,UP ^{= X} G,DJ' ¹ G,DO ^{+ _ X} G,Do)^L,DO where XG, UP and XG, DO indicate gas quality respectively at the inlet and at the outlet of the valve, h _L , _UP and h _L , _D0 indicate the specific enthalpy of the liquid respectively at the inlet and at the outlet of the valve, h _G , _UP and h _G , _D0 indicate the specific enthalpy of the gas respectively at the inlet and at the outlet of the valve.

As corrective coefficients can be selected, for example, multiplicative coefficients C _D , _, v r that correct :

- the total mass flow rate G (in fact, the corrective coefficient is the "discharge coefficient", known to the persons skilled in the art) ; in this case, the following correspondence applies G→ C _D G ;

gas quality x _G , a function of pressure p and of temperature T; in this case, the following correspondence applies x _G (p,T)—λ _χν χ ₀ (ρ,Τ) ·

the derivative of gas quality with respect to pressure xG, a function of pressure p and of temperature T; in this case, the following correspondence applies (?,TJ—>· A _{a F} p,T);

dp ' dp

The simplified valve model can then be written in representation in the following way:

where u _v is the vector of the "forcing" variables, e.g.:

- pressure loss between upstream and downstream Δρ;

- upstream pressure P _UP ;

- upstream temperature T _UP ;

- mass flow rate of the upstream gas G _G , _UP

- liquid mass in the upstream conduit element (or in the reservoir) per unit of volume m _L , _UP;

gas mass in the upstream pipeline (or in the reservoir) per unit of volume m _G , up ;

- percentage of valve opening k _v ,

yv is the variable output vector, e.g.:

- total mass flow rate G

- mass flow rate of the downstream gas G _G , _DO

- downstream temperature T _D0 .

λ _ν is the vector containing the aforementioned corrective coefficients.

Preferably, the simplified analytical mathematical reservoir model can comprise the following constitutive equations :

- an equation that ties the total mass flow rate G to the difference between static pressure and inflow pressure :

^{G= (} 1>IPIIPRES-PW)

where Φι _ΡΚ ( ) is a function that represents the IPR curve known to the persons skilled in the art, p _RE s and p _w indicates respectively the pressure in the reservoir and in the well;

- an equation that describes the mass flow rate of gas G _G as a function of the total mass flow rate G:

where RS,RES indicates the mass flow rate percentage of gas in the reservoir;

- an equation that describes the temperature interval betwee the reservoir and the inlet of the well ΤΊ _Ν .

where XG,RES and x _G , indicates gas quality respectively in the reservoir and in the well, h _L , _RES and h _L , _w indicate the specific enthalpy of the liquid respectively in the reservoir and in the well, h _G , _RES and h _G , _w indicate the specific enthalpy of the gas respectively in the reservoir and in the well.

As corrective coefficients can be selected, for example, a single multiplicative coefficient _χ ^ which corrects gas quality x _G as a function of pressure p and of temperature T; in this case, the following correspondence applies x _G (p,T)—λ _χΚ χ _β (ρ,Τ) .

The simplified reservoir model can then be written in representation in the following way:

where u _R is the vector of the "forcing" variables, e.g.: - static pressure of the reservoir p _RE s,

- inflow pressure p _w ,

- temperature of the reservoir T _RES ,

y _R is the variable output vector, e.g.:

- total mass flow rate G

- mass flow rate of gas G _G

- temperature at the inlet of the well T _w ,

λ _κ is the aforementioned corrective coefficient.

Preferably, the simplified analytical mathematical separator model can comprise the following constitutive equations :

- an equation that ties the mass flow rate of liquid at the outlet G _L , _OUT to the mass flow rate of liquid at the inlet G _L ,IN: where G _DRAIN indicates the maximum drain rate of the separator;

- an equation that ties the mass flow rate of gas at the outlet G _G , _OUT to the mass flow rate of gas at the inlet G _G , IN:

- an equation that describes the volume of liquid V _L inside the separator: V — ^ ^L _' ^{IN ~} ^ ^L > ^0UT

where G _L , _I N and G _L , _OUT are respectively the mass flow rate of liquid flowing in and out and p _L is the density of the liquid at the separation conditions of temperature SEP and pressure PSEP- As corrective coefficients can be selected, for example, a single multiplicative coefficient λ s _ep which corrects the density of the liquid gas p _L as a function of the pressure and of the temperature at the separator.

The simplified separator model can then be written in representation in the following way:

-½ ⁼ fsi ^x s ^■ > ^u s ^{■ '} > ) where u _s is the vector of the "forcing" variables, e.g.: - separation pressure P _SEP ,

- separation temperature T _SEP ,

- mass flow rate of liquid flowing in G _L , _IN ,

- mass flow rate of gas flowing in G _G , _I N,

y _s is the variable output vector, e.g.:

- mass flow rate of liquid flowing out G _L , _OUT/

- mass flow rate of gas flowing out G _G , _OUT ,

- internal volume of liquid V _L ,

x _s is a state variable coinciding with one of the outputs, i.e. the volume of liquid V _L ,

λ ₃ is the aforementioned corrective coefficient.

The constitutive equations of the above simplified analytical mathematical models are derived from fluid dynamics law that are known in themselves; the parameters and the variables contained in these equations can be calculated in a known manner.

As stated above, the corrective coefficients of each simplified model are estimated in such a way as to adapt the results obtained by the model to reference date .

For this purposes, the step of modeling the component blocks 120 preferably comprises a learning step 200 in which the corrective coefficients are estimated for each simplified model for at least one stationary or transient flow regime (stratified, dispersed bubbly, slug and so on) and/or for at least one stationary or transient heat regime.

This learning step 200 can be carried out with different mathematical methods for the minimization of the discrepancy between the data obtained from the simplified model and the reference data.

For example, the least squares method can be used, considering as reference data the values of interest of pressure, temperature, flow rates, etc. measured in the field or calculated by means of finite volume thermo- fluid dynamic simulators.

Let y _x and u _x be the vectors of the reference data to be used respectively as reference output variables and forcing variables of the simplified model X.

Consider the case in which the simplified model X is the one relating to a conduit; in this case x _P is the vector of state variables of reference, which can be calculated starting from y _P and u _P with an equation

The system to be solved with the least squares technique to estimate the corrective coefficients λ _Ρ of the conduit model is for example

where to is an instant in stationary state conditions. Solving the system indicated above is equivalent to imposing that, at the instant to, the state x _P is stationary and that the output variables of the simplified model are equal to the reference variables (measured or simulated with a reference model) .

Let us consider the case in which the simplified model X is the one relating to a valve.

The system to be solved with the least squares technique to estimate the corrective coefficients λ _ν of the valve model is for exam le

where to...t _N is the sequence of instants belonging to the time interval of interest. Solving the system indicated above is equivalent to imposing that the output variables of the simplified model are equal to the reference variables (measured or simulated) .

Let us consider the case in which the simplified model X is the one relating to a reservoir.

The system to be solved with the least squares technique to estimate the corrective coefficients λ _κ of the valve model is for exam le

where to...t _N is the sequence of instants belonging to the time interval of interest. Solving the system indicated above is equivalent to imposing that the output variables of the simplified model are equal to the reference variables (measured or simulated) .

Consider the case in which the simplified model X is the one relating to a separator.

The system to be solved with the least squares technique to estimate the corrective coefficients λ ₃ of the separator model is for example where to...t _N is the sequence of instants belonging to the time interval of interest. Solving the system indicated above is equivalent to imposing that the output variables of the simplified model are equal to the reference variables (measured or simulated) .

Figure 7 is a schematic representation of the learning step 200 relating to the simplified analytical mathematical models used to describe the system shown in figure 2. In particular, in an exemplifying manner, figure 7 relates to a learning step 200 carried out for the slug regime. This learning step 200 comprises considering 210 a particular operating condition, e.g. the one represented by the charts 300. This operating conditions corresponds to a sudden reduction in flow rate, obtained by partially closing both wellhead valves (turn down) , in accordance with the curves of the charts 300. For this operating condition, the reference data for learning are determined 220 performing a simulation by means of a known finite volume thermo-fluid dynamics simulator or carrying out measurements in the field. Thereafter, the equation systems are determined 230 and solved 240 to estimate the corrective coefficients λ _Ρ of the conduit model, the corrective coefficients λ _ν of the valve model, the corrective coefficients X _R of the reservoir model on the basis of the previously obtained reference data. In particular, the equation systems for estimating the corrective coefficients λ _Ρ of the conduit model are as many as the component blocks modeled as a conduit; in the example in figure 7 they are five. Similarly, the equation systems for estimating the corrective coefficients λ _ν of the valve model are as many as are the component blocks modelled as a valve, two in the example in figure 7; lastly, the equation systems for estimating the corrective coefficients λ _κ of the valve model are as many as are the component blocks modeled as a valve, two in the example in figure 7. Solving 240 the aforesaid equation systems equations, the set of corrective coefficients λ _Ρ , λ _ν , λ _κ are obtained.

In any case, the simulation method can comprise estimating the corrective coefficients for each flow and/or thermal regime, selecting as a reference the one corresponding to corrective coefficients that more closely approach the unit.

As stated above, a production and transport system can be described by a schematic representation like the one in figure 2.

The simulation method 100, according to the present invention, advantageously comprises a step in which an oriented graph is generated 130 on the basic of the aforesaid schematic representation. An example of oriented graph is shown in figure 4 and it relates to a series connection between a reservoir, a well and a valve .

The graph presents a plurality of edges j, a plurality of nodes I and a plurality of loops k. The loop are the cyclical paths of the graph, or a set of adjacent edges in which the last edge ends in the node from which the first one starts. It can be demonstrated that their number is M-N+l, where M is the number of edges and N the number of nodes.

To each node I is associated a set of scalar values, e.g. values of pressure, temperature or mass of liquid or gas per unit of volume.

To each edge j is associated a set of pairs of scalar values :

- a flow ξ-j, which can be a value of mass flow rate; the orientation of the individual edge j matches the direction of the flow ξ-,;

- a gradient ψ-j , obtained from the difference between the corresponding scalar values of the connected nodes, e.g. a difference in pressure, in temperature, and so on .

These quantities are not equal in each point of the conduit and of the valve, therefore they can only be associated to inlets and outlets thereof.

Each edge j represents a component block of the set comprising at least the following blocks:

- a reservoir in zero flow conditions similar to a

"generator" of imparted pressure and temperature;

the IPR of the reservoir i.e. the relationship between flow rate and pressure at the inlet of a well;

- the inlet of a conduit or of a well;

- the outlet of a conduit or of a well;

- the inlet of a valve;

- the outlet of a valve;

- the inlet of a separator;

- the outlet of a separator.

Each edge j, with the exception of those representing the IPR of the reservoir, connects one of the nodes I of the graph to an individual reference node (ref) which has no corresponding "physical" node in the system and to which are associated zero values of pressure, temperature and so on. In this way, each edge has a uniquely defined gradient ψ. This mathematical representation provides the algebraic instruments to write in analytical matrix form the relations of continuity of multiphase flow rate, of pressure, of temperature, etc. described below.

The simulation method 100 thus comprises a topological description step in which the topology of the system is described 140 by means of a set of topological equations derived from the aforementioned oriented graph and called topological system ∑ _T .

In particular, in the first place a first incidence matrix A=[aij] is generated that has as many rows as are the nodes of the graph and as many columns as are the edges of the graph. In detail, the element aij=l if there is an edge that connects the nodes i and j ; otherwise, the element aij=0.

From the first incidence matrix A is obtained a second incidence matrix B=[bij] and a loop matrix C=[Ckj] .

In detail, the second incidence matrix takes into account the orientation of the edges of the graph and therefore the element b±j=l if the flow ξj is flowing out of the i ^th node, the element bij=-l if the flow ξj is flowing into the i ^th , bi =0 if aij=0.

Once the second incidence matrix B is obtained, it is possible to write the equations Βξ=0 that impose the continuity of the mass flow rates of the liquid and gaseous phases at each node.

With regard to the loop matrix C, the element <¾=1 if the flow ξj is oriented like the k ^th loop, the element C _k =-1 if the flow ξ _j is not oriented like the k loop, C _k =0 if the vertex j is not part of the loop k.

With the loop matrix C it is possible to write the equations at the loops Ci|r=0, which impose the continuity of pressure, temperature and mass of the liquid and gaseous phases at each node of the system. The set of equations:

is indicated as topological system ∑ _T . The set of equations of the simplified models of valves, reservoirs and separators and of the non-differential equations of the simplified conduit models is indicated as static system ∑ _s . The set of differential equations of the simplified conduit models is indicated as thermo-fluid dynamic system ∑ _D .

The simulation method 100 according to the present invention, thus, comprises the step of solving 150 the complete system of equations comprising the plurality of topological equations and of the constitutive equations. The variables obtained by solving the aforesaid equations describe the thermo-fluid dynamic behavior of the production and transport system to be investigated .

Preferably, the step of solving 150 the complete system of equations can be represented schematically as shown in figure 8 and it comprises the use of an ODE (Ordinary Differential Equation) solver for solving the differential equations of the thermo-fluid dynamic system ∑ _D .

In detail, consider that x _n is the set of state variables of all the conduits at the generic instant t _n , u _n is the set of forcing variables and boundary conditions of the system at the same instant t _n , i.e.

- pressures and temperatures at the separators

- static pressures of the reservoirs

- temperatures of the reservoirs

- percentage of opening of the valves.

The ODE solver performs the step-by-step numerical integration of the time derivative of the state variables t _n +i , calculating the state variables at the next instant: where f _N ( ) is the set of all equations of the static system ∑ _s , of the thermo-fluid dynamic system ∑ _D and of the topological system ∑ _T and Λ is the set of corrective coefficients of all the simplified models.

The complete system f _N ( ) comprises "stiff" equations, i.e. whose numeric resolution with explicit integration methods is unstable. Therefore, a possible ODE solver that provides a good compromise between calculation rapidity and accuracy is the Trapezoidal Rule with the 2 ^nd order Backward Differentiation Formula, i.e. an implicit method that provides a first trapezoidal step and a second "2 ^nd order backward differentiation" step, methodologies known in themselves to the persons skilled in the art.

The step-by-step integration carried out by the ODE solver expressed by the equation (1) needs, at the first step, knowledge of the initial state xo . Preferably, the initial state xo is estimated imposing that said initial state xo be a system equilibrium state, i.e.

f_ _¥ (x _8. . u ₅ ;A) = o

The step of solving 150 the complete system of equations, with reference to figure 8, preferably comprises the following operations:

solving the system (∑ _T , ∑ _s ) i.e. the topological system ∑ _T and the static system ∑ _s considering as inputs the forcing variables u _n and the state variables x _n calculated at the n ^th step and obtaining the output variables y _n ;

- among the output variables, selecting the auxiliary variables a _n which are the following:

• pressure at the outlets of the conduits ρ ₀ υτ;

• liquid mass flow rate at the inlets of the conduits G _L ,IN;

• liquid mass flow rate at the inlets of the conduits GG, IN;

• temperatures at the inlets of the conduits ΤΊ _Ν ;

- solving the thermo-fluid dynamic system ∑ _D by means of the ODE solver considering as inputs the auxiliary variables a _n and the state variables x _n at the n ^th step. Figure 9 shows an exemplifying comparison of the results obtainable by the simulation method 100 according to the present invention and by a known reference finite volume thermo-fluid dynamic model illustrated by way of example in Figure 2 in the turndown operating condition (decrease of the flow rate in the line through partial closure of the well head valves) . It can be noted that the learning process corrects the discrepancies of the simulation method 100 bringing the error below 2% at steady state. The main differences between simulation method 100 and method used by the reference simulator consist of the reduction of the number of segments per individual conduit (from 100 to 1) and in of the reduction of the number of layers of insulation for the thermal description (from 21 to 1) .

The above description clearly illustrates the characteristics of the simulation method of the present invention, as well as its advantages.

In fact, the simulation method according to the present invention, using simplified analytical mathematical models to describe the various components of a hydrocarbons production and transport system, makes it possible to simulate in a simple and accurate manner the thermo-fluid dynamic behavior of the system.

The simplification of the present invention assures a significant performance increase in terms of calculation time, while maintaining adequate adherence to the physics of the system and sufficient reliability of the simulation.

The simplified models provided by the simulation method of the present invention comprise a set of corrective coefficients that are estimated through learning processes based on results obtained with finite volume thermo-fluid dynamic simulators or with actual measurements carried out on the field. In this way, the simulation method can be "trained" to describe the behavior of a specific geometry of the system in one or more operating regimes. Subsequently, it can be used to carry out simulations of the specific geometry of the system in other operating regimes or even partly modifying the geometry of the system, since the physics of the system is also represented in the simulation method. This technique can be extended to the simulation of the entire productive life of the system. This type of simulation requires markedly shorter times than finite volume thermo-fluid dynamic simulators and it makes it possible to identify very rapidly potentially critical situations for flow assurance (e.g. danger of formation of waxes, hydrates, etc.). The instrument of the patent also allows to select the situations in which it may be advisable to rely on a complete simulation for a more detailed and accurate analysis of the phenomena.

The combined use of finite volume thermo-fluid dynamic simulators and of a simulator based on the simulation method of the present invention makes it possible to reduce total simulation times and the computational burden for the analysis of a system, compared to the use of finite volume thermo-fluid dynamic simulators alone .

Lastly, it is clear that the simulation method thus conceived can be subject to numerous modifications and variations, without departing from the scope of the invention; moreover, all details can be replaced by technically equivalent elements.