Improvement of a multiphase flow model for wellhead chokes under critical and subcritical conditions using field data

The characterization of the multiphase flow through valves and orifices is a problem yet to be solved in engineering design, and there is a need for a prediction model able to simulate the complexity of this kind of flow in relation to fluid thermodynamic behaviour, and applicable to different incoming stream conditions and compositions. The present paper describes the development of a global model for the calculation of the discharge coefficient of orifices and choke valves operating under two- and three-phase flow as well as critical and subcritical conditions. The model generalizes the hydrovalve model developed by Selmer-Olsen et al. (in: Wilson (ed) Proceedings of 7th international conference on Multiphase Production, BHR Group, pp 441–446, 1995) and the Henry–Fauske (J Heat Transfer 93: 179–187, 1971. https://doi.org/10.1115/1.3449782) non-equilibrium model on the basis of an updated definition of the discharge coefficient. The model has been adapted to real choke valve geometries, by fitting the discharge coefficient and model parameters using field data from three production wells. The model developed is a global quartic function with different constants for the different valve geometries. The new discharge coefficient allows to simulate field data with high accuracy.


Introduction
Despite the energy transition from fossil fuels to renewable energy is started in the last decades, fossil fuels will continue to play an important role in a mix energy scenario; the optimization of the exploitation of old and new wells, including the transport of the operating fluids, became a vital issue for 1 3 the whole oil and gas industry. A common problem to be faced in order to optimize the production is a proper design of equipment with restrictions operating under multiphase flows, such as flow metering systems, choke valves, orifices, etc. (Dehkordi et al. 2017). In the present paper, the multiphase flow through orifices and wellhead chokes is studied. Wellhead chokes are widely used in the oil and gas to control the production flow rate from wells, to maintain stable pressure downstream the choke, and to provide the necessary backpressure to a reservoir to avoid formation damage from excessive drawdown (Sachdeva et al. 1986). Consequently, it is extremely important for petroleum engineers to be able to select the correct choke size for a given application (Richardson et al. 2006). This is only possible through accurate modelling of choke performance (Al-Attar 2008).
The multiphase flow through chokes is a complex process characterized by fluid properties, flow conditions, properties of the restriction, etc. From a practical point of view, this complex system could be simplified: for the purpose of modelling, a wellhead choke can be treated as an orifice in a pipe characterized by a coefficient which takes into account differences between real and theoretic discharge (Sachdeva et al. 1986).
Different methods and approaches for a proper design of a multiphase choke can be found in the literature. A lot of experiments have been carried out in order to understand the complex thermodynamics of the fluid in restrictions. For example, Fossa and Guglielmini (2002) experimentally investigated the two-phase flow pressure drop through singularities such as thin and thick orifices. Balakhrisna et al. (2010) conducted experiments to study the change of flow patterns during the simultaneous flow of high viscous oil and water through the sudden contraction, such as Colombo et al. (2013). Another recent approach to proper design multiphase valves is to conduct computational fluid dynamic (CFD) simulations: in particular, these types of simulations can be used to gain confidence and verify the performance of equipment (Malavasi and Messa 2014). However, these are simulations very onerous from a computational time point of view. In the flow assurance sector, there is need to be competitive in addressing sizing problems with simplified numerical approaches and this can be obtained through analytical models. In fact, in the last 40 years, several tens of models related to the calculation of a multiphase mixture flow rate flowing across orifices have been published (Chisholm 1983;Leung 1986Leung , 1990Leung , 1995Leung , 1996Darby 2000;Darby et al. 2001), each of which is based on specific set of assumptions that may be valid for certain specific conditions, but may not be accurate for others. However, it seems that no model exists, able to accurately predict the relevant flow characteristics for the whole possible operating envelope of typical oil and gas industry working conditions (Boccardi et al. 2005). An inaccurate model could lead to undersizing the valve (the valve would have insufficient capacity) or to its oversizing (the flow rate and the pressure drop through the valve would be higher than expected). Both these situations can negatively affect the stable operation of the valve, with possible damages.
The sizing of multiphase choke valves is quite complex because of: • phase change due to pressure drop. In particular, the liquid can evaporate and the vapour quality increment brings to a further increase of the pressure drop; on the other side, a dry gas can condensate due to the cooling across the valve (Joule-Thomson effect); • slip between gas and liquid phases; • possible marked non-equilibrium conditions between phases; • the velocity of sound does not follow a linear trend (in particular the sonic velocity in multiphase mixtures can be very much smaller than that of either of its components).
As reported by Giacchetta et al. (2014), the models available in the literature can be classified as homogeneous and non-homogeneous models. In the homogeneous models, the mixture is considered as a single fluid with thermodynamic and physical properties obtained by averaging different phase characteristics and considering the same speed for all phases. In the non-homogeneous approach, these assumptions are not valid. The homogeneous class can be further divided into equilibrium and non-equilibrium models. In the first one, a thermodynamic equilibrium between the phases during the passage through the orifice is considered. In the second, the non-equilibrium effects derive from the possible phase change delay with respect to the fast fluid transition time across the orifice, where pressure and temperature can exhibit huge variations in a short distance. In the oil and gas field, where the operating fluids are characterized by very different thermodynamic properties, it is normal practice to use non-equilibrium models.
The 1D thermo-hydraulic simulation tool OLGA (by Schlumberger) includes a multiphase flow choke model developed by Selmer-Olsen et al. (1995). The choke model in OLGA calculates the pressure drop and the critical flow rate, considering a circular-symmetric flow geometry. The pressure drop model is called hydrovalve. The choke model uses mixture balance equations for mass, momentum, and energy. Compression of gas into the narrow throat is accounted for in the model. Two thermodynamic nonequilibrium models are available in OLGA to calculate the phase fractions in the choke: • Frozen flow model • Henry-Fauske model In the Frozen model, it is assumed that there is no mass transfer between two phases during the efflux, i.e. the title is "frozen" to the valve upstream conditions value. This flow model is based on the following assumptions: equal velocities of phases, no heat or mass transfer between the phases and the gas (or vapour) is modelled as a perfect gas (Chisholm and Watson 1966).
In the Henry and Fauske Model (Henry and Fauske 1971), a phase change from upstream to throat conditions is considered, by assuming a partial thermodynamic equilibrium.
This paper presents a new numerical model addressing the hydrocarbon two-phase choked flow through orifices and its extension to three-phase flows (including water) and subcritical conditions. It is based on a previous model developed in the past by some of the authors (Giacchetta et al. 2014), the Henry-Fauske model and the hydrovalve model. The new model is based on an algorithm presented below and it allows to establish multiphase valve sizing criteria on the basis of simplified numerical approaches.

Flow equations
Henry and Fauske proposed that during two-phase critical flow through nozzles, orifices and short tubes, there is no time for a complete change in quality to equilibrium during the acceleration to the throat, so that the gas fraction at the throat x t is expressed as a function of stagnation gas fraction x 0 (corresponding to upstream restriction) by: where with x E [-] = equilibrium void fraction, dx E dp 1 Pa = change in equilibrium gas fraction with pressure, N[-] = equilibrium parameter, Δp[Pa] = pressure drop from upstream conditions to the throat.
If N = 1, the model corresponds to the homogeneous equilibrium model (HEM), while if N = 0, it coincides with the homogeneous frozen model (El-Wakil 1971). Therefore, the quantity N describes the partial phase change occurring at the throat. In the low-quality region (1 − x 0 ≈ 1) for which this formulation was intended, it can be shown numerically that the mass transfer is dominated by the behaviour of the liquid phase and the correlation essentially describes the flashing of the liquid. The experimental results indicated that the critical flow rates are in relatively good agreement with (1) Δp the homogeneous equilibrium model for stagnation qualities greater than 0.1, so that N is set equal to unity when x 0 ≥ 0.1. A stagnation quality of 0.1 corresponds to throat equilibrium qualities ranging from 0.125 to 0.155 depending on the pressure. An average value of 0.14 was chosen, thus: Evaluating the liquid as incompressible and the vapour as an ideal gas and considering the flow to be homogeneous, Henry and Fauske obtained Eq. 4 for the specific flow rate: where n is the expansion exponent that, if the phases maintain equal temperatures, can be expressed by (Tangren et al. 1949): where Cp and Cv represent the specific heat capacity at constant pressure and at constant volume, respectively.
For what concerns the throat pressure p t , it can be evaluated starting from the inlet pressure p 0 (Henry and Fauske 1971): x Et 0.14 for x Et < 0.14 N = 1 for x Et ≥ 0.14 (4) The most important feature of the Henry-Fauske model is that this model requires only the knowledge of the inlet (stagnation) conditions. This model is used in the present paper as a background of a new model to properly simulate choke valves operating in multiphase flow through the 1-D thermo-hydraulic simulator OLGA by Schlumberger.
Of course, a model for the calculation of the pressure drops must be coupled to the Henry-Fauske model in order to completely characterize the flow across the restriction. This is done by the hydrovalve model implemented in OLGA and developed by Selmer-Olsen et al. (1995). For a fluid that flows from position 1 through the throat to position 2 (see Fig. 1), it describes the pressure-drop to throat by the Bernoulli and continuity equations (Eq. 12).
In the above, u m is the mixture velocity and the momentum density, ρ m , is calculated by Eq. 13: where k is the slip ratio, u g the gas velocity, u l the liquid velocity, x the gas mass fraction. The recovery after the throat is described by the momentum and continuity equations: where M is overall mass flow through the choke. Unlike other models, the hydrovalve model does not assume a frictionless and adiabatic flow through the choke and, in order to take into account the difference between an ideal mass flow rate and the real mass flow rate, it introduces a parameter called discharge coefficient C d . The critical flow through the choke is found at the maximum of Eq. 12. Differentiating Eq. 12 with respect to the pressure and combining with Eq. 15 yields the following relation for the critical flow M c , where the throat area, A t , is corrected with the choke discharge coefficient, C d , to find the minimum flow area

Material and methods
The innovative approach developed to simulate the hydrocarbon multiphase flow through orifices and choke valves is based on the combination of flow equations as per "Flow equations" section, and on a new correlation for the calculation of the discharge flow coefficient, required as input by OLGA, developed by the authors. In particular, the new correlation is an extension of the orifice model developed by Giacchetta et al. (2014), for the simulation of hydrocarbon two-phase choked flow through orifices, to three-phase flows (including water) and subcritical conditions. The extension of the previous orifice model has been carried out by using experimental data presented in Schüller's papers (2003Schüller's papers ( , 2006, and through comparison with numerical calculations performed by OLGA. As a final step, the model has been fitted to obtain real choke valves flow performance prediction, by introducing appropriate adjustments and corrections in order to obtain a new general mathematical model for the complete characterization of multiphase flow through valves and orifices. The new model is based on an algorithm presented below.
The work has been divided into two phases: Phase A and Phase B. During the first phase of the work, Phase A, the

Fig. 1 Choke scheme Throat Section t Usptream Section 1 D ownstream Section 2
High Turbulence Region 1 3 orifice discharge flow coefficient, to be used in the hydrovalve model, has been calculated at first by Chisholm's model (1983) and then corrected by the model developed by Giacchetta et al. (2014). It was found that the model developed by Giacchetta et al. (2014) for choked two-phase flow across orifices is not able to simulate also subcritical and threephase flows. Thus, a new discharge flow coefficient model, suitable also for subcritical and three-phase flows across orifices, has been determined starting from the previously calculated Chisholm's value; the new discharge coefficient allows to simulate three-phase flow with higher accuracy, obtaining an average error in the flow rate prediction of 6%. At the end of this first stage, a second part of the work (Phase B) has been accomplished where the new Phase A model has been tested against field data relevant to operating choke valves, provided by an oil company, and a further development has been carried out to obtain a new general model suitable for multiphase flow through different valve geometries.
Concerning the reference database for Phase A, by Schüller et al. (2003Schüller et al. ( , 2006, the fluids used were recombined oil from the Njord field in the North Sea, natural gas from the Kaarstoe terminal in Norway, and water with added salts to give typical produced-water properties. The reference work consists of two parts. In the first part (Schüller et al. 2003), a total number of 367 single-, two-, and three-phase tests were performed. Two different flow geometries resembling an orifice and cage choke and three different choke diameters were tested (11, 14, and 18 mm). To extend these data into the critical flow region, after the installation of a particular multiphase flow pump, 142 new single-, two-, and three-phase tests were performed in the second part (Schüller et al. 2006), making a total data set of 509 data points.
However, only the results related to the 11-mm diameter orifice were presented. The gas/oil/water mixtures composition shown in Table 1 (Schüller et al. 2003) has been used in this study. The fluid has been characterized using Soave-Redlich-Kwong (SRK) EoS (Equation of State) with the Peneloux correction factor (Table 2).
During Phase B, the data relevant to three production wells belonging to an oil production field located in the South of Europe, named W1 (Well 1), W2 (Well 2) and W3 (Well 3) were used. The relevant hydrocarbon gas and oil phase compositions are shown in Tables 3, 4 and 5. The fluids have been thermodynamically characterized using the PR (Peng-Robinson) EoS.

Model development
The development procedure of the general model splitted in two work phases is described in the following "Phase A" and "Phase B" sections.

Phase A
During Phase A, a total of 51 steady-state experimental data have been simulated according to the computational procedure described here below. All data are relevant to the abovementioned critical and subcritical two-or three-phase (gas, oil, water) flow collected by Schüller et al. (2003Schüller et al. ( , 2006. The purpose of the computational procedure is to determine the mass rate flowing through the orifice operating in two-or three-phase critical or subcritical conditions, in order to match the measured values. The algorithm is composed by an iterative calculation procedure, as shown in Fig. 2. The calculation of the discharge coefficient is started from the Chisholm's value (C d,Ch ), as indicated by Giacchetta et al. (2014); then this is iteratively adjusted to match the calculated mass flow rates with the experimental values. The steps of the procedure represented in the diagram of Fig. 2, as previously defined by Giacchetta et al. (2014), are the following: 1. input data needed to define the problem are the orifice inlet pressure and temperature, P 1 and T 1 , respectively, the discharge pressure, P 2 , the pipe diameter, the orifice diameter, fluid properties, pipe spatial discretization, pipe material properties, time integration parameters; 2. on the basis of upstream temperature, T 1 , pressure, P 1 , and mass fractions, x 1 , determine gas and liquid specific volumes, v g and v l , heat capacities, Cp g and Cp l through the Equation of State, and then the fluid expansion exponent, n, and the velocity ratio, K; 3. calculate the orifice discharge coefficient by using Chisholm's model; 4. set the calculated C d,Ch value in OLGA and run the simulation; 5. the simulator calculates the mass flow rate on the basis of hydrovalve and Henry-Fauske models; 6. compare the mass flow rate estimated by OLGA against the experimental results; 7. if M OLGA ≠ M exp , adjust the discharge coefficient with further simulations in such a way to match the calculated mass flow rates with the experimental values; 8. if M OLGA = M exp , then C d,New is found.
The orifice discharge coefficient determined in the previous study by Giacchetta et al. (2014), as depicted by Eq. 17 and Table 5, has been modified in this study to include the extended flow conditions as said above.

Phase B
In the second part of the model development, the Phase B, a new general discharge flow coefficient model, suitable for subcritical and critical two-and three-phase flow across chokes valves, has been developed.
Initially, data provided by the oil company have been used to assess if the model proposed at the end of Phase A is suitable to simulate real choke geometries. It was found that the above model is not able to properly simulate the hydrocarbon multiphase flow through real choke valves. Thus, the same numerical procedure implemented in Phase A has been adopted to develop a new model suitable for real choke geometries. Two different flow geometries of choke valves were tested: the first one is a cage type geometry and it is installed on well W1, while the second one is a duse choke geometry and is installed on wells W2 and W3.
Field operating data, relevant to the above-mentioned oil production gathering system located in the South of Europe, are summarized in Table 6 and they are related to 21 measurements. Measurements are given as follows: choke upstream and downstream fluid pressure (P1 and P2), choke upstream and downstream fluid temperature (T1 and T2), gas, oil and water flow rates. The choke upstream pipe size (17)   Giacchetta et al. (2014) START Input data to define the problem: geometry, thermodynamic properties, pressure and temperature boundary conditions.
Set P 2 / P 1 and T 1 The new choke model has the following characteristics: • it is able to simulate two-and three-phase flows; • it is able to simulate critical and subcritical flow conditions; • it provides the discharge flow rate as main output result, given upstream and downstream pressure.

Results and discussion
Phase A Table 7 provides the comparison between the first iteration mass flow rate estimated by OLGA, M OLGA , and the experimental results, M exp , as well as pressures, temperature, mass fractions and mass flow rates for each analysed flow case of Phase A. M OLGA is calculated by using C d,Ch as input for the simulator. From Table 7, it can be seen that, under the conditions investigated, the discharged mass from the orifice is underestimated for some cases and overestimated for other ones. The maximum error in the flow rate prediction is equal to 18.5%. In order to evaluate the discrepancy between experimental and numerical results, a parameter R, defined as the ratio between the first iteration mass flow rate provided by OLGA when using C d,Ch and the experimental mass flow rate, has been introduced and analysed as function of the main operating parameters. In Figs. 3, 4, 5 and 6, this ratio is represented as a function of the pressure drop across the choke (ΔP), the inlet gas mass fraction (x g1 ), the inlet oil mass fraction (x o1 ) and the inlet water mass fraction (x w1 ), respectively, for all the runs. It can be concluded that the ratio R is almost independent from the operating conditions. This means that a correlation on the discharge coefficient can be determined and applied within a typical orifice operating envelope, disregarding the particular operation characteristics. Then, the solving procedure described in "Phase A" section has been implemented and the discharge coefficient providing the best matching with experimental data, C d,New , has been found. These new values are presented in Table 8.  Among the flow scenarios of Table 7, 20 cases with threephase or subcritical flow conditions were selected, and the relevant flow discharge coefficient determined by the current study C d,New has been compared with the discharge coefficient C d,Cr calculated by the correlation as per Giacchetta et al. (2014). The comparison is shown in Table 9, where it can be seen that the model developed by Giacchetta et al. (2014) for choked two-phase flow across orifices is not able to simulate also subcritical and three-phase flows. Thus, a new discharge flow coefficient model, suitable also for subcritical and three-phase flow across orifices, is necessary.
From Table 8, it can be noticed that a linear dependence of C d,New on C d,Ch is not observed and a more complex form of C d,New is likely to be expected. To obtain a new correlation for this parameter, a nonlinear regression was applied to the actual data set. In order to develop a new discharge flow coefficient C d,New, suitable also for subcritical and threephase flow across orifices, the data have been correlated by means of the same quartic function reported in Eq. 17, where the new constants a, b, d, c, d and e are reported in Table 10.

Phase B
During Phase B, a total of 21 steady-state field data measurements on operating choke valves, have been investigated according to the computational procedure presented in Fig. 2.
All data are relevant to two-and three-phase (gas, oil, water) flow through choke valves. Table 11 provides the comparison between the first iteration mass flow rate estimated by OLGA, M OLGA , and the experimental results, M exp , as well as pressures, temperature, mass fractions and mass flow rates for each analysed flow case measured by the oil company. M OLGA is calculated by using C d,Ch as input for the simulator. From Table 11, it can be seen that, under the conditions investigated, the discharged mass from the choke valve is overestimated for all cases except for case W1-8 for well W1. For this case, the flow rate is underestimated also using a discharge flow coefficient equal to 1. The strange behaviour of the case W1-8 suggests that measurement errors may be present and   therefore it has not been considered for the fitting of the discharge coefficient model. As done for Schüller's data in Phase A, in order to evaluate the discrepancy between experimental and numerical results, the parameter R has been introduced and analysed as function of the main operating parameters. In Figs. 7, 8, 9 and 10, this ratio is represented as a function of the pressure drop across the choke (ΔP), the inlet gas mass fraction (x g1 ), the inlet oil mass fraction (x o1 ) and the inlet water mass fraction (x w1 ), respectively, for all the runs.
From the above Figures, it can be concluded that the ratio R is almost independent from the operating conditions. This means that a correlation for the discharge coefficient can be determined and applied within a typical choke geometry, disregarding the particular operation characteristics.
Then, the solving procedure shown in Fig. 2 has been applied and the discharge coefficient which allows to match experimental data, C d,New , has been found. These new values are presented in Table 12. Table 12 shows that the dependence between C d,New and C d,Ch is not linear; in fact, the quarter function expressed by Eq. 17, was found adequate for the relationship being sought. Different regression parameters are expected for the two geometries of chokes. In fact, the operating and geometrical characteristics of the first type of choke (cage choke valve, installed on well W1) are quite different with respect to the second one (duse choke valve, installed on wells W2 and W3). This agrees with what reported by Schüller et al. (2003Schüller et al. ( , 2006: by using the hydrovalve model to simulate two choke geometries, they found that each type of geometry is characterized by a specific discharge coefficient. The constants a, b, d, c, d and e of Eq. 17 are presented in Table 13 for each type of choke geometry.
In order to establish the accuracy of the new model, it has been applied to the available set of measurements relevant to the wells W1, W2 and W3; new choke valves discharge coefficient values, C d , have been calculated for inclusion in the OLGA input of updated simulations. Table 14 provides the comparison between the new iteration mass flow rate estimated by OLGA, M OLGA2 , and the experimental results, M exp also showing the values of the discharge coefficient calculated by the new model.
In this case also, the results presented are in good agreement with the results obtained by Schüller et al. (2003Schüller et al. ( , 2006: in fact, for the simulation of the flow through a simplified choke geometry, they found that a discharge coefficient of the order of 0.5 must be used; this is in agreement with the results shown in Table 14 for wells W2 and W3. For the valve installed on W1 instead, the correct discharge coefficient is lower; this is probably due to the specific geometry of the valve, showing a more complicated arrangement.

Conclusions
In the present study, a model for the calculation of the flow discharge coefficient of orifices and choke valves operating under hydrocarbon multiphase flow conditions has been developed. In a first phase of the work, a previous model for critical two-phase flow through orifices by Giacchetta et al. (2014) has been extended to simulate three-phase and subcritical flow conditions also. The relevant method is based on the application of the combined Henry-Fauske and hydrovalve models as implemented in the multiphase flow simulation tool OLGA, and by determining the flow discharge coefficient based on matching laboratory data for flow through orifices. Then, in the second phase of the work, the same model has been applied to two different choke   Test point (-) valve geometries, referred to real operating choke valves installed in an oil production field in South Europe.
A new flow discharge coefficient has been developed starting from the orifice model, using field operating data related to the two above-mentioned choke valve geometries, duse and cage type, respectively; data are relevant to three production wells called W1, W2 and W3.
The model developed is a global quartic function with different constants for the different valve geometries. Similar results were found by Schüller et al. (2003Schüller et al. ( , 2006: by using the hydrovalve model to simulate two choke geometries, they concluded that each type of geometry is characterized by a particular discharge coefficient. For the simulation of the flow through a simplified choke geometry, they found that a discharge coefficient of the order of 0.5 must be used; this is in agreement with the results shown in Table 14 for the type of chokes (duse choke valve) installed on wells W2 and W3. On the contrary, a much lower discharge coefficient has been found appropriate for the simulation of flow performance of cage choke valve installed on well W1, probably due to a more complex valve internal arrangement and associated flow path.
The new discharge coefficient allows to simulate field data with high accuracy, obtaining an average error equal to 6% on predicted flow rate (maximum absolute error equal to 25% only for one test case).
The new choke model has the following characteristics: • it is able to simulate two-and three-phase flows; • it is able to simulate both critical and subcritical flow conditions; • it provides the discharge flow rate as main output result, given upstream and downstream pressure and temperature.  Due to the strong dependency of discharge flow coefficient from the choke valve geometry, it is pointed out that the new developed correlations are referred to the types of valves investigated. The following recommendations are provided for future activities: a. check the new developed correlations with data on similar valves; b. for other different types of valves, an updated analysis should be carried out.
A spreadsheet associated with the model has been developed to calculate the discharge coefficient, used as input value to the OLGA choke valve model. The OLGA model calculates as output: • total mass flow rate; • mass flow rate for each phase; • choke outlet temperature.
Funding The authors received no specific funding for this work.

Compliance with ethical standards
Conflict of interest On behalf of all the co-authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.