Field of invention

This invention relates to an autonomous flow management system for regulating a multiphase flow in a pipeline-based transport system which utilises a novel computer-implemented method for predicting the multiphase fluid behaviour in the pipeline-based transport system.

Background

New oil and gas reserves are often discovered a considerable distance apart from existing oil and gas extraction facilities leading to a need for transporting unprocessed fluids over long distances, maybe hundred kilometres or more. Multiphase fluid transport over such distances are often challenging due to complex non-linear interactions between the fluid phases resulting in an irregular and unstable behaviour causing pressure drops, deposit formations, and/or liquid accumulations which may manifest themselves as a range of problematic flow phenomena such as slug flow, bubbly flow, stratified flow, annular flow, chum flow, etc.

The ability to correctly predict multiphase fluid flow in transportation systems, especially when long travelling distances are involved, is vital for the economy and operability for e.g. hydrocarbon production facilities, carbon capture and storage facilities etc. Accurate simulation models enabling predicting fluid flow behaviour consequences of different designs and approaches of the transportation system may save considerable investment costs.

Furthermore, due to the irregular and unstable nature of multiphase flows, it may and occasionally does arise flow instabilities which should be mitigated to avoid overflooding of slug-capturers, overloading of flow separators or other flow receiving equipment/flow handling facilities, and to maintain and stabilise the flow as near as possible to an intended plateau production level.

For example, liquid slugs are often problematic for pipelines by causing load variations and subsequent vibrations that reduce the lifetime of pipe fittings and other vulnerable components, and/or it may lead to need for inadvertent flaring. The ability to correctly simulate the slug flow characteristics of a given pipeline flow is therefore important for investigating the effect of and implementation of mitigating actions, such as e.g. topside choking, gas lift etc. Here, both simulation time (CPU) and accuracy of the predictions are essential.

Prior art

Prediction of the characteristics of multiphase flows requires specific computer- based simulations known as computational fluid dynamics (CFD). The computer codes of CFD-software usually consist of three main elements: (i) a pre-processor,

(ii) a solver, and (iii) a post-processor [Ref. 1, pp. 1 - 4],

The pre-processor element concerns the definition/input of the fluid flow problem to be simulated by the operator via an operator interface and the transformation of the definition/input to a computer readable form applicable for the solver. The preprocessor element typically comprises:

- input/definition of the geometry of the flow region of interest, often denoted as the computational domain,

- division of the computational domain into a set of non-overlapping subdomains, often denoted as grid generation,

- input/defining the physical and chemical phenomena to be simulated,

- input/definition of the fluid properties, and

- input/specification of boundary conditions at cells coinciding with the boundary of the computational domain.

The post-processor element concerns the output of the simulation/simulation results and may include data-base modules for storing data, data visualisation modules, graphic presentation modules etc.

The solver element concerns the numerical solution of the natural laws/equations governing transport phenomena, convection, diffusion and source terms (creation of destruction of an entity of the flow. The numerical algorithm for solving the equations typically consists of three steps: i) Integration of the governing equations of the fluid flow over the computational domain. ii) Discretisation - converting the integral equations into a system of algebraic equations, and iii) Solving the algebraic equations by iteration.

The governing equations of the fluid flow are mathematical statements expressing the conservation laws of physics to ensure conservation of mass, momentum and energy of the fluid. The fluid is treated as a continuum where its behaviour is described in terms of macroscopic properties such as velocity, pressure, density and temperature; see e.g. [Ref 1, pp. 9-10], The conservation laws of physics lead to a transport equation describing the motion of an arbitrary scalar fluid property, , which in a general formulation may be given on the form: (1) where “p“ is density, is the velocity vector, “D” is the diffusion coefficient, “ S<|>“ is a source term for flow property Φ|>, and “t” is the time coordinate. The first term (from left) describes the rate of increase of property (|) of the fluid element, the second term describes the net rate of flow of fluid property (|) out of the fluid element, the third term describes the net rate of increase of fluid property (|) of the fluid element due to diffusion, and the fourth term describes the rate of increase of quantity (|) due to sources.

The general formulated transport equation, eqn. (1), is usually a starting point for computational procedures in computational fluid dynamics by being applied to derive specific, often simplified, partial differential equations for the conservation of mass, momentum and energy for all fluid fields/phases of the fluid flow to be simulated in the specific geometry of the computational domain. CFD-simulations of gas-liquid flows in long pipelines usually requires using one-dimensional (ID) averaged conservation equations to achieve acceptable simulation times. For example, the mass conservation of a fluid phase, k, of a one-dimensional, unsteady and compressible multiphase flow with no source terms, eqn. (1) becomes reduced and may be given on the form: '

Here “ is the field mass of fluid phase k as a function of x, “w(x)k“ is the flow velocity for fluid phase k as a function of x, and p _k = a _k p _k where at is the volume fraction of fluid phase k and pk is the mass density of fluid phase k.

The computational domain specific derived conservation equations, such as e.g. eqn. (2), are integrated and discretised and solved over the computational domain by the solver element of the CFD-software. A widespread and much applied numerical solution technique is the finite volume method which is distinguished by treating the sub-domains of the computational domain as finite control volumes and integrating the governing equations of the fluid flow over each finite control volume (grid cell) of the computational domain. Applied on e.g. eqn. (2), the finite control volume integrates the eqn. over a finite control volume, AVi, and over a time step, A tn, and may be given on the form (3)

The next step in the finite volume method is transforming the continuous integral equation to a discrete algebraic equation. There are known several discretisation schemes for transforming the continuous integral equations into the discrete algebraic relations suitable for being solved numerically. The different difference schemes and their pros and cons in view of numerical stability and prediction accuracy are well known to the person skilled in the art. For flows dominated by convection, a well-known and much applied discretisation method for transforming the continuous integral equations into discrete algebraic relations is the upwind differencing scheme, see e.g. [Ref. 1, pp. 146-149], which applied on the mass conservation eqn. (2) may be given on the discrete form: (4)

Here, is the mass flux on the right internal face of the i’th finite control volume, Ai+1/2 is the cross-sectional area of the right internal face of the i’th finite control volume, Pk,i+i/2 is the field mass of fluid phase k at position xi+1/2, and uk,i+i/2 is the fluid flow velocity at position xi+1/2. Depending on the choice of grid (e.g. collocated or staggered), the values at position xi+1/2 may have to be interpolated using methods well known to the person skilled in the art, for example using a first-order upwind scheme. Correspondingly, is the mass flux on the left internal face of the i’th finite control volume. By “marching through” the entire computational domain, i.e. defining a discrete mass conservation equation for each of the entire set of i=l, ..., N finite control volumes, the continuity equation for the conservation of mass is transformed to a set of algebraic relations of discrete values which may be solved by an explicit (direct) or implicit (indirect, i.e. iterative) solution algorithm with given boundary conditions and initial flow parameters.

Implicit numerical solution schemes determine the solution for the state of the system at a given time-step forward in time by solving an equation involving both the current state of the system (which is known) and the state of the system forward in time (which is unknown). This requires use of a matrix or iterative solution which gives somewhat extra computation but has the advantage that the numerical solution scheme is unconditionally stable and thus enables employing long timesteps which saves computation. However, at the cost of less accurate predictions.

Explicit numerical solution schemes determine the solution for the state of the system at a given time-step forward in time directly from the current state of the system (which is known). This has the advantage that the solution is significantly more accurate but is encumbered by being conditionally stable. A necessary condition for the stability of an explicit numerical solution scheme is that the numerical domain of dependence bounds the physical domain of dependence, i.e. the “physical” velocity u(x) should be less than the “marching velocity” Ax/ At of the numerical method. This condition is known as the Courant-Friedrichs-Lewy (CFL) condition, which for a one-dimensional first order upwind scheme may be given as: Here At _n is the n’th time step increment in the calculation, Ax is the length of the i’th finite control volume and u(x _i+1 / ₂ , t _n +i) is the velocity of the fluid property at the upwind boundary of the i’th finite control volume at the end of the n’th time step. The CFL-condition should hold for every finite control volume of the computational domain and for every time step applied in the numerical calculations.

The limitation of the CFL-condition may be presented graphically such as e.g. given in figures 2a) and 2b). Both figures illustrate a section of a computational domain including five neighbouring finite control volumes having nodes at positions xi-3, xi- 2, ... xi+i, respectively. The mass inside each finite control volume at the beginning of a n’th time step is indicated graphically by the vertical height of the node points at the centre of each finite control volume, i.e. the total amount of mass present in the finite control volume at the beginning of the time step equals the area of a rectangle being defined by the west and east internal faces, the abscissa and a horizontal line (parallel with the abscissa) crossing through the node point such as illustrated in figure 2a) by the shaded area in the i-2’ th finite control volume. In the case of a first order upwind scheme, the mass passing through e.g. the internal face at position xi+1/2 during time step Atn is estimated as the area of the rectangle below the nodal point having a width determined as the product u(xi+i/2, tn+i) Atn, where u(xi+i/2, tn+i) is the x-component of the velocity vector at the end of the n’th time step. If this product is equal to or less than the spatial increment Axi (i.e. the width of the i’th finite control volume), the CFL -condition is fulfilled. The amount of mass estimated passing out of the i’th finite control volume (over the internal face at position xi+1/2) during the n’th time step is less than the total amount of mass being present in the finite control volume.

However, if the product u(xi+i/2, tn+i)- Atn is larger than the spatial increment Axi, a non-physical situation arises because the numerical calculations estimate that more mass will flow out of the finite control volume than it contains as illustrated in figure 2b), which leads to estimates including “negative” masses. Thus, a violation of the CFL-condition may cause the estimation of negative masses in the finite control volumes which may lead to numerical instabilities in the solver.

The CFL-condition may impose a rather harsh numerical load in CFD-calculations of multiphase flows involving rapidly moving phases and/or flow phenomena requiring a fine-masked grid to be predicted acceptably accurate. Then the CFL- condition may dictate use of very small time-increments leading to unacceptably extensive and heavy computational loads to make CFD-simulations of multiphase fluid flows a practically applicable tool for designing and/or trouble-shooting transport systems for multi-phase flows.

Griebel and Klitz [Ref. 4] discloses a fast and mass-conserving method for simulation of two-phase flow problems. One popular method is the coupling of level-set and volume-of-fluid (CLSVOF), which benefits from the advantages of both approaches and results in improved mass conservation while retaining the straightforward computation of the curvature and the surface normal. Despite its popularity, details on the involved complex computational algorithms are hard to find and if found, they are mostly fragmented and inaccurate. In contrast, this article can be used as a comprehensive guide for an implementation of CLSVOF into the existing level-set Navier-Stokes solvers on Cartesian grids in three dimensions.

US 2019/228124 discloses an improved computer implemented method for modelling transport processes in fluids is disclosed. Instead of modelling based on using an infinitesimal fluid element of a continuous medium, the method approximates fluid flow in a fluid system as a model gas flow in a model gas system being identical to the fluid system. The method is adapted to model gas flow including dilute gas flow for high Knudsen numbers (Kn). The method delivers a new basis for prediction of dynamic evolution of the model gas system by considering a pre-established or known dynamic history of the system during a preinitial period. A new generation of Computational Fluid Dynamics software products, which are based on the disclosed analytical tools and methods, are anticipated having capability to modelling gases from the continuum flow regime (Kn<0.01) to the free molecular flow regime (Kn>10), considerably higher accuracy of prediction, and lower computation cost.

From WO 2017/174532 it is known a simulation method (100) for simulating the thermo-fluid dynamic behaviour of multiphase fluids in a hydrocarbons production and transport system, said method comprising the following steps: outlining (110) the hydrocarbons production and transport system as a plurality of interconnected component blocks, thus creating a schematic representation; modelling (120) each component block with a simplified analytical mathematical model selected from the group of models comprising at least one conduit model, a valve model, a reservoir model and a separator model, each simplified analytical mathematical model comprising a plurality of constitutive equations adapted to describe the thermo-fluid dynamic behaviour of the corresponding component block; generating (130) an oriented graph on the basis of the schematic representation; determining (140) a plurality of topological equations on the basis of the oriented graph; determining (150) a plurality of output variables adapted to describe the thermo-fluid dynamic behaviour of the system by solving the set of the plurality of topological equations and of the constitutive equations.

WO 2019/164859 discloses a compositional simulator that is fully compositional and more robust and accurate is disclosed. Relative permeability (kr) and capillary pressure (Pc) are modelled as state functions, making them unique for a given set of inputs, which can include Euler characteristic, wettability, pore connectivity, saturation, and capillary number. All of these are made to be a function of composition, T, and P or rock properties. These state function kr-Pc models are fully compositional and can fit experimental data, including complex processes such as hysteresis. The models can be tuned to measured relative permeability data, and then give consistent predictions away from that measured data set. Phase labelling problems are eliminated. Flux calculations from one grid block to another are based on phase compositions. Simulations for three or four-phase hydrocarbon phases are possible. Time-step sizes increase to stability limits of implicit-explicit methods. Fully implicit methods are possible and significant improvements are expected.

From US 2018/260499 it is known a fluid is modelled as a set of discrete particles. Each of the particles is associated with one or more properties, and a spatial distance comprising a smoothing length over which the one or more properties are to be smoothed. A corresponding trajectory is simulated for each of the particles. The corresponding trajectory is used to formulate a first solution for simulating transport within the fluid. A first predicted error is determined for the first solution. An iterative adjustment is performed to at least one of a quantity of particles, the smoothing length, or the one or more corresponding properties, to formulate a second solution for simulating transport with the fluid, and a second predicted error is determined for the second solution, until the second predicted error is within a predefined boundary.

Objective of the invention

The main objective of the invention is to provide an autonomous flow management system for regulating a multiphase flow in a pipeline-based transport system which utilises a computer-implemented method for predicting fluid behaviour of a multiphase flow in the pipeline-based transport system having a solver enabling using an explicit numerical solution scheme being stable independently of the Courant-Friedrichs-Lewy (CFL) condition.

A further objective of the invention is to provide an autonomous flow management system for regulating a multiphase flow in a pipeline-based transport system which utilises a computer-implemented method for predicting fluid behaviour of a multiphase flow in the pipeline-based transport system having a solver enabling using an explicit numerical solution scheme being stable independently of the Courant-Friedrichs-Lewy (CFL) condition, and which preserves the positivity of the mass everywhere.

Another objective of the invention is to provide an autonomous flow management system utilising a computer-implemented method for predicting fluid behaviour of a multiphase flow in the pipeline-based transport system to trouble-shoot and regulate the multiphase flow to mitigate flow instabilities. Description of the invention

The autonomous flow management system of the present invention is based on utilising a specific computer-implemented method believed to be novel for simulating a multiphase flow flowing in a pipeline-based transport system to regulate the flow with the aim to stabilize the multiphase flow at an intended flow volume plateau and/or to mitigate eventual flow instabilities/irregularities which may arise.

The autonomous flow management system according to the invention is not tied to any specific configuration of the flow management system but may have any configuration known to the skilled person as long as it comprises one or more sensor(s) monitoring one or more characteristic flow parameter(s) of a (physical) multiphase flow flowing in a (physical) pipeline-based transport system being managed by the autonomous flow management system, computer-based means for executing the specific computer-implemented method applied on a virtual model of the (physical) pipeline-based transport system with the monitored characteristic flow parameter(s) as boundary condition(s), and one or more flow regulating actuator(s) located in the (physical) pipeline-based transport system and adapted to regulate the multiphase flow, and a production management unit applying the simulation results to adjust/regulate the (physical) flow of multiphase fluid in the (physical) pipeline-based transport system by activating the one or more (physical) flow regulating actuators. The monitoring/registering of the one or more characteristic flow parameters should preferably be applied at an upstream end and at a downstream end of the (physical) pipeline-based transport system.

The specific computer-implemented method according to the invention for simulating the multiphase flow utilises the possibility of preserving the positivity of the mass in the numerical calculations independently of the Courant-Friedrichs- Lewy (CFL) condition by applying an approach for estimating the mass flux which comprises for a n’th time step, Atn; applying a polynomial to spatially reconstruct the mass, _k (x), for the k’th fluid phase being present in each of the N finite control volumes of the computational domain, and then applying the discretised flow velocity, uk (which can be located either at the centre point of the mass cells in a collocated grid, or at the internal faces in a staggered grid), to determine a domain of dependence, D _fcJ -, representing the distance the k’th fluid phase has travelled from its starting point to a j'th internal face, where j e 1/2, 3/2,..., N-l/2, during the n’th time step for each internal face j of the computational domain, and sum the spatially reconstructed mass being present in the domain of dependence, D _fcJ - and assume it passes through the j’th internal face during the n’th time step. The spatial reconstruction of the mass of the k’th fluid phase over the computational domain may be obtained by any manner known and conceivable to the skilled person, provided the following conditions are satisfied: ^ The reconstructed mass function should have, within each finite control volume, an average value equal to the discretised mass at the cell's centre point. In other words, it should conserve the mass present in each cell In one example embodiment, the spatial reconstruction of the mass of the k’th fluid phase over the computational domain may also satisfy the condition that there must not be any negative value anywhere in the reconstruction. Figure 3a) illustrates schematically an example of using a first order polynomial to represent the mass present in the finite control volumes as linear interpolation lines over each finite volume. Use of higher order polynomials gives a more accurate spatial reconstruction. The range and length of the domain of dependence, depends on the flow direction, the flow velocity, and the duration of the n’th time step. If the flow direction is positive, the domain of dependence stretches from the j’th internal face in a westward direction as shown in figures 3b) and 3c), first in cell i (the j’th internal face is between control volumes i and i+1). As seen on the figures, the domain of dependence, ^^ _^,^ , may be longer than the i’th finite control volume such that it extends into one or more neighbouring finite control volumes. In the example shown schematically in figures 3a) and 3b), the domain of dependence stretches into two and a half neighbouring finite control volumes from position being inside the i-3’th finite control volume to position x j. If the flow direction is negative, the domain of dependence stretches in eastward direction (not shown on the figures). Figure 3c) illustrates the mass being present in the domain of dependence, ^^ _^,^ , by the shaded area between the abscissa and the graphic representation of (in this example, linear interpolation lines). The shaded area over the domain of dependence represents the sum of mass which is assumed passing through the j’th internal face during the n’th time step. Since all mass present in the domain of dependence is summed and assumed passing through the corresponding internal face, the above approach ensures that the mass passing through the internal face actually exists, thus not pulling more mass out of a control volume than is present. Furthermore, since the domain of dependence may extend over more than one finite control volume, the numerical method is stable also when the “physical” velocity, u(x) becomes larger than the “marching velocity”, Δx/Δt, of the numerical scheme. By the above solution, it is obtained that the 1D-simulation of the multiphase flow becomes more accurate/realistic at similar or smaller computer loads as compared to prior art 1D-multiphase simulations. This enables the flow management system according to the invention to obtain a more stable mode of operation, which in turn can give an improved plateau production level increasing the flow volumes of multiphase fluid which may be transported through the transport system and/or improved mitigation of flow irregularities which may arise.

Thus, in a first aspect, the invention relates to an autonomous flow management system 100 for regulating flow-through of a multiphase fluid containing a number of k, where k is a positive integer, fluid phases, through a pipeline-based transport system 10, comprising:

- a flow simulation unit 30,

- a sensor configuration comprising at least a first sensor 51 located at an upstream end 11 and a second sensor 51 located at an downstream end 12 of the pipeline-based transport system 10 and measuring one or more characteristic flow parameter(s) of the multiphase fluid flowing through the pipeline-based transport system 10,

- an actuator configuration comprising at least one actuator 61 adapted to regulate the flow of fluid through the pipeline-based transport system 10, and

- a control unit 20 adapted to:

- receive signals 53, 54 from the sensor configuration measuring one or more characteristic flow parameter(s) and transferring the signals to one or more boundary conditions 21 passed on to the flow simulation unit 30, and

- receive simulation results 32 from the flow simulation unit 30 and transferring the simulation results to set point values 22 passed on to the actuator(s) 61 of the actuator configuration , wherein

- the flow simulation unit 30 comprises a computer loaded with a software, which when executed performs a computer-implemented method simulating the fluid behaviour of the multiphase flow flowing in the pipeline-based transport system 10 with the boundary condition(s) 21 from the control unit 20, characterised in that the computer-implemented method comprises: applying a one-dimensional (ID) computational fluid dynamic (CFD) model describing the geometry of the pipeline-based transport system and the multiphase flow flowing therein, solving the ID CFD model to simulate the fluid behaviour of the multiphase flow flowing in the pipeline-based transport system, and describing the determined fluid behaviour as macroscopic fluid properties such as flow velocity, pressure, density, temperature, etc., and the ID CFD model applies a finite volume method to solve the model, wherein the geometry of the pipeline-based transport system is defined as a computational domain extending along an axis represented by the cartesian coordinate x and being divided into a set of N, where N is a positive integer, non-overlapping finite control volumes separated by an internal face between adjacent finite control volumes, and wherein the one-dimensional computational fluid dynamic model is adapted to estimate the mass flow of a k’th fluid phase out of a i’th finite control volume during a n’th time step by applying a polynomial to spatially recon- struct the mass, ^^^ _^ ^{^} ^^ ^{^} , of the k’th fluid phase being present in each of the N finite control volumes of the computational domain, and then for each j’th internal face, where j є 1/2, …, N+1/2, of the computational domain, execute the following steps i) and ii): i) reconstruct the flow velocity, uk(x), of the k’th fluid phase as a function of position x and apply the reconstruction to determine a domain of dependence, ^^ _^,^ , representing the distance the k’th fluid phase has travelled during the n’th time step, and ii) sum the spatially reconstructed mass being present in the domain of dependence, ^^ _^,^ , and assume the summarised mass passes through the j’th internal face over the n’th time step, into the i’th finite control volume when u _k (x _j ) < 0 or into the i+1’th finite control volume when uk(xj) > 0, where uk(xj) is the flow velocity at the j’th internal face. As used herein, the term “computational domain” is a virtual representation of the geometry of a section of a pipe containing a multiphase flow which is to be simulated. The computational domain is divided into a number of N non-overlap- ping finite control volumes where each i’th, i = 1,…,N, finite control volume represents a separate “slice” of the computational domain. As used herein, index “i” is applied when the focus is on the finite control volumes (i є 1, …, N) and index “j” is applied when the focus is on the internal faces j є 1/2, …, N+1/2. The notations “i-1/2” and “i+1/2” as used herein relate to the internal faces of the finite control volume i (left and right internal faces, respectively), and “j-1/2” and “j+1/2” are used to relate finite control volumes to the internal face j (left and right control volumes, respectively). The terms “left” and “west” are used herein interchangeably to indicate a direction opposite to the grid marching direction (i.e. towards decreasing index “i” or “j”), and consequently, the terms “right” and “east” are used herein interchangeably to indicate a direction in the grid marching direction (i.e. towards increasing index “i” or “j”). As used herein, the subscript n denoting the n’th time step is omitted to simplify the notation since all relations relate to the n’th time step. The concepts of the finite volume method, explicit or implicit numerical solution algorithms, upwind differencing scheme, staggered or collocated grids, etc. and how to implement these in CFD-models are well known and mastered by persons skilled in the art.

The present invention is not limited to any choice of discretisation scheme or numerical solution algorithm except that the CFD-model shall apply a finite volume method approach. That is, the numerical scheme according to the first aspect of the invention may be applied in any ID CFD-model regardless of which differencing scheme and/or explicit or implicit numerical solution algorithm being implemented in the CFD-model. However, the numerical scheme according to the invention is especially beneficial for explicit numerical schemes by enabling violating the CFL- condition without excessively compromising the numerical stability, making CFD- models having an explicit numerical scheme a preferred example embodiment. Explicit formulated numerical schemes have a relatively high numerical accuracy and are well suited for multicore CPUs and computational cluster machines frequently used in cloud computing.

The numerical scheme according to the first aspect of the invention relates to the solution of the transport equations for the mass transfer of each fluid phase and field present in the multiphase flow. The transport equations governing momentum and energy of the fluid phases and fields of the multiphase flow may be solved by the CFD-model in any known manner. Thus, the present invention encompasses any CFD-model having a solver for solving the transport equations where the transport equations governing the mass conservation/flow of mass is solved according to the first aspect of the invention, and where the transport equations governing the momentum and energy are solved in any known manner. The transport equations governing momentum and energy may be subject to a CFL-condition.

The reconstruction of the flow velocity, uk(x), of the k’th fluid phase as a function of position x to determine a domain of dependence, D _fcJ -, may be obtained in several ways known to the skilled person. The invention according to the first aspect may apply any conceivable manner known to the skilled person.

As an example, linear interpolation of the discretised flow velocity within each finite control volume may be applied to reconstruct the flow velocity, uk(x), as a function of the position x, here given for the i’th finite control volume and a staggered mesh arrangement: (6)

Eqn. (6) may be rearranged to the form: (7)

Equations (6) and (7) can be written for all the finite control volumes. For example, if they are written for the i+l’th finite control volume, they will be valid in x ∈ - By joining all the finite control volumes, a reconstruction for any x in the computational domain is obtained.

The domain of dependence, for the j’th internal face at the n’th time step may be determined by applying and integrating the equation of the flow characteristics: (8) with the velocity as expressed in eqn. (7). This gives the following relation:

(9) where is the starting position of the flow characteristic of the k’th fluid phase at the beginning of the applied time step, crossing the position at a time t. The finite control volume index i denotes the cell in the upwind direction and is equal to j - 1/2 when the flow direction is positive, i.e. when uk(xj) > 0. When the flow direction is negative, i.e. when uk(xj) < 0, eqn. (11) still applies if the index i is set to be j + 1/2. Eqn. (9) may be solved for , under the requirement that x _k (At) = x _; -, which means that at the end of the time step At, the characteristic will cross the internal face j . This gives:

(10)

If the starting position lies within the i’th finite volume, the domain of depend- ence, D _kJ -, for the k’th fluid phase and the j ’th internal face at the n’th time step, will extend from the starting position to the j'th internal face. In this case it may be advantageous to denote to indicate that it is the starting point of the domain of dependence,

( H )

One advantage of the invention according to the first aspect is that it enables applying time steps being larger than the time the flow characteristic uses to cross the i’th finite control volume without compromising the numerical stability. However, if the n’th time step is larger than the time the flow characteristic uses to cross the i’th finite control volume, the starting position determined by eqn. (10) will be lying outside of the i’th finite control volume. To accommodate for this situation, the invention according to the first aspect may further comprise a procedure to determine the domain of dependence which takes into account the time needed for the k’th flow characteristic to cross the finite control volumes. The time needed for the flow characteristic to cross a finite control volume may be determined by solving eqn. (10) for At, which when applied in the i’th finite control volume may be given as:

(12)

Here At _{c k} is the time needed for the k’th flow characteristic to cross the i’th finite control volume. Eqn. (12) may also be applied to determine the time needed for the k’th flow characteristic to cross the i-l ’th finite control volume by applying the variables related to the i-l’th finite control volume: ^an d Eqn. (12) may be applied similarly to determine the crossing time for any other finite control volume. If the ratio of velocities is negative, Eqn. (12) becomes invalid. It happens when the velocity changes sign in the cell. In this case, no characteristic can cross the position where u _k (x) = 0. Thus, the crossing time is infinite, and the characteristic will start within the cell.

Close to the boundary nodes, the domain of dependence may extend out of the computational domain P. This means that the total crossing time from the computational domains' west or east end to the internal face j is smaller than the time step At. Then, there is a fraction of the time step for which the mass crossing internal face j will originate from outside the computational domain. To accomplish for this situation, the method may be extended with a special treatment of the nodes. Depending on the type of boundary condition that we want to achieve, two quantities may be useful to calculate, the remainder of the time step ^t _{r k} , and the extension of the domain of dependence outside of the computational domain.

When the domain of dependence would stretch outside of the computational domain, P, the remainder of the time step At _{r k} is defined as

(12a)

For the extension of the domain outside of the computational domain, P, the velocity outside of the pipe must be extrapolated. As an example, it can be constantly extrapolated and equal to the velocity in the west- or eastmost cell, for the west or east node, respectively. Then, equation (10) can be applied with the remainder of the time step ^t _{r k} , and u _{k i} set equal to the extrapolated velocity. Equation (10) in this case will have a numerical issue, since there is at the denominator due to the constant extrapolation. This issue is solved f urther down. The notation for intervals as used herein follows the international standard ISO 80000-2, where the brackets “[“ and “]” indicate a closed interval border, and the parenthesises “(“ and “)”indicate an open interval border. For example, [a, b] is the closed interval containing every real number from a included to b included: [a, b] = {x ∈ ℝ ∣ a ≤ x ≤ b}, while (a, b] is the left half-open interval from a excluded to b included: (a,b] = {x ∈ ℝ ∣ a < x ≤ b}. An example embodiment of a procedure for determining the domain of dependence of the j’th internal face and k'th fluid phase when ^^ _^,^ ് 0, may be given as: Initialisation: i) if u _k,j > 0: set the upwind cell index i = j-1/2 and the direction index s = -1 if uk,j < 0: set the upwind cell index i = j+1/2 and the direction index s = 1 (ii) set the cell counter m = 0 (iii) set Δt _r , the remainder of the time step after crossing cells ^^ to ^^ ^ ^^ ^^, equal to the full time step Δt Step 1: (iv) apply eqn. (12) to determine the crossing time of the (i+sm)’th finite control volume, Δtc,k,i+sm (v) if Δt _c,k,i+sm > Δt _r (the characteristic will not cross the whole (i+sm)’th finite control volume: Go to step 2 or else (the characteristic will cross the whole (i+sm)’th finite control volume and enter next one): if i+s(m+1)<1 or i+s(m+1)>N: (the end of the pipe's domain has been reached) set Δt _r = Δt _r – Δt _c,k,i+sm , and terminate the procedure, or else: set m = m + 1 and Δt _r = Δt _r – Δt _c,k,i+sm , and go to step 1 Step 2: (vi) Apply eqn. (10) with the flow velocity shorthand’s and as defined below eqn. (7), and timestep Δt=Δt _r to determine the characteristic starting position, which defines the domain of dependence ^^ _^,^ according to eqn. (11) and term inate the procedure. The above procedure for determining the domain of dependence, may be presented graphically as shown in e.g. figures 3b) and 3c), which illustrates a domain of dependence stretching from the j ’th internal face across three neighbouring finite control volumes and a distance into a fourth neighbouring finite control volume. The flow direction is indicated by the black arrows located at the internal faces. The figure illustrates an example with a positive flow direction such that the domain of dependence stretches from the internal face at position xj to the starting position lying inside the i-3 ’th finite control volume.

The above procedure for determining the domain of dependence is described for either a negative or a positive flow direction. In the case of zero flow velocity there is no flow to be modelled and the domain of dependence will be zero. This situation may be taken care of by setting the domain of dependence, to zero when uk,j = 0. The particular case where the velocity change sign from one internal face to the next is naturally handled, since the condition to leave step 1, where will be true. Step 2 has no particular issue with velocity changing sign.

The above procedure for determination of the domain of dependence may have an issue when the velocities in cells i-1/2 and i+1/2 are equal to each other, since the denominator u _{k i} of eqns. (10) and (12) goes to towards zero. However, taking the limit of constant velocity we have that: (13) which shows that the Au _{k £} in the denominator is a numerical issue rather than a fundamental problem with the above approach. The above approach may thus be made numerically robust by e.g. defining the change in the Courant number over the i’th finite control volume, ACi, to be: (14) and apply the expression of eqn. (14) to rewrite the problematic term in eqn. (10) as: < ¹⁵ ) and expand the expression in a Taylor series: (16)

If the second order term in eqn. (16) is less than the e value returned from the Fortran function EPSILON which returns a positive real value being the smallest possible value of e making 1+e not equal to 1 on the computer system, then the numerical value of the function on the left side of eqn. (16) is numerically indistinguishable from the two first terms of the Taylor expansion, such that a numerically robust form of eqn. (10) may be given as: (17) where the index i=j -1/2 when the flow direction is positive, i.e. when

Uk(xj) > 0. When the flow direction is negative, i.e. when uk(xj) < 0, eqn. (17) still applies if the index i is set to be: i=j+ 1/2. At is the applied time step.

Similarly, the time needed for the flow characteristic to cross a finite control volume may be made numerically robust by defining r as the ratio:

(18) and rewrite eqn. (12) on the form:

(19)

The Taylor expansion of the right-hand side of eqn. (19) around rk,i = 1 is:

If the second order term of eqn. (20) is less than the e value returned from the Fortran function EPSILON, then the numerical value of the function on the left of eqn. (20) is indistinguishable from 1, such that a numerically robust expression for the flow characteristic to cross a finite control volume, may be given as: ( ₂₁₎

In an example embodiment, the invention according to the first aspect may further comprise determining the domain of dependence, D by executing the following steps i) to vi):

Step 1 (initialisation): i) if uk,j > 0: set an upwind cell index i = j-1/2 and a direction index s = -1, if uk,j < 0: set an upwind cell index i = j + 1/2 and a direction index s = 1, ii) set a cell counter m = 0, 111) set Atr, the remainder of the time step after crossing cell i to i + sm, equal to the full time step At,

Step 2: iv) otherwise: where 1 and e is a positive real number obtained from executing a Fortran

EPSILON computer implemented function, if Ate, k, i+sm > Atr: go to step vi), or else if i+s(m+l)<l or i+s(m+l)>N: set Atr = Atr - Ate, k, i+sm, and terminate the procedure, or else: set m = m + 1 and Atr = Atr - Ate, k, i+sm, and go to step iv),

Step 3): vi) determine a characteristic starting position, by solving the following eqn.; and go to step vii),

Otherwise: determine a characteristic starting position, by solving the following eqn.; and go to step vii), where

Step 4:

This terminates the procedure.

As mentioned above, the spatial reconstruction of the mass in the finite control volumes may be obtained in several ways known to the person skilled in the art. The invention according to the first aspect may apply any mathematical technique for spatially reconstructing the mass in the finite control volumes known and conceivable to the skilled person, provided the following condition is satisfied:

• The reconstructed mass function should have, within each finite control volume, an average value equal to the discretised mass at the cell's centre point. In other words, it should conserve the mass present in each cell

In one example embodiment, the spatial reconstruction of the mass of the k’th fluid phase over the computational domain may also satisfy the condition that there must not be any negative value anywhere in the reconstruction.

The probably simplest spatial reconstruction is a piecewise constant function, equal within each finite control volume to the discretised mass at the finite control volume centre point. Another example is using a piecewise linear function composed of first order polynomial, one per finite control volume, whose average over the volume is equal to the discretised mass at the finite control volume centre point. Higher order polynomials can also be used and will result in a higher convergence order of the scheme. Especially even-order polynomials are suited since they allow expanding the stencil equally to both sides of the finite control volume being subject for mass reconstruction. Thus, the spatial reconstruction of the mass as a function of position variable x may be obtained by applying an even numbered polynomial: (22) and then, for each i’th finite control volume, where i e 1, ..., N, of the computational domain, and each k'th phase, define a set of coefficients, c _a , through the condition: ( ²³ ) and solve for the coefficients co, ci, ..., cp to express the spatial reconstruction of the mass of k'th phase present in the i’th finite control volume as p There is one set of coefficients co, ci, ... , cp for each phase k in each cell i, but the k and i indices are dropped for readability. By repeating this procedure for every finite control volume of the computational domain, the mass of each finite control volume is spatially reconstructed over the entire computational domain.

The simplest spatial reconstruction is to set P = 0 such that p(x) becomes the average mass of phase k in the i’th finite control volume at the beginning of the n’th time step

An example embodiment employing a second order polynomial; p(x) = c ₀ + c-px + c ₂ x ² is as follows. The coefficients c ₀ , q, c ₂ for phase k in the i'th cell are calculated by requiring that the integrals of the polynomial in the three ranges x e [xi-3/2, xi-1/2], x e [xi-1/2, xi+1/2] and x e [xi+1/2, xi+3/2], be equal to the masses in the i-l’th, the i’th, and the i+l’th finite control volumes, respectively. Using that: (24) the system to solve becomes: (25) (26) (27)

Solving eqns. (25) to (27) determines the coefficients co to C2 which may be used to reconstruct the spatial reconstruction of the mass present in the i’th finite control volume as: Combining the polynomial pieces for all the control volumes in the domain give the spatial reconstruction of the mass over the whole domain,

Combining the spatially reconstructed mass for the whole domain with the domains of dependence at all internal faces, the mass fluxes between the finite control volumes in Eqn. (4) may be calculated in a new manner. The mass flux over the j'th internal face is calculated by summing the mass over its associated domain of dependence by evaluating the integral < ²⁸ ) where is the spatial reconstruction of the mass of the k’th fluid phase over the whole computational domain, composed of all the pieces for i G [1, IV], which is integrated over the domain of dependence, , and k,j is the mass flux across the j’th internal face. Aj is the cross-sectional area at the position of the j'th internal face. Higher order polynomials may have regions where they have negative values leading to summation of negative masses. Thus, in an example embodiment, the invention according to the first aspect may further comprise a rescaling of the spatial reconstruction of the mass which may by a limiter function known from [Ref. 2, slide 21] which limits the polynomial to zero or positive values while maintaining the condition defined in eqn. (23): where if m < 0 or otherwise 0 = 1, and where

Close to the pipe nodes, the domain of dependence may reach outside of the computational domain P. Then, part of the mass flowing through will originate from the outside of the computational domain. Depending on the type of boundary condition that we want to achieve, the following can be done:

In one example embodiment applying an imposed mass flow rate boundary condition, the mass of phase k flowing through internal face j during the remainder of the time step, Atr,k, that was stored in the algorithm to determine the domain of dependence, will be: where F _{in k} is the imposed mass flow rate. Then, the mass flow rate through the internal face during time step At will be determined by integral (28) to which is summed the flow rate during the remainder of the time step, as (28a) where denotes the part of the domain of dependence which is within the computational domain.

In another example embodiment applying an imposed pressure boundary condition, the velocity may be applied to decide the mass flow rate, given a user-defined phase split of mass entering the pipe. In turn, the difference between the user-imposed pressure and the pipe's pressure decides the velocity update during the time step, through the pressure gradient term. This may be handled by constantly extrapolating the velocity and applying Equation (17) with AC _{k £} = 0, and being the velocity in the first or last cell, for the west or east boundary condition, respectively. Equation (17) gives in this case a outside of the computational domain. Then, the mass flow rate through the internal face during time step At may be determined by integral (28), in which for the part of the domain of dependence outside the computational domain it is applied a mass defined as where the index in denotes the volume fraction imposed at the node and the density corresponding to the imposed pressure.

Numerical schemes which are higher order in space are known to generate spurious oscillations at discontinuities and extrema. On the other hand, numerical schemes which are first order in space do not. Thus, in one embodiment, the rescaling of the spatial reconstruction of the mass further comprise a step to modulate the order of the scheme locally using limiters, based, in each cell, on the value of the solution in the cell and its immediate neighbours. This gives a combined spatially first- and higher-order numerical scheme.

The spatial reconstruction of the mass may be made locally zeroth order by setting 0 = 0 in eqn. (29), that is, make the reconstruction constant in a cell. A zeroth-order mass reconstruction gives a spatially first-order numerical scheme. Thus, the limiter is defined to prevent the reconstructed polynomial from overshooting or undershooting the neighbouring cell values, by using a limiting condition of the same form as the one limiting negative values, i.e. applying eqn. (29) where, 0 is determined as follows:

At the left internal face of the i’th internal control volume, set: and at the right internal face, set: and then set 6 = MIN(0 _L , 0 _R f and in addition, at extrema, that is, for cells satisfying the following condition

Finally, we keep the most restrictive 0 (closest to 0) between the presently defined limiter and the limiter for the positivity preserving property in eqn. (29).

As used herein, the term “pipeline-based transport system” encompasses all components of the transport system necessary to transport the fluid including pipeline segments, splits, joins, valves, pumps, sources, sinks, etc. An example of a pipeline-based transport system for produced liquids in oil and gas extraction is shown in figure 1 which illustrates schematically an example embodiment of such transportation system. This example embodiment comprises a plurality of tie- backs/pipelines (2) connecting a production well (1) to a nearby satellite hub (3) which collects the produced fluid in a region and passes the produced fluid in a satellite pipeline (4) to a common hub (5). The example embodiment comprises four satellite hubs (3) connected to the common (5) by a satellite pipeline (4) each. The common hub (5) passes the produced fluid to a processing facility located either offshore on the sea surface via a riser (not shown in this embodiment) or to an onshore production facility via fluid transportation pipeline (6). The transport system usually involves one or more fluid pumps (7) to provide the necessary flow pressure to move the fluids through the transport system. The above example embodiment should not be interpreted narrowly. The pipeline-based transport system may have any conceivable configuration from a single pipeline for fluid transport, to interconnected networks pipelines for fluid transport in e.g. chemical process industry plants, for connecting offshore production facilities to onshore produced fluid receiving facilities etc.

The configuration of an example embodiment of the flow management system according to the invention is schematically illustrated in the diagram shown in figure 11. The pipeline-based transport system being managed by the flow management system 100 is shown schematically on the figure as a box 10 having an upstream end 11 receiving a fluid to be transported through the transport system to a downstream end 12 where the fluid is delivered to a fluid receiving facility. The flow management system comprises further control unit 20, a flow simulation unit 30, a sensor configuration comprising at least a first sensor 51 located at the upstream end 11 and a second sensor 52 located at the downstream end 12 of the pipeline-based transport system 10, and an actuator configuration comprising at least one actuator 61 adapted to regulate the through flow of one or more fluid phases of the multiphase flow.

In another example embodiment of the flow management system according to the invention, the flow management system further comprises a second actuator 62 located at the downstream end 12 and/or a third actuator (not shown in the figure) located anywhere in-between the upstream 11 and downstream 12 end of the pipeline-based transport system 10. In this embodiment, the set point-value for the second actuator 62 is transferred from the control unit as signal 23 and the set pointvalue for the third actuator is transferred from the control unit as signal 24.

In one embodiment, the control unit 20 may be a Distributed Control System, a Programmable Logic Controller, an Edge Gateway, a SCADA system or a Historian System or Timeseries Database being implemented to covering automation layers 0, 1, and 2 according the standard: ANSI/ISA-95.00.01-2010 (IEC 62264-1 Mod) Enterprise-Control System Integration - Part 1 : Models and Terminology.

The computer-implemented method for predicting the fluid behaviour need information of the gas and liquid phase ratios and the temperature of the multiphase flow entering the transport- system at its upstream end to predict the fluid behaviour. In some appliances, the flow rates and thus gas and liquid phase ratios of the flow entering the pipeline-based transport system is constant or practically constant. In such cases the information of the gas and liquid phase ratios may be entered as an input variable for the computer-implemented method. Thus, in one embodiment., the sensor configuration of the flow management system comprises at least a temperature sensor located at the upstream end of the pipeline-based transport system. In other appliances, the gas and liquid phase ratios may vary. This, in one embodiment, the sensor configuration of the flow management system comprises at least a flow sensor and a temperature sensor, both located at the upstream end of the pipeline-based transport system.

The control unit 20 receives sensor signals 53, 54 from the sensor configuration , which typically is electric signals, and transforms them into one or more measured flow parameter(s) such as e.g. flow velocity of one or more fluid phases, pressure, temperature, density of one or more fluid phases, volume or mass fraction of one or more fluid phases, etc. These measured one or more flow parameter(s) are passed on to the flow simulation unit (30) and applied as boundary condition(s) in the simulation of the multiphase flow.

The computer-implemented method for predicting the fluid behaviour need information of the gas and liquid phase ratios and the temperature of the multiphase flow entering the transport- system at its upstream end to predict the fluid behaviour. In some appliances, the flow rates and thus gas and liquid phase ratios of the flow entering the pipeline-based transport system is constant or practically constant. In such cases the information of the gas and liquid phase ratios may be entered as an input variable for the computer-implemented method. Thus, in one embodiment, the first sensor 51 of the sensor configuration of the flow management system comprises at least a temperature sensor located at the upstream end 11 of the pipeline-based transport system 10. In other appliances, the gas and liquid phase ratios may vary. Thus, in one embodiment, the first sensor 51 of the sensor configuration of the flow management system comprises at least a flow sensor and a temperature sensor, both located at the upstream end 11 of the pipeline-based transport system 10.

In one embodiment, the first sensor 51 of the sensor configuration of the flow management system comprises a temperature sensor located at the upstream end 11 of the pipeline-based transport system 10, and: either:

- the first sensor 51 further comprises a pressure sensor and the second sensor 52 comprises a pressure sensor,

- the first sensor 51 further comprises a pressure sensor and the second sensor 52 comprises a flow sensor,

- the first sensor 51 further comprises a flow sensor and the second sensor 52 comprises a pressure sensor, or

- the first sensor 51 further comprises a flow sensor and the second sensor 52 comprises a flow sensor.

The simulation results from the flow simulation unit are applied to regulate the flow in the pipeline-based transport system 10 by adjusting the actuator(s) 61 of the actuator configuration to set-point values determined by the control unit 20 taking the flow simulations results into account. The set point values may be determined by using one or several of the following algorithms that should be well known to those proficient in the art: PID control loop, Pre-trained machine learning algorithm, and/or Global or local optimum search algorithm.

In one embodiment, the actuator 61 of the actuator configuration may be either a control valse, a drum separator, a compressor, a gas injector, a pump etc., and the set-points, depending on which actuator(s) being applied, may be one or more of a setpoint for speed or power of rotating or reciprocating process equipment, as for example a compressor or a pump, a set point for a valve, a binary start/stop or open/close signal, etc.

In a third aspect, the invention relates to a method for trouble-shooting flow problems during operation of a pipeline-based fluid transportation system for transporting a multiphase fluid flow, wherein the method comprises:

- applying the computer-implemented method of the flow management system according to first aspect of the invention to predict the effect on the fluid behaviour in the transport system from a set of possible mitigation actions, determine which mitigation action which is to be physically implemented on the transport system having flow problems, and implementing the mitigation action by engaging the actuator configuration .

In a fourth aspect, the invention relates to a method for avoiding flow problems during operation of a pipeline-based fluid transportation system for transporting a multiphase fluid flow, wherein the method comprises applying the computer- implemented method of the flow management system according to anyone of claims 1 to 16 to predict the fluid behaviour in the transport system for a timespan into the future and determine the need for initiating mitigation actions during the timespan to define an operational schedule to be implemented by the actuator(s) of the actuator configuration.

The mitigation actions may be regulating the flow volumes in the transport system, topside choking, artificial lift (e.g. gas lift), and others.

List of figures

Figure 1 illustrates schematically an example embodiment of a pipeline system for transporting processed fluids in oil and gas extraction. Figures 2a) and 2b) illustrate graphically the limitation of the CFL-condition when estimating mass fluxes out of control volumes.

Figure 3a) is a schematic representation of an example of a spatial reconstruction of the mass of the k’th fluid phase in the finite control volumes using a first order polynomial to represent the mass as linear interpolation lines over each finite volume.

Figure 3b) is a schematic representation of the same spatial reconstruction as shown in figure 3a) and which also illustrates an example of a domain of dependence stretching in counter-flow direction from the internal face at position xj to a position %o _k being inside the i-3 ’th finite control volume.

Figure 3c) is a schematic representation of the same spatial reconstruction as shown in figures 3a) and 3b) and which also illustrates the mass being present in the domain of dependence, by the shaded area.

Figure 4 is a graphical presentation of a comparison of predicted fluid behaviour in a pipeline made with a prior art solver and a solver according to the computer- implemented method applied by the invention.

Figure 5 is a graphical presentation of a comparison of simulated liquid level (continuous liquid plus bubbles) in a second comparison test of a stretch from 1000 to 2000 meter of a pipeline.

Figure 6 is a graphical presentation of a comparison of simulated liquid level (continuous liquid plus bubbles) in a third comparison test of a stretch from 1000 to 2000 meter of a pipeline.

Figure 7 is a graphical presentation of negative volume fractions caused by the minmod limiter used instead of the PPS, in the third comparison test.

Figure 8 is a graphical presentation of a comparison of simulated volume fraction of continuous liquid, plotted on the top row, and volume fraction of bubbles, plotted on the bottom row, of a fourth comparison test.

Figure 9 is a graphical presentation of a predicted flow using a solver not using the numerical scheme according to the invention and which the source term for the hydraulic force is deactivated to give the same velocity for the gas and liquid phase.

Figure 10 is a graphical presentation of a predicted flow using a solver using the numerical scheme according to the invention and which the source term for the hydraulic force is deactivated to give the same velocity for the gas and liquid phase.

Figure 11 is a drawing schematically illustrating an example embodiment of a flow management system according to the invention. Figure 12 is a drawing schematically illustrating another example embodiment of a flow management system according to the invention.

Verification of the invention

The numerical scheme applied by the computer-implemented method applied by the invention according to the invention enables predicting multiphase flows in pipelines with an improved (lower) numerical diffusivity as compared to prior art numerical models. Simulations of hydrodynamic wave growth leading to plug flow (when the wave crests reach the upper parts of the pipe wall and slows the liquid flow) are particularly sensitive to numerical diffusion. At the same time, slug flows involve relatively rapid moving fluid phases involving relatively high pressure and mass gradients which requires relatively fine mesh sizes and relatively small time steps making such simulations particularly computational heavy.

Comparison test 1

The fluid behaviour in a 400 m long (linear) pipe having an internal diameter of 0.194 m and inclined 1° downward, and which is supplied with gas corresponding to a superficial velocity of 3.8 m/s and a liquid corresponding to a superficial velocity of 1 m/s at pressure 20 bar is predicted with a prior art solver utilising an implicit numerical scheme in which the mass fluxes are made partially explicit and compared with a similar prediction made with an implicit-explicit solver (implicit in pressure and velocity, explicit in mass and momentum convection) having incorporated the numerical scheme according to the first aspect of the invention. These flow conditions are known to cause wave growth towards the end of the pipe.

By performing the predictions with a relatively coarse grid, the numerical diffusion will dominate the physical models of the CFD-code and thus enabling assessing which meshes on the prior art and the present solver yield similar results by comparing the distance required for onset of the wave build-up as a measure of the accuracy of the numerical schemes.

The prior art solver was run three times with a mesh size (Axi) equal to 0.5, 1, and 2 timed the pipe diameter D, respectively. The result is shown graphically in figure 4 under the header “old solver”. The predictions gave a stratified flow at the downstream end of the pipe and an onset of wave build-up at 200 to 300 metres downstream into the pipe. The solver according to the first aspect of the invention obtained the same onset of the wave build-up when applying a mesh size of 7, 10, and 14 times the pipe diameter. These results are also shown in figure 4 under the header “new solver”.

From figure 4 we see that a solver according to the invention may use a mesh in the order of 10 times coarser to yield a similar onset of wave growth as a prior art solver, which represents a significant reduction in the computational load. The comparison is rather qualitative, but still gives an order of magnitude of the reduction in numerical diffusion, which gives increased numerical accuracy. Note that this case only involves hydrodynamic wave growth, which is particularly sensitive to numerical diffusion. One cannot necessarily expect such mesh ratios in all types of flows, since for example terrain slugging is not particularly sensitive to numerical diffusion. Still, even in cases dominated by terrain slugging, it might still be important to resolve hydrodynamic slugs, thus setting a minimum requirement on the refinement of the mesh.

Comparison test 2

This test investigates the numerical load when using an explicit solver having implemented the numerical scheme according to the first aspect of the invention to simulate the fluid behaviour in a 2 km long stretch and duration of 500 seconds of the same pipe fed with the same multiphase flow as described above for comparison test 1. The solver applies a second order polynomial for spatial reconstruction of the mass in the finite control volumes and a second order limiter function to rescale the spatial reconstruction to preserve positivity of the mass over the computational domain.

As a comparison, the same explicit solver is applied without the numerical scheme according to the first aspect of the invention but instead applies a higher-order upwind scheme using the minmod limiter function [Ref.3, page 110], A second comparison is also made with the same explicit solver without the numerical scheme according to the first aspect of the invention but a higher-order upwind scheme using the monotonized central-difference (MC) limiter function [Ref. 3, p.112]. MC limiter is expected to give more accurate results than minmod, at the risk of being numerically unstable. Amongst other, positivity of mass is never assured with higher-order limiters, but minmod is very mild and safe. MC is recommended in the same reference to be "a good default choice for a wide class of problems".

Without the numerical scheme according to the first aspect of the invention, it is necessary to apply a CFL condition defined with the true numerical velocities (the actual u _k j used in the scheme), set to be CFL < 0.8. In addition, the CFL may be violated during or at the end of the time step. This cannot be allowed, and if CFL > 1 during or at the end of the time step, the time step is restarted with a smaller At. Time step restarts can happen even with the numerical scheme according to the first aspect of the invention, for other reasons. Even though the numerical scheme according to the first aspect of the invention removes the CFL restriction related to mass convection, other physical phenomena which may be implemented in the numerical model may require strict time step restriction to preserve numerical stability (e.g. for surface waves), or less strict restriction to preserve modelling accuracy (e.g. friction of phase change). Thus, a CFL-condition is still applied in these example numerical calculations, but less strictly and only at the beginning of the time step. A local violation of the CFL condition is acceptable - avoiding restriction of the time step for the whole pipe due to one too high velocity - as is also a violation of the CFL condition at the end of the time step, and as such a time step restart can be avoided.

The resulting computational load is shown in table 1. Table 1 compares the time steps applied in the numerical simulations, number of restarts during the simulation, total used CPU time and time used to solve the mass transport and pressure equations (marked as “Time in mass conv.”). The column marked “PPS” is the simulation with the numerical scheme according to the first aspect of the invention. The column marked “NO PPS - minmod” is the first comparison simulation without the numerical scheme according to the first aspect of the invention but applying the minmod limiter function and the column marked “No PPs - MC” is the second comparison simulation without the numerical scheme according to the first aspect of the invention but applying the MC limiter function. The mesh size in all three simulations was 5 times the pipe diameter.

Figure 5 is a graphical presentation of the simulated liquid level (continuous liquid plus bubbles) in the stretch from 1000 to 2000 meter of the pipe for all three simulations. The “left box” marked “PPS” presents the simulated liquid level along the pipe segment obtained from the explicit solver having implemented the numerical scheme according to the first aspect of the invention. The “middle box” marked “No PPS - minmod” presents the simulated liquid level along the pipe segment obtained from the explicit solver using the minmod limiter function without the numerical scheme according to the first aspect of the invention, while the “right box” marked “No PPS - MC” presents the same for the simulation with the explicit solver using the MC limiter function without the numerical scheme according to the first aspect of the invention.

Table 1 Comparison of numerical workload The explicit solver having implemented the numerical scheme according to the first aspect of the invention (PPS) runs with the same average time step as the explicit solver using the minmod limiter function (minmod). This means that the ability to avoid time step restriction due to mass convection is not playing a role in this comparison. Most of the time, however, this ability is useful in case of velocity spikes at the slug/bubble transition, which is an artifact of the underlying physical models. The time step is considerably lower with the explicit solver using the MC limiter function (MC), as well as the number of time step restarts, indicating that there are probably velocity spikes with this scheme.

The numerical scheme according to the first aspect of the invention is relatively computationally heavy, as can be observed in the lower row of Table 1. Despite running with the same time step, the CPU time of the numerical scheme according to the first aspect of the invention (PPS) is 18% higher than with the explicit solver using the minmod limiter function (minmod). The difference in run time is coming from the execution time of the PPS, which is what makes the difference in the row 'Time in mass conv' of Table 1. The comparison of the CPU time with MC is however in favor of the PPS, due to the lower average time step, resulting in a higher total number of time steps to solve, as well as to the number of time step restarts, which requires to start over with a new smaller time step.

Comparison test 3

This test is similar to comparison test 2 above except for applying the same three numerical simulations on a 2000 meters long linear pipe having an internal diameter of 0.290 m and inclined 1° upward, and which is supplied with gas corresponding to a superficial velocity of 1.5 m/s and a liquid corresponding to a superficial velocity of 0.075 m/s at pressure 20 bar. The results are given in table 2 and in figure 6, respectively.

In this case, the explicit solver using the MC limiter function (MC) was not able to complete the simulation, due to instability issues. The explicit solver using the minmod limiter function ran with a time step close to twice as small as the explicit solver having implemented the numerical scheme according to the first aspect of the invention. The consequence is that the explicit solver with the PPS in this case run 25% faster, i.e. a total CPU time of 94 s versus 124 s, even though it is spending 20 s more, 26 s versus 6 s, in run time solving the mass transport equations as compared to the No PPS run. Table 2

These results indicate that “No PPS - minmod” seems to be more diffusive than PPS, as suggested by the fact that many waves are not reaching the top of the pipe and becoming slugs. In addition, negative masses are observed with minmod as limiter, even though the minmod limiter is a mild higher-order limiter. Figure 7 shows an example of this in a zoom on the simulated gas volume fraction using the No PPS - minmod scheme (middle plot in Figure 6). With the PPS, no negative mass can be observed at any time step.

Comparison test 4

This test is a 1000 m pipe of diameter 0.1 m, initialised with single-phase gas at rest. Then, liquid is injected at the left node at a superficial velocity of 1 m/s. Both liquid and gas are incompressible. We expect to see a liquid front between singlephase gas and single-phase liquid propagate at a velocity of 1 m/s. The same numerical schemes as in tests 2 and 3 are applied (PPS, No PPS - minmod and No PPS - MC, referring to a solver having implemented the numerical scheme according to the first aspect of the invention, a solver using the minmod limiter function without the numerical scheme according to the first aspect of the invention, and a solver using the MC limiter function without the numerical scheme according to the first aspect of the invention, respectively).

After 900 s, the results are plotted in Figure 8. The volume fraction of continuous liquid is plotted on the top row, and the volume fraction of bubbles is plotted on the bottom row. With the PPS applied, there are no volume fractions below 0 or above 1, while both minmod and MC limiters cause negative volume fractions (the sum of all volume fractions - continuous gas and liquid, bubbles and droplets at a given position in the pipe are always equal to 1, plus or minus the convergence criterion - equal to le-3 in the present case). This is a well-known behaviour of higher-order limiters to cause oscillations at discontinuities if they do not degenerate to first order fast enough. Minmod being more cautious than MC, it causes milder oscillations. The result of oscillations is that negative volume fractions are predicted. We can see here that the PPS keeps all volume fractions positive. Comparison test 5

For this test, the source term for the hydraulic force (pressure gradient due to a gradient in liquid level) is deactivated. Thus, if the gas and liquid velocities are exactly equal, the liquid-gas interface is expected to be advected as is at the same velocity. It is run with a solver having implemented the numerical scheme according to the first aspect of the invention, with a mass reconstruction using second order polynomials. A 300 m flat pipe is meshed with 318 cells, with an initial velocity of 5.57 m/s both in gas and in liquid, and an initial crenel-shaped gas-liquid interface as shown in Figure 9, left plot. The liquid and gas densities are respectively 818.7 kg/m ³ and 64 kg/m ³ . The inlet boundary condition is defined to inject the phases with exactly the right mass rate to keep the same volume fractions and velocities. In Figure 9, right plot, the result after advection during 30 s is plotted. The crenel-shaped interface is advected and somewhat similar to the initial shape, but spurious oscillations at both discontinuities have appeared. In Figure 10, the result is plotted after running the test with a solver having implemented the numerical scheme according to the first aspect of the invention, with a mass reconstruction using second order polynomials, and a limiter to avoid spurious oscillations. On the right plot, the spurious oscillations have disappeared. Numerical diffusion has smoothed the discontinuities to some extent, as expected. The crenel shape is otherwise advected as is. The small waves at ca. 150 m are small disturbances caused by the inlet node at the start of the case, which have been advected.

References

1 Veersteg, H. K. and Malalasekera, W., sec. Ed. (2007), “An Introduction to Computational Fluid Dynamics”, Pearson Education Limited,

ISBN: 978-0-13-127498-3

2 Chi-Wang Shu, “Positivity-preserving high order schemes for convection dominated equations”, presentation held at Brown University and retrievable on the internet: https://cfd.ku.edu/JRV/Shu.pdf

3 Randall J. Leveque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002, ISBN 978-0-521-00924-9

4 M. Griebel and M. Klitz, “CLSVOF as a fast and mass-conserving extension of the level-set method for the simulation of two-phase flow problems”, NUMERICAL HEAT TRANSFER, PART B 2017, VOL. 71, NO. 1, 1-36 http://dx.doi.org/10.1080/10407790.2016.1244400