ESP Lifted Oil Field: Basic Model for Control, and Comparison of Simulation Tools

1. Introduction

1.1. Background

Petroleum products have been key energy carriers for more than a century. Current focus on climate1 implies a change towards sustainable energy carriers. To succeed in this change, a transition period from the use of fossil fuel is necessary. In the transition, improved operation of petroleum production through model based optimal operation will be necessary. Petroleum production entails slow (reservoir; months) and fast (reservoir-to-separator; seconds) subsystems; this is a focus of on-going research project “DigiWell”2. Vertical transport of petroleum from oil well to surface requires sufficient pressure to counteract gravitational and friction forces. If the oil-well heel pressure is insufficient for such transport, either (i) gas is injected in the vertical pipe to “blow” the petroleum fluids to the surface [“gas lifted”], or (ii) an Electric Submersible Pump [ESP] is installed in the vertical pipe to increase the pressure [“ESP lifted”] sufficiently. Here, the focus is on the dynamics of transport from the reservoir formation to a surface manifold via an ESP, and further horizontal transport from the manifold to a separator.

Industrial simulation tools typically put main emphasis on the dynamics of the reservoir (time constant: months) and use steady state models for the reservoir-to-surface transport. This emphasis is inadequate for daily operation and control. Here, a simplified, yet complete, dynamic model for oil transport from reservoir to separator is discussed. The model provides an understanding of the dynamic behavior of such systems, and is suitable for industrial control design, as well as for control and petroleum production studies. Emphasis is put on a simple, yet stringent model development, while avoiding variable unit complexities.

1.2. Previous work

Sharma and Glemmestad [1] (see also Sharma [3]) provide a dynamic model of oil transport from reservoir to separator suitable for control design. Binder et al. [4] discuss an older model; other models typically are CFD models, which are too complex for control design.

Sharma’s model considers a case with 4 vertical pipes from oil reservoirs to a single manifold, with 2 horizontal pipes from the manifold to a single separator. Each vertical pipe has an ESP and a choke valve at a common manifold entrance; the pump speeds can be manipulated individually. The horizontal pipes have booster pumps to counteract friction effects. The original ESP model includes induction motors, but the dynamics of the pump actuator is fast, and is neglected in later work. Sharma and Glemmestad [1] provide a novel ESP model, a simple model for a booster pump, and use a valve model based on on the ANSI/ISA S75.01 standard3. The model in Sharma and Glemmestad [1] was re-structured and simplified in Lie [2], emphasizing dimensionless equipment models, and thereby eliminating some level of complexity in common industry models.

The model with ESP in Sharma and Glemmestad [1] is mainly relevant for the production of heavy oil. Several papers use this model in advanced industrial control studies, Krishnamoorthy et al. [5], Delou et al. [6], Santana et al. [7].

Mixtures of liquid oil and water form an emulsion when stirred (e.g., in a multi-stage ESP); for such emulsions, the viscosity — and hence the friction — varies dramatically with water content, Justiniano and Romero [8]. Sharma and Glemmestad [1] assume an unrealistic linear viscosity dependence on water fraction.

1.3. Structure of paper

Section 2 gives an overview of the transport system from oil reservoir via manifold to a separator, and key equipment models. Section 3 develops a simple mechanistic model of the system. Section 4 contrasts two modeling languages for simulation: Modelica and Julia’s ModelingToolkit. Section 5 illustrates model behavior and the use of modeling/simulation tools for analysis and control. Finally, Section 6 provides some conclusions.

2. System description

Production of a mixture of water and crude oil in liquid phase is considered, where evaporation of liquids is assumed negligible.

2.1. System topology

Oil production systems merge several boreholes from the same or different reservoirs through vertical pipes into a manifold. Normally, more than one horizontal transport pipe is needed from the manifold to a separator for sufficient transport capacity. Water is commonly injected into the manifold to reduce friction loss in the horizontal pipes; the added water is typically recycled from the separator, and is at close to production temperature. Figure 1 shows a system with $n_{\mathrm{w}}$ wells/vertical pipes via a common manifold to $n_{\mathrm{t}}$ transportation/horizontal pipes leading to the separator.

In Figure  Figure 1 , $p_{\mathrm{f}}^j$ is formation pressure for well j, $p_{\mathrm{h}}^j$ is bore-hole heel pressure, ${\dot{V}}_{\mathrm{h}}^j$ is volumetric heel flow rate, $A_{\mathrm{v}}^{-,j}$ is vertical cross sectional area below the ESP pump, ${\ell }_{\mathrm{v}}^{-,j}$ is the length of the vertical pipe below the ESP pump, $p_{\mathrm{p}}^{\mathrm{i},j}$ is inlet pressure to the (ESP) pump, $p_{\mathrm{p}}^{\mathrm{e},j}$ is effluent pressure after the (ESP) pump, $A_{\mathrm{v}}^{+,j}$ is vertical cross sectional area above the ESP pump, ${\ell }_{\mathrm{v}}^{+,j}$ is the length of the vertical pipe above the ESP pump, ${\dot{V}}_{\mathrm{c}}^{\mathrm{i},j}$ is volumetric flow rate at the inlet to the choke valve, $p_{\mathrm{c}}^{\mathrm{i},j}$ is pressure at the inlet to the choke valve, $p_{\mathrm{c}}^{\mathrm{e},j}$ is pressure at the effluent from the choke valve. Next, $p_{\mathrm{m}}$ is manifold pressure, while ${\dot{V}}_{\mathrm{w}}$ is water added to the manifold to reduce viscosity. At the outlet from the manifold, $p_{\mathrm{bp}}^{\mathrm{i},j}$ is the influent pressure to the booster pump (BP), $p_{\mathrm{bp}}^{\mathrm{e},j}$ is effluent pressure after the booster pump, ${\dot{V}}_{\mathrm{t}}^j$ is volumetric flow rate in a transport pipe from manifold to separator, $p_{\mathrm{s}}^{\mathrm{i},j}$ is pressure at the inlet to the valve into the separator, and $p_{\mathrm{s}}^{\mathrm{e},j}=p_{\mathrm{s}}$ is the separator pressure.

For simplicity, it is assumed that $A_{\mathrm{v}}^{-,j}=A_{\mathrm{v}}^{+,j}=A_{\mathrm{v}}$. All vertical pipes are assumed connected to the same manifold pressure $p_{\mathrm{m}}$; hence effluent choke pressure satisfies $p_{\mathrm{c}}^{\mathrm{e},j}=p_{\mathrm{c}}^{\mathrm{e}}=p_{\mathrm{m}}$ for all j. The influent pressure to the booster pumps, $p_{\mathrm{bp}}^{\mathrm{i},j}$ are all assumed to be equal to the outlet pressure from the manifold, and have the same value, $p_{\mathrm{bp}}^{\mathrm{i},j}=p_{\mathrm{bp}}^{\mathrm{i}}=p_{\mathrm{m}}$. Likewise, all transport pipes end up in the same separator: $p_{\mathrm{s}}^{\mathrm{e},j}=p_{\mathrm{s}}$ for all j.

2.2. Fluid properties

The petroleum fluid properties are important. Density $\rho$ varies with pressure p and temperature T, $\rho \left(p,T\right)$. Neglecting temperature dependence, and assuming constant isothermal compressibility ${\beta }_T$,4 $\rho \left(p\right)$ is given as

$$\rho ={\rho }_{0}exp\left({\beta }_T\left(p-p_0\right)\right)$$

(1)

where $\left({\rho }_0,p_0\right)$ is some reference state.

Defining water cut ${\chi }_{\mathrm{w}}$ as ${\chi }_{\mathrm{w}}\triangleq {\dot{V}}_{\mathrm{w}}/\dot{V}$: volumetric flow rate of water divided by total flow rate of the fluid, total density $\rho$ can be expressed as

$$\rho ={\chi }_{\mathrm{w}}{\rho }_{\mathrm{w}}+\left(1-{\chi }_{\mathrm{w}}\right){\rho }_{\mathrm{o}};$$

(2)

here, ${\rho }_{\mathrm{w}}$ and ${\rho }_{\mathrm{o}}$ are constant densities of pure water and crude oil, respectively.

In reality, water and crude oil have different isothermal compressibilities. Here, we simplify and assume an overall value for ${\beta }_{\mathrm{T}}$. Using data in Appendix A1, density $\rho$ varies ca. $10\mathrm{kg}/{\mathrm{m}}^3$ with pressure variation in the range 25–$225\mathrm{bar}$, Figure 2.

We thus assume constant density in pipes, but a pressure-dependent density will be assumed in the manifold.

Sharma and Glemmestad [1] propose a simple linear mixing rule for kinematicviscosity $\nu$:

$$\nu ={\chi }_{\mathrm{w}}{\nu }_{\mathrm{w}}+\left(1-{\chi }_{\mathrm{w}}\right){\nu }_{\mathrm{o}}.$$

(3)

With $\nu$ known, dynamic viscosity $\mu$ can be computed (if needed) as

$$\mu =\nu \rho .$$

The linear interpolation model of Eq. 3 is used here to ease comparison with results in Sharma and Glemmestad [1], even though it is not physically realistic [8].

2.3. Well-bore production

Total production from the reservoir (formation pressure $p_{\mathrm{f}}$) relates volumetric petroleum fluid rate ${\dot{V}}_{\mathrm{h}}$ at the well-bore heel as ${\dot{V}}_{\mathrm{h}}\propto p_{\mathrm{f}}-p_{\mathrm{h}}$, where $p_{\mathrm{h}}$ is heel pressure and the proportionality constant $C_{\mathrm{pi}}$ is the productivity index,

$${\dot{V}}_{\mathrm{h}}=C_{\mathrm{pi}}·\left(p_{\mathrm{f}}-p_{\mathrm{h}}\right);$$

$C_{\mathrm{pi}}$ is unit-dependent. Here, we instead propose a dimensionless form,

$${\dot{V}}_{\mathrm{h}}={\dot{V}}_{\mathrm{pi}}^{\mathrm{c}}\frac{p_{\mathrm{f}}-p_{\mathrm{h}}}{p_{\mathrm{pi}}^{\varsigma }}$$

(4)

where ${\dot{V}}_{\mathrm{pi}}^{\mathrm{c}}$ is the productivity index capacity in the same unit as ${\dot{V}}_{\mathrm{h}}$, and $p_{\mathrm{pi}}^{\varsigma }$ is scaling pressure with the same unit as $p_{\mathrm{f}},p_{\mathrm{h}}$.

2.4. Pump models

Electric Submersible Pump 

Pump models are often given as

$$\Delta p_{\mathrm{p}}=\rho g{h}_{\mathrm{p}};$$

(5)

here, $h_{\mathrm{p}}=h_{\mathrm{p}}\left(\dot{V},f_{\mathrm{p}}\right)$ is pump head with volumetric flow rate $\dot{V}$ and control input $f_{\mathrm{p}}$ — rotational pump frequency $\mathrm{Hz}$.

Sharma and Glemmestad [1] provide values for minimal, maximal, and best-efficiency-point flow rates,

$$\begin{array}{cc}\hfill \frac{{\dot{V}}_{min}}{{\dot{V}}_{min,0}}& =\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{{\dot{V}}_{max}}{{\dot{V}}_{max,0}}& =\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{{\dot{V}}_{\eta }}{{\dot{V}}_{\eta ,0}}& =\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}.\hfill \end{array}$$

(8)

In Sharma and Glemmestad [1], a comprehensive model for the pump head of a multi-stage ESP is developed. To ease change of units, their model is here rewritten in dimensionless form

$$\begin{array}{cc}\hfill \frac{h_{\mathrm{p}}\left(\dot{V},f_{\mathrm{p}}\right)}{h_{\mathrm{p},0}}& ={\left(\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\right)}^2+\sum _{j=1}^3{a}_j{\left(\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\right)}^{2-j}{\left(\frac{\dot{V}}{{\dot{V}}^{\varsigma }}\right)}^j.\hfill \end{array}$$

(9)

In Eq. 9, $h_{\mathrm{p},0}$ is a nominal scaling head, $f_{\mathrm{p}}$ is the pump rotational frequency in the same unit as that of the nominal rotational frequency $f_{\mathrm{p},0}$, $\dot{V}$ is the actual volumetric flow rate out of the pump, ${\dot{V}}^{\varsigma }$ a scaling flow rate, and $a_1,\dots ,a_3$ are dimensionless model parameters5. Sharma and Glemmestad [1] include a head curve plot; the result in Figure 3 based on a dimensionless model is identical to their plot.

In addition, Sharma and Glemmestad [1] provide a model for the mechanical power requirement ${\dot{W}}_{\mathrm{p}}^{\mathrm{m}}={\dot{W}}_{\mathrm{p}}^{\mathrm{m}}\left(\dot{V},f_{\mathrm{p}}\right)$ for operating the pump6, again rewritten in dimensionless form,

$$\begin{array}{cc}\hfill \frac{{\dot{W}}_{\mathrm{p}}^{\mathrm{m}}}{{\dot{W}}_{\mathrm{p},0}^{\mathrm{m}}}& ={\left(\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\right)}^3+\sum _{j=1}^4{b}_j{\left(\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\right)}^{3-j}{\left(\frac{\dot{V}}{{\dot{V}}^{\varsigma }}\right)}^j.\hfill \end{array}$$

(10)

In Eq. 10, ${\dot{W}}_{\mathrm{p},0}^{\mathrm{m}}$ is a nominal scaling power consumption to operate the pump, $b_1,\dots ,b_4$ are dimensionless model parameters, while $f_{\mathrm{p}}$ and $\dot{V}$ are as above.

The actual power added to the fluid is

$${\dot{W}}_{\mathrm{p}}=\Delta p_{\mathrm{p}}\dot{V},$$

(11)

which gives the efficiency as

$$\eta =\frac{{\dot{W}}_{\mathrm{p}}}{{\dot{W}}_{\mathrm{p}}^{\mathrm{m}}}=\frac{\Delta p_{\mathrm{p}}\dot{V}}{{\dot{W}}_{\mathrm{p}}^{\mathrm{m}}}$$

(12)

where it is assumed that ${\dot{W}}_{\mathrm{p}}$ and ${\dot{W}}_{\mathrm{p}}^{\mathrm{m}}$ have the same units. Sharma and Glemmestad [1] include an efficiency curve plot; the result in Figure 4 based on a dimensionless model is identical to their plot.

Booster pump 

For the booster pump in horizontal pipes, a simpler model is suggested in Sharma and Glemmestad [1], rewritten in dimensionless form as

$$\frac{\Delta p_{\mathrm{bp}}\left(f_{\mathrm{bp}}\right)}{\Delta p_{\mathrm{bp},0}}={\left(\frac{f_{\mathrm{bp}}}{f_{\mathrm{bp},0}}\right)}^2$$

(13)

Here, $\Delta p_{\mathrm{bp}}\left(f_{\mathrm{bp}}\right)$ is the pressure increase at the given pump frequency/speed $f_{\mathrm{bp}}$, in the same unit as $\Delta p_{\mathrm{bp},0}$ — which is the pressure increase at the nominal pump frequency $f_{\mathrm{bp},0}$.

Pump control input 

In reality, the pump speed ($f_{\mathrm{p}}$, $f_{\mathrm{bp}}$) is not a control input. Instead, a motor is used to control the torque applied to the aggregate of motor and pump.

In Sharma [3], a model of the induction motor driving the ESP is developed. The experience [3] is that the motor dynamics is much faster than that of the mechanical system; hence in most of his work, Sharma [3] neglects the motor dynamics. However, it is not clear whether Sharma [3] considers the mechanical dynamics of accelerating the pump itself. This dynamics would be described by the kinetic energy balance in rotational form (the “swing equation”), which can be written as

$$\frac{\mathrm{d}K}{\mathrm{d}t}=P_{\mathrm{i}}-{\dot{W}}_{\mathrm{p}}^{\mathrm{m}}$$

with kinetic energy K

$$K=\frac{1}{2}{J}_{\mathrm{t}}{\omega }_p^2,$$

$J_{\mathrm{t}}$ is the total moment of inertia for the pump, the motor, and a possible flywheel, while $P_{\mathrm{i}}$ is the input power from the motor (control input), and ${\dot{W}}_{\mathrm{p}}^{\mathrm{m}}$ is as in Eq. 10.

In Sharma [3], the moment of inertia is approximately given as $J\approx 71\mathrm{kg}{\mathrm{m}}^2$. It is not clear whether this is the motor moment of inertia or the total moment of inertia. However, using such a moment of inertia leads to a pump time constant which is still much faster than the dynamics of the flow ${\dot{V}}_{\mathrm{v}}$, etc., hence the pump/motor dynamics is neglected here, and for simplicity, it is assumed that $f_{\mathrm{p}}$ is a control input.

2.5. Valve models

Sharma and Glemmestad [1] base their valve models on the ANSI/ISA S75.01 standard7. Here, instead a dimensionless description is proposed with extension to a control input,

$$\dot{m}={\dot{m}}_{\mathrm{v}}^{\mathrm{c}}·f\left(u_{\mathrm{v}}\right)\frac{{\rho }_{\mathrm{i}}}{{\rho }_{\mathrm{e}}}\sqrt{\frac{\left(p_{\mathrm{i}}-p_{\mathrm{e}}\right)/p^{\varsigma }}{{\rho }_{\mathrm{i}}/{\rho }^{\varsigma }}}$$

(14)

where ${\dot{m}}_{\mathrm{v}}^{\mathrm{c}}$ is the valve mass flow rate capacity, $u_{\mathrm{v}}\in \left[0,1\right]$ is the valve control signal, $f:\left[0,1\right]\to \left[0,1\right]$ is the valve characteristics, ${\rho }_{\mathrm{i}},{\rho }_{\mathrm{e}}$ are influent and effluent densities, respectively, $p_{\mathrm{i}},p_{\mathrm{e}}$ are influent and effluent pressures, respectively, while ${\rho }^{\varsigma },p^{\varsigma }$ are scaling density and pressure, respectively.

2.6. Friction loss

The friction drop along the pipe can be given by the Darcy-Weisbach model8 as

$$\frac{\Delta p_{\mathrm{f}}}{\ell }=f_{\mathrm{D}}\frac{\rho }{2}\frac{v^2}{D}$$

(15)

where $f_{\mathrm{D}}$ is Darcy’s friction factor given by Colebrook’s9 implicit expression. One explicit approximation to Colebrook’s expression is due to Swamee and Jain [9],

$$\frac{1}{\sqrt{f_{\mathrm{D}}}}=-2·{log}_{10}\left(\frac{5.74}{N_{\mathrm{Re}}^{0.9}}+\frac{ϵ/D}{3.7}\right),$$

(16)

where $N_{\mathrm{Re}}$ is the Reynolds number,

$$N_{\mathrm{Re}}=\frac{\rho vD}{\mu }=\frac{vD}{\nu },$$

(17)

$\mu$ is dynamic viscosity, $\nu$ is kinematicviscosity, and $ϵ$ is the “roughness height” of the pipe internal surface. Linear velocity v is related to volumetric flow rate $\dot{V}$ by

$$\dot{V}=vA$$

(18)

where A is the cross-sectional area of the pipe.

2.7. Why dimensionless models?

Example: ESP pump model 

As a first example, consider the ESP model in Eq. 9. In the original formulation in Sharma and Glemmestad [1],10

$$h_{\mathrm{p}}={\overline{a}}_0\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}+{\overline{a}}_1\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\dot{V}+{\overline{a}}_2{\dot{V}}^2+{\overline{a}}_3\frac{f_{\mathrm{p},0}}{f_{\mathrm{p}}}{\dot{V}}^3,$$

(19)

where parameters ${\overline{a}}_j$ have rather complicated units and the equation is hard-coded to assume11 $\left[\dot{V}\right]=\mathrm{gal}/\min$, while $\left[h_{\mathrm{p}}\right]=\mathrm{ft}$. In practice, either the rest of the model has to be posed using these units, or one has to operate with several copies of variables, e.g., ${\dot{V}}_{\mathrm{g}/\mathrm{m}}$ and $\dot{V}$, and remember to correctly convert between these versions of the flow rate. Both of these approaches are error-prone, and also require several versions of variables.

A much better solution is to write the model in dimensionless form. The simplest way to do this for the model in Eq. 19, is as

$$\begin{array}{cc}\hfill \frac{h_{\mathrm{p}}}{h_{\mathrm{p},0}}& ={\tilde{a}}_0\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}+{\tilde{a}}_1\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\frac{\dot{V}}{{\dot{V}}^{\varsigma }}+{\tilde{a}}_2{\left(\frac{\dot{V}}{{\dot{V}}^{\varsigma }}\right)}^2+{\tilde{a}}_3\frac{f_{\mathrm{p},0}}{f_{\mathrm{p}}}{\left(\frac{\dot{V}}{{\dot{V}}^{\varsigma }}\right)}^3.\hfill \end{array}$$

(20)

If we choose $h_{\mathrm{p},0}\equiv 1\mathrm{ft}$ and ${\dot{V}}^{\varsigma }\equiv 1\mathrm{gal}/\min$, then ${\tilde{a}}_j\equiv {\overline{a}}_j$. Suppose we want to generate the plot in Figure 3. Because that figure plots $h_{\mathrm{p}}$ in $\mathrm{ft}$ (the “native” unit), while the flow rate $\dot{V}$ is given in $\mathrm{bbl}/\mathrm{d}$ (“native” unit is $\mathrm{gal}/\min$), this result is produced by choosing ${\dot{V}}^{\varsigma }=1\mathrm{gal}/\min =34.29\mathrm{bbl}l/\mathrm{d}$, which can easily be found using the WolframAlpha app12, Figure 5.

In practice, it may be better to choose a more natural scaling unit, e.g., SI units. In that case, it is necessary to change parameters ${\overline{a}}_j$; for the parameters in Eq. 9, the scaling parameters are in SI units, where ${\overline{a}}_j\to a_j$ and $a_j$ is given in Table A2. To find $a_j$, choose $h_{\mathrm{p},0}=1\mathrm{ft}=0.3048\mathrm{m}$, ${\dot{V}}^{\varsigma }=1\mathrm{gal}/\min =6.309·{10}^{-5}{\mathrm{m}}^3/\mathrm{s}$, and write Eq. 20 as

$$\begin{array}{cc}\hfill \frac{h_{\mathrm{p}}}{\underset{\to h_{\mathrm{p},0}}{\underbrace{h_{\mathrm{p},0}{\tilde{a}}_0}}}& =\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}+\underset{a_1}{\underbrace{\frac{{\tilde{a}}_1}{{\tilde{a}}_0{\dot{V}}^{\varsigma }}}}\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\dot{V}+\underset{a_2}{\underbrace{\frac{{\tilde{a}}_2}{{\tilde{a}}_0{\left({\dot{V}}^{\varsigma }\right)}^2}}}{\dot{V}}^2+\underset{a_3}{\underbrace{\frac{{\tilde{a}}_3}{{\tilde{a}}_0{\left({\dot{V}}^{\varsigma }\right)}^3}}}\frac{f_{\mathrm{p},0}}{f_{\mathrm{p}}}{\dot{V}}^3,\hfill \end{array}$$

where the new $h_{\mathrm{p},0}$ is given in unit $\mathrm{m}$, $a_1$, $a_2$, $a_3$ are the new, dimensionless parameters in Eq. 9, while new scaling flow rate $\to {\dot{V}}^{\varsigma }=1{\mathrm{m}}^3/\mathrm{s}$. With this modified correlation ($a_1,\dots ,a_3$) for $h_{\mathrm{p}}$ using SI units, the dimensionless form is as in Eq. 9.

Example: control valve 

The ANSI/ISA S75.01 standard13 for compressible (i.e., ${\dot{m}}_{\mathrm{i}}={\dot{m}}_{\mathrm{e}}=\dot{m}$, ${\rho }_{\mathrm{i}}\ne {\rho }_{\mathrm{e}}$), non-choked fluids without fitting is

$$C=\frac{\dot{m}}{N_6\frac{{\rho }_{\mathrm{i}}}{{\rho }_{\mathrm{e}}}\sqrt{\frac{p_{\mathrm{i}}-p_{\mathrm{e}}}{{\rho }_{\mathrm{i}}}}}.$$

(21)

Here, C is the valve coefficient, $\dot{m}$ is the mass flow rate through the valve, ${\rho }_{\mathrm{i}}$ is the influent density, ${\rho }_{\mathrm{e}}$ is the effluent density, $p_{\mathrm{i}}$ is the influent pressure, $p_{\mathrm{e}}$ is the effluent pressure, $N_6$ is used to handle unit conversion. Typically, tabular values for $N_6$ are given which are valid for different combinations of units for $\dot{m}$, $\rho$, and p. This makes change of units rather complicated. A dimensionless formulation as in Eq. 14 greatly simplifies the use of the valve model in different units, and also includes a control valve characteristic.

3. Dynamic model

3.1. Balance laws

The model is based on the total mass balance (manifold) and the linear momentum balance (pipes). The total mass balance is expressed as

$$\frac{\mathrm{d}m}{\mathrm{d}t}={\dot{m}}_{\mathrm{i}}-{\dot{m}}_{\mathrm{e}}$$

(22)

where m is accumulated mass in the system, t is time, $\dot{m}$ is mass flow rate, and indices $\left(\mathrm{i},\mathrm{e}\right)$ denote influent and effluent, respectively.

The linear momentum balance is

$$\frac{\mathrm{d}\mathfrak{m}}{\mathrm{d}t}={\dot{\mathfrak{m}}}_{\mathrm{i}}-{\dot{\mathfrak{m}}}_{\mathrm{e}}+F,$$

(23)

where $\mathfrak{m}$ is linear momentum given as $\mathfrak{m}=mv$ with linear velocity v, $\dot{\mathfrak{m}}$ is momentum flow rate given as $\dot{\mathfrak{m}}=\dot{m}v$, and F is total force. With constant fluid density, ${\dot{\mathfrak{m}}}_{\mathrm{i}}={\dot{\mathfrak{m}}}_{\mathrm{e}}$, and the momentum balance reduces to Newton’s law, $\frac{\mathrm{d}\mathfrak{m}}{\mathrm{d}t}=F$.

3.2. Vertical pipes with ESP

We assume constant density in the pipes, causing volumetric vertical flow rate ${\dot{V}}_{\mathrm{v}}$ to be the same everywhere: ${\dot{V}}_{\mathrm{h}}={\dot{V}}_{\mathrm{c}}^{\mathrm{i}}={\dot{V}}_{\mathrm{v}}$. Furthermore, Eq. 23 reduces to Newton’s law. Momentum is given as $\mathfrak{m}=\dot{m}v$ with $\dot{m}=\rho {\dot{V}}_{\mathrm{v}}$, and v related to ${\dot{V}}_{\mathrm{v}}$ by Eq. 18. The total force is $F=F_{\mathrm{p}}+F_{\mathrm{b}}-F_{\mathrm{f}}-F_{\mathrm{g}}$, with

• Pressure forces at inlet and outlet of the pipe,

  $$F_{\mathrm{p}}=p_{\mathrm{h}}A-p_{\mathrm{c}}^{\mathrm{i}}A$$

  (24)
• Possible pressure boost due to a pump,

  $$F_{\mathrm{b}}=\Delta p_{\mathrm{p}}A,$$

  (25)

  with $\Delta p_{\mathrm{p}}$ given by Eqs. 5, 9,
• Friction loss,

  $$F_{\mathrm{f}}=\Delta p_{\mathrm{f}}A,$$

  (26)

  with $\Delta p_{\mathrm{f}}$ given by Eqs. 15, 16, 17, 18,
• Flow against gravity, with a vertical height h,

  $$F_{\mathrm{g}}=\Delta p_{\mathrm{g}}A,$$

  (27)

  with

  $$\Delta p_{\mathrm{g}}={\rho }_{\mathrm{v}}gh.$$

  (28)

In addition, we need information about how flow rate ${\dot{V}}_{\mathrm{v}}$ relates to the bottom hole pressure via the productivity index, Eq. 4, and how the flow rate ${\dot{V}}_{\mathrm{v}}$ relates to the choke valve flow, Eq. 14.

The most structured formulation would be to pose the momentum balance (here: Newton’s law) as the differential equation, and add all necessary algebraic equations. However, the OpenModelica DAE solver struggles with such a formulation: the valve equation Eq. 14 is implicit in pressure difference; in the iteration to find $\Delta p_{\mathrm{v}}=p_{\mathrm{i}}-p_{\mathrm{e}}$, if $\Delta p_{\mathrm{v}}$ becomes negative, the square root gives a complex number, and the simulation crashes.14 Instead, the differential variable has been changed to ${\dot{V}}_{\mathrm{v}}$; then the valve equation can be inverted and expressed as $\Delta p_{\mathrm{v}}\propto {\dot{V}}_{\mathrm{v}}^2$.

The following formulation is used in OpenModelica and ModelingToolkit:

$$\frac{\mathrm{d}{\dot{V}}_{\mathrm{v}}}{\mathrm{d}t}=\frac{p_{\mathrm{h}}-p_{\mathrm{i}}^{\mathrm{c}}+\Delta p_{\mathrm{p}}-\Delta p_{\mathrm{f}}-\Delta p_{\mathrm{g}}}{{\rho }_{\mathrm{v}}\ell /A}$$

(29)

$$\begin{array}{cc}\hfill {\rho }_{\beta }^0& ={\chi }_{\mathrm{w}}{\rho }_{\mathrm{w}}+\left(1-{\chi }_{\mathrm{w}}\right){\rho }_{\mathrm{o}}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \nu & ={\chi }_{\mathrm{w}}{\nu }_{\mathrm{w}}+\left(1-{\chi }_{\mathrm{w}}\right){\nu }_{\mathrm{o}}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \mu & ={\rho }_{\beta }^0\nu \hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{v}}& ={\rho }_{\beta }^{0}exp\left({\beta }_T\left(p_{\mathrm{c}}^{\mathrm{i}}-p_{\beta }^0\right)\right)\hfill \end{array}$$

(33)

$$p_{\mathrm{h}}=p_{\mathrm{f}}-p_{\mathrm{pi}}^{\varsigma }\frac{{\dot{V}}_{\mathrm{v}}}{{\dot{V}}_{\mathrm{pi}}^{\mathrm{c}}}$$

(34)

$${\dot{m}}_{\mathrm{v}}={\rho }_{\mathrm{v}}{\dot{V}}_{\mathrm{v}}$$

(35)

$$p_{\mathrm{i}}^{\mathrm{c}}=p_{\mathrm{m}}+p_{\mathrm{v}}^{\varsigma }\frac{{\rho }_{\mathrm{v}}}{{\rho }_{\mathrm{v}}^{\varsigma }}{\left(\frac{{\dot{m}}_{\mathrm{v}}}{{\dot{m}}_{\mathrm{v}}^{\mathrm{c}}}\right)}^2\frac{1}{f_{\mathrm{c}}^2\left(u_{\mathrm{v}}\right)}$$

(36)

$$\begin{array}{cc}\hfill h_{\mathrm{p}}& =h_{\mathrm{p},0}\left({\left(\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\right)}^2+a_1\frac{f_{\mathrm{p}}}{f_{\mathrm{p},0}}\frac{{\dot{V}}_{\mathrm{v}}}{{\dot{V}}^{\varsigma }}\right.\hfill \\ \hfill & \left.+a_2{\left(\frac{{\dot{V}}_{\mathrm{v}}}{{\dot{V}}^{\varsigma }}\right)}^2+a_3\frac{f_{\mathrm{p},0}}{f_{\mathrm{p}}}{\left(\frac{{\dot{V}}_{\mathrm{v}}}{{\dot{V}}^{\varsigma }}\right)}^3\right)\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill \Delta p_{\mathrm{p}}& ={\rho }_{\mathrm{v}}g{h}_{\mathrm{p}}\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill v_{\mathrm{v}}& =\frac{{\dot{V}}_{\mathrm{v}}}{A}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill N_{\mathrm{Re}}& =\frac{v_{\mathrm{v}}{d}_{\mathrm{v}}}{{\nu }_{\mathrm{v}}}\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill f_{\mathrm{D}}^{\mathrm{v}}& =\frac{1}{4{\left({log}_{10}\left(\frac{5.74}{N_{\mathrm{Re}}^{0.9}}+\frac{{ϵ}_{\mathrm{v}}/d_{\mathrm{v}}}{3.7}\right)\right)}^2}\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill \Delta p_{\mathrm{f}}& =\ell ·f_{\mathrm{D}}^{\mathrm{v}}\frac{{\rho }_{\mathrm{v}}}{2}\frac{v_{\mathrm{v}}^2}{d_{\mathrm{v}}}\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill \Delta p_{\mathrm{g}}& ={\rho }_{\mathrm{v}}gh.\hfill \end{array}$$

(43)

If we only consider the model of a single vertical pipe, we need to specify (i) initial state (e.g., ${\dot{V}}_{\mathrm{v}}$), (ii) all “input” variables, i.e., $p_{\mathrm{f}}$, $f_{\mathrm{p}}$, $p_{\mathrm{m}}$, and possibly water cut ${\chi }_{\mathrm{w}}$, and (iii) all parameters, i.e., ${\rho }_{\mathrm{w}}$, ${\rho }_{\mathrm{o}}$, ${\nu }_{\mathrm{w}}$, ${\nu }_{\mathrm{o}}$, $p_{\beta }^0$, ℓ, A, $p_{\mathrm{pi}}^{\varsigma }$, ${\dot{V}}_{\mathrm{pi}}^{\mathrm{c}}$, $p_{\mathrm{v}}^{\varsigma }$, ${\rho }_{\mathrm{v}}^{\varsigma }$, ${\dot{m}}_{\mathrm{v}}^{\mathrm{c}}$, $h_{\mathrm{p},0}$, $f_{\mathrm{p},0}$, ${\dot{V}}^{\varsigma }$, $a_1,$$a_2$, $a_3$, g, $d_{\mathrm{v}}$, ${\nu }_{\mathrm{v}}$, ${ϵ}_{\mathrm{v}}$, h.

3.3. Manifold

We assume a perfectly mixed manifold. Assuming constant manifold volume $V_{\mathrm{m}}$, and adding water at flow rate ${\dot{V}}_{\mathrm{w}}$ to dilute the fluid to a specified manifold water cut ${\chi }_{\mathrm{w}}^{\mathrm{m}}$, thus reducing friction loss in the pipe towards separator, ${\dot{V}}_{\mathrm{w}}$ must be approximately

$${\dot{V}}_{\mathrm{w}}=\frac{{\chi }_{\mathrm{w}}^{\mathrm{m}}-{\chi }_{\mathrm{w}}}{1-{\chi }_{\mathrm{w}}^{\mathrm{m}}}{\dot{V}}_{\mathrm{v}}.$$

(44)

Total mass balance for the manifold can then be expressed as

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{p}_{\mathrm{m}}}{\mathrm{d}t}& =\frac{1}{{\rho }_{\mathrm{m}}{V}_{\mathrm{m}}{\beta }_T}\left({\rho }_{\mathrm{v}}{\dot{V}}_{\mathrm{v}}+{\rho }_{\mathrm{w}}{\dot{V}}_{\mathrm{w}}-{\rho }_{\mathrm{m}}{\dot{V}}_{\mathrm{t}}\right)\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\rho }_{\beta }^0& ={\chi }_{\mathrm{w}}^{\mathrm{m}}{\rho }_{\mathrm{w}}+\left(1-{\chi }_{\mathrm{w}}^{\mathrm{m}}\right){\rho }_{\mathrm{o}}\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{m}}& ={\rho }_{\beta }^{0}exp\left({\beta }_T\left(p_{\mathrm{m}}-p_{\beta }^0\right)\right)\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\dot{V}}_{\mathrm{w}}& =\frac{{\chi }_{\mathrm{w}}^{\mathrm{m}}-{\chi }_{\mathrm{w}}}{1-{\chi }_{\mathrm{w}}^{\mathrm{m}}}{\dot{V}}_{\mathrm{c}}^{\mathrm{i}}\hfill \end{array}$$

(48)

In practice, the water cut ${\chi }_{\mathrm{w}}$ and flow rate ${\dot{V}}_{\mathrm{c}}^{\mathrm{i}}$ are not known perfectly, and it is necessary to use a feedback control system to manipulate ${\dot{V}}_{\mathrm{w}}$ instead of using Eq. 44.

For the manifold model, we must know (i) the initial manifold pressure, (ii) the vertical inflow ${\dot{V}}_{\mathrm{v}}$ and the horizontal transport flow ${\dot{V}}_{\mathrm{t}}$ from manifold to separator, as well as manifold water cut ${\chi }_{\mathrm{w}}^{\mathrm{m}}$, and (iii) parameters.

3.4. Transport pipe

For simplicity, we will neglect the separator inlet valve, and assume that $p_{\mathrm{s}}^{\mathrm{i},j}\equiv p_{\mathrm{s}}$. It is straightforward to reverse this assumption.

The model of the horizontal pipe from manifold to separator is almost identical to the vertical pipe from reservoir to manifold. The essential differences are (i) no gravity pressure drop, (ii) simpler booster pump model, (iii) neglecting pressure drop from pipe into separator, (iv) no need for a production index model. The complete model is

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\dot{V}}_{\mathrm{t}}}{\mathrm{d}t}& =\frac{p_{\mathrm{m}}-p_{\mathrm{s}}+\Delta p_{\mathrm{bp}}-\Delta p_{\mathrm{f}}^{\mathrm{t}}}{{\rho }_{\mathrm{t}}{\ell }_{\mathrm{t}}/A_{\mathrm{t}}}\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill {\rho }_{\beta }^{0,\mathrm{t}}& ={\chi }_{\mathrm{w}}^{\mathrm{m}}{\rho }_{\mathrm{w}}+\left(1-{\chi }_{\mathrm{w}}^{\mathrm{m}}\right){\rho }_{\mathrm{o}}\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill {\nu }_{\mathrm{t}}& ={\chi }_{\mathrm{w}}^{\mathrm{m}}{\nu }_{\mathrm{w}}+\left(1-{\chi }_{\mathrm{w}}^{\mathrm{m}}\right){\nu }_{\mathrm{o}}\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{t}}& ={\rho }_{\beta }^{0,\mathrm{t}}{\nu }_{\mathrm{t}}\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{t}}& ={\rho }_{\beta }^{0}exp\left({\beta }_T\left(p_{\mathrm{m}}-p_{\beta }^0\right)\right)\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill \Delta p_{\mathrm{bp}}& =\Delta p_{\mathrm{bp},0}{\left(\frac{f_{\mathrm{bp}}}{f_{\mathrm{bp},0}}\right)}^2\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill v_{\mathrm{t}}& =\frac{{\dot{V}}_{\mathrm{t}}}{A_{\mathrm{t}}}\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill N_{\mathrm{Re},\mathrm{t}}& =\frac{v_{\mathrm{t}}{d}_{\mathrm{t}}}{{\nu }_{\mathrm{t}}}\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill f_{\mathrm{D}}^{\mathrm{t}}& =\frac{1}{4{\left({log}_{10}\left(\frac{5.74}{N_{\mathrm{Re},\mathrm{t}}^{0.9}}+\frac{{ϵ}_{\mathrm{t}}/d_{\mathrm{t}}}{3.7}\right)\right)}^2}\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill \Delta p_{\mathrm{f}}^{\mathrm{t}}& ={\ell }_{\mathrm{t}}·f_{\mathrm{D}}^{\mathrm{t}}\frac{{\rho }_{\mathrm{t}}}{2}\frac{v_{\mathrm{t}}^2}{d_{\mathrm{t}}}.\hfill \end{array}$$

(58)

Again, we need to know the initial condition of the differential variable (${\dot{V}}_{\mathrm{t}}$), the inputs (${\chi }_{\mathrm{w}}^{\mathrm{m}}$, $f_{\mathrm{bp}}$, $p_{\mathrm{m}}$, $p_{\mathrm{s}}$), and the parameters.

3.5. Combined model

For illustration, we use two vertical pipes, one manifold, and one horizontal transport pipe from manifold to separator; Sharma and Glemmestad [1] use 4 vertical pipes, one manifold, and two horizontal transport pipes. Both Modelica and Julia’s ModelingToolkit have support for building classes/reusable models. Because of the similarity between the models for vertical and horizontal pipes, it would be possible to collect these in the same class/constructor and just differentiate between them with a function argument. The manifold model should be a separate class, though.

With re-usability of such classes/constructors, modeling of the combined system simply consists of (i) instantiating one model per unit (2 vertical pipes, one horizontal transport pipe, and the manifold), and (ii) connecting the various instances. Specifically, the vertical pipes should see the same manifold pressure $p_{\mathrm{m}}$, the vertical transport pipe should have the same inlet pressure as the manifold pressure $p_{\mathrm{m}}$, the influent volumetric flows to the manifold should be the sum of the flows from the vertical pipes and the viscosity diluting water feed ${\dot{V}}_{\mathrm{w}}$ now being

$${\dot{V}}_{\mathrm{w}}=\frac{{\sum }_{i=1}^2\left({\chi }_{\mathrm{w}}^{\mathrm{m}}-{\chi }_{\mathrm{w}}^i\right){\dot{V}}_{\mathrm{v}}^i}{1-{\chi }_{\mathrm{w}}^{\mathrm{m}}};$$

(59)

the effluent volumetric flow from the manifold is still ${\dot{V}}_{\mathrm{t}}$.

For a proper re-usable implementation, connections should be done using connectors (supported by both Modelica and ModelingToolkit). Connectors are not implemented here.

4. Simulation tools

Modelica is a mature language dating back to the 1990s; ModelingToolkit [MTK] is some 4–5 years old and is still evolving rapidly. MTK is more general than Modelica, and is also integrated in the larger eco-system of Julia. Currently, MTK does not support a graphical flow-sheeting tool, and it is unclear whether MTK allows for as large models as OpenModelica. Both tools have extensive support for building libraries.

The combined model has been solved using the free languages/tools OpenModelica [10,11] and ModelingToolkit [12] for Julia. To illustrate the similarity between OpenModelica code and a current formulation using ModelingToolkit, the following listing shows parts of the Modelica code for the reservoir heel–to–manifold; to save space, description of quantities is only included for constant $\pi$ to illustrate how it is done:

Next, the following listing shows similar parts of the ModelingToolkit code for the reservoir heel–to–manifold:

These listings show that Modelica code and ModelingToolkit code have a high degree of similarity. A few things to note:

• In Modelica, the independent temporal variable has a fixed name (time), and the time differentiation operator has a fixed name (der). In ModelingToolkit, both of these can be freely named by the user. In order to make unit models work together (e.g., in a standard library), it is, however, necessary to standardize on a name for time (commonly t); differentiation can be given a name as, e.g., Dt = Differential(t) or similar.
• In Modelica, quantities need to be specified with a type (e.g., Real), and are prepended with a qualifier (e.g., constant, parameter) — except for variables. For Julia and ModelingToolkit, the data type is inferred, unless explicitly stated. In the code above, quantities in MTK are grouped within begin...end blocks in macros (identifiers prepended by @, e.g., @parameters).
• Modelica has a simple way to handle implicit algebraic equations, and in many cases an initial guess of the algebraic variable is not required (see variable p_c_ _i in the Modelica code). In ModelingToolkit, initial values for unknowns after structural simplification (“states”) must be provided with numeric values (see variable p_c_ _i in the MTK code).
• In ModelingToolkit, initial values of differential variables can be changed outside of the code, hence default values can be written as Vd_v(t)=23.15e-3. In Modelica, only parameters can be changed outside of the code (after compilation), hence a parameter has been defined to hold the default initial value Vd_v(start = Vd_v0, fixed = true.
• Modelica uses symbol = for mathematical equality; MTK uses symbol ~ since Julia already uses symbol = for assignment.

The default solver in OpenModelica is excellent, although here it struggled if the model is posed as a DAE formulation with momentum as differential variable. ModelingToolkit can use solvers from the large, high quality DifferentialEquations.jl package [13]. With ModelingToolkit, more thought is currently required when choosing solver, accuracies, etc., compared to OpenModelica. Also, OpenModelica handles step-changes in inputs well, while for the DifferentialEquations.jl solvers, it is often necessary to specify the time points where step changes occur. On the other hand, the solutions from ModelingToolkit include interpolation functions, which yields smooth solutions with considerably fewer data points than for Modelica.

Results presented in Section 5 compare numerical solutions for the Reservoir heel–to–Manifold system for ModelingToolkit vs. OpenModelica.

OpenModelica’s support for linearization and plotting can be accessed from Julia via the OMJulia API [14]. ModelingToolkit is integrated in the Julia eco-system, with support for linearization, plotting, control systems analysis, random variables, etc., and has overall more possibilities that OpenModelica if further analysis is required.

OpenModelica is currently reported to handle models up to approximately ${10}^6$ variables/equations; various conference presentations indicate that ModelingToolkit currently can solve models of up to approximately ${10}^5$ variables/equations.15 The model above (Reservoir_2_Manifold) is reported by OpenModelica to have 15 variables (differential+algebraic) and 15 equations. In ModelingToolkit, linear equations with “observed” variables are stripped off from the model (function structural_simplify()) before solving the model. The above Reservoir_2_Manifold model in the listing is reported to have 2 “states” and 2 equations by ModelingToolkit. It is not clear whether the ModelingToolkit claim of ${10}^5$ variables/equations is before or after the “observed” variables are stripped off.

Other commonly used languages for scientific computing are MATLAB (commercial) and Python (free). Compared to both of these languages, Julia (free) has a more extensive set of differential equation solvers16. Neither MATLAB nor Python offer equation based modeling languages with library/re-use support such as Modelica or ModelingToolkit; MathWorks do offer Simscape17 (commercial) with MATLAB integration for such use, though.

5. Results

5.1. Reservoir heel to manifold

Parameters, initial conditions, and system inputs are given in Appendix A. Figure 6 shows the input variation for the Reservoir heel–to–manifold (R2M) case.

A step change in the formation pressure (red curve) as in Figure 6 is not very realistic; such changes are normally slow. The manifold pressure (blue curve) is not normally an input function, but rather a dependent variable in the overall system as in Section 5.2, and thus also normally varies slowly. The pump speed (green curve) is, however, a control variable, and can change fast. Still, the inputs in Figure 6 will help provide useful information about time constants in the system.

Figure 7 shows the response in (vertical) volumetric flow rate ${\dot{V}}_{\mathrm{v}}$, with comparison between Julia (red, solid) and OpenModelica (blue, dash-dot).

An important observation is that the volumetric flow rate is continuous under the step changes in Figure 6. This makes sense, in that the momentum of the fluid (oil-water) is substantial. Time constants are in the range of $0.2$ –$0.5\mathrm{s}$.

Figure 8 shows the response in the choke valve inlet pressure, $p_{\mathrm{c}}^{\mathrm{i}}$.

Observe that the choke valve inlet pressure normally is continuous under step changes, but that it changes discontinuously upon a step change in manifold pressure at $t=1.5\mathrm{s}$. Again, this makes sense: when the manifold pressure drops, Figure 6, the choke valve inlet pressure must also drop in proportion so that the flow through the valve changes continuously, see Eq. 14.

Figure 9 shows the response in the reservoir pressure heel, $p_{\mathrm{h}}$.

As noted, the formation pressure can not change in a step, but if it does, the well heel pressure must also change in proportion (i.e., discontinuously) to maintain continuity in the flow rate.

Figure 10 shows the response in the ESP pump pressure head, $\Delta p_{\mathrm{p}}$.

Again, the discontinuous change in pump pressure head is due to the step change in the manifold pressure, and the result makes sense.

So far, ModelingToolkit/Julia and OpenModelica have given (seemingly) identical simulation results, Figure 7, Figure 8, Figure 9 and Figure 10. Figure 11 shows the pipe friction loss, $\Delta p_{\mathrm{f}}$.

It is interesting to contrast the magnitude of the friction loss in Figure 11 and how much smaller it is than the pressure boost in the ESP, Figure 10. Obviously, if a more realistic viscosity model had been used, and in particular if emulsification occurs [8], the friction pressure drop might increase considerable with an unfortunate mixing fraction of oil and water.

Apart from this, it is interesting to observe a slight discrepancy between the result from ModelingToolkit/Julia and OpenModelica in this case. This discrepancy is maintained during a multitude of tests with different solvers and accuracy for the DifferentialEquation.jl solvers [13] and the solvers supported by OpenModelica. It seems like there is a minor discrepancy at $t=0$ for $\Delta p_{\mathrm{f}}$, which propagates throughout the solution. Because the codes and initial conditions are the same for both implementations, a natural suspicion is that the difference is due to different handling of the implicit algebraic equation at initial time, and that $\Delta p_{\mathrm{f}}$ is rather sensitive to such an inaccuracy.

5.2. Reservoir heel to separator

Parameters, initial conditions, and system inputs are given in Appendix A. For vertical pipe #2, scaling pump head $h_{\mathrm{p},0}$ is set to 80% of the value suggested in Appendix A. For this more complete system (2 vertical pipes, one horizontal transport pipe), there is no observed difference between the solution from ModelingToolkit/Julia and OpenModelica.

Figure 12 shows the input variation for the Reservoir heel–to–separator (R2S) case; because of slower dynamics for this larger system, the locations of the step changes have been changed, and the step change in the manifold pressure (Figure 6) has been replaced by a step change in the separator pressure (Figure 12).

For comments on formation pressure and pump speed inputs, see Section 5.1. Although the separator pressure is not normally an input function, there may be action applied to the separator that may create relatively quick changes in the separator pressure.

Figure 13 shows the pressures in front of the choke valves for the vertical pipes, as well as the manifold pressure.

Figure 13 demonstrates the positive pressure drop over the choke valves, and that they are different for the two valves ($\Delta p_{\mathrm{c}}^j=p_{\mathrm{c}}^{\mathrm{i},j}-p_{\mathrm{m}}$). Therefore, one should expect different flows through the two valves. Because the manifold pressure is a dependent (dynamic) variable in this case, there is no discontinuity in the pressures of Figure 13.

Figure 14 displays the pressure increase over the ESP pumps in the two vertical pipes.

In this figure, the sudden drop in $\Delta p_{\mathrm{p}}$ is due to a sudden change in the pump speed $f_{\mathrm{p}}$, and is thus realistic.

Figure 15 shows the friction pressure drops in the two vertical pipes and in the horizontal pipe.

Here, the interesting thing is that there is no visible difference between OpenModelica and Julia solution of $\Delta p_{\mathrm{f}}$ for the three pipes; confer Figure 11. Of course, this could be due to the zoomed-out view, but also on close inspection of $\Delta p_{\mathrm{f}}$ for the initial few seconds for vertical pipe 1, there is virtually no difference.

The resulting time constants and overall behavior in Figure 13, Figure 14 and Figure 15 are similar to those in Sharma [3]. Particularly in Figure 13 and Figure 14, some oscillatory behavior/overshoot is noticeable. This is to be expected due to the elasticity of the oil/water mixture with the given non-zero value of isothermal compressibility ${\beta }_T$.

Figure 16 shows vertical flow rates from reservoir to manifold in the two pipes, as well as the flow from manifold to separator (thick, solid lines), and the effect of uncertain productivity indices in Well 1, ${\dot{V}}_{\mathrm{pi}}^{\mathrm{c},1}\sim \mathcal{N}\left(7·{10}^{-4},{10}^{-4}\right)$, and uncertain isothermal compressibility in the petroleum fluid, ${\beta }_T\sim {\mathcal{U}}_{\left[0.3/1.5·{10}^9,3/1.5·{10}^9)\right]}$.

ModelingToolkit has support for efficient Monte Carlo studies; this is comparatively more complicated using Modelica + OMJulia.

5.3. Linearized model

ModelingToolkit.jl has good support for linearization of models. To linearize a system, it is necessary to provide a model where the inputs have not been defined (sys_p in Figure 17), a vector of input variables (sys_in), a vector of output variables (sys_in), and an operating point (keyword op, value op_0). If the ModelingToolkitStandardLibrary.jl (similar to Modelica’s Standard Library) is used, alternatively, some “virtual” inputs/outputs can be added in the form of Analysis Points which simplifies linearization; this is not discussed further here.

In Figure 17, linearization is performed using the named_ss function in ControlSystemsMTK.jl, which has similar arguments as function linearize in ModelingToolkit.jl.

Figure 17 suggests 6 poles/states in the system, and 3 transmission zeros. Obviously, the system has 4 states/differential variables (flow rates ${\dot{V}}_{\mathrm{v}}$ for each of the vertical pipes, pressure $p_{\mathrm{m}}$ for the manifold, and flow rate ${\dot{V}}_{\mathrm{t}}$ for the horizontal transport pipe). Comparing poles and zeros, one might suspect that two spurious/“infinitesimal” poles have been added together with two spurious/“infinitesimal” zeros, and that canceling out the tiny poles and zeros should give the correct transfer function. Observe also that there is a finite zero at $-7.27$ that cancels out one of the 4 true eigenvalues.

It is possible to write one’s own linearization code using a (symbolic) Jacobian function applicable to ModelingToolkit.jl models. If one assumes that the original model is a DAE of index 0 or 1, this will indeed give the correct transfer function with 4 states (and one zero that cancels out one of the eigenvalues). Why does ModelingToolkit.jl produce spurious additional poles/zeros? Possibly because the linearization algorithm in ModelingToolkit.jl makes no assumption of the index of the DAE, thus producing the two spurious “infinitesimal” poles/zeros.

Figure 18 shows the Bode plot of the transfer function from ESP rotational speed ($f_{\mathrm{p}}$)18 to flow rate into the separator (${\dot{V}}_{\mathrm{t}}$), as found using Julia + ModelingToolkit.jl based on the transfer function provided by the system in Figure 17. Package ControlSystems.jl also has a function balance_statespace() which in combination with minimal realization provides a transfer function with 3 poles or 1 time constant and one damped resonator,

$$P\left(s\right)\approx \frac{1.094·{10}^{-3}}{\left(1+0.15s\right)\left(1+2·0.489\frac{s}{1.06}+{\left(\frac{s}{1.06}\right)}^2\right)};$$

(60)

see pale curves in Figure 18.

It is possible to instead use Modelica+OMJulia for linearization, and then use ControlSystems.jl for Julia [similar capabilities as MATLAB’s Control Toolbox] for plotting and analysis/design. However, control analysis and design is simpler to do if a Julia set-up is used also for modeling and simulation.

5.4. Single-loop controller tuning

Figure 19 shows a unit step response for the linearized model using convenience function step() in package ControlSystems.jl.

A crude approximation of the system is read off Figure 19 as the plant transfer function $P_{\approx }\left(s\right)$

$$P_{\approx }\left(s\right)=K\frac{exp\left(-\tau s\right)}{s}$$

(61)

with

$$\begin{array}{cc}\hfill K& \approx 6.27·{10}^{-3}\hfill \\ \hfill \tau & \approx 0.5.\hfill \end{array}$$

A PI controller

$$C\left(s\right)=K_{\mathrm{p}}\frac{1+T_{\mathrm{i}}s}{T_{\mathrm{i}}s}$$

(62)

based on Skogestad’s method [15,16,17] is used with

$$\begin{array}{cc}\hfill T_{\ell }& =1.3·\tau \hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill \kappa & =4\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill K_{\mathrm{p}}& =\frac{1}{K\left(\tau +T_{\ell }\right)}\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill T_{\mathrm{i}}& =\kappa \left(\tau +T_{\ell }\right);\hfill \end{array}$$

(66)

where tuning parameters are: $T_{\ell }$, closed loop time constant; $\kappa$, integral time modifier.

Figure 20 shows a unit reference step response for the closed loop system [blue line]; here, utility function feedback(P*C) has been used to construct the closed-loop system. Observe that the closed loop time constant $T_{\ell }$ is not achieved, and that the resulting system is rather oscillatoric.

The reason for the oscillatoric behavior is that single-loop tuning methods are not designed for oscillatoric systems as the one in Eq. 60.

5.5. Double-loop controller tuning

In order to better handle the oscillatoric behavior of the system, which has generic form

$$P_{\zeta }\left(s\right)=\frac{K}{\left(1+T_{1}s\right)\left(1+2\zeta \frac{s}{{\omega }_0}+{\left(\frac{s}{{\omega }_0}\right)}^2\right)},$$

(67)

an idea of Nandong [18] is pursued, where the control signal is split into an “inner” control signal

$$u_{\mathrm{i}}=-C_{\mathrm{i}}\left(s\right)y$$

(68)

and an “outer” control signal $u_{\mathrm{o}}$,

$$u=u_{\mathrm{i}}+u_{\mathrm{o}}$$

(69)

where first $C_{\mathrm{i}}$ is designed to reduce the oscillations in the system.

Consider the following proper inner controller

$$C_{\mathrm{i}}\left(s\right)=K_{\mathrm{c}}^{\mathrm{i}}\frac{1+T_{1}s}{1+\alpha T_{1}s};$$

(70)

Inserting controller Eq. 68 into $y=P_{\zeta }\left(s\right)u$ with u as in Eq. 69 and $P_{\zeta }\left(s\right)$ as in Eq. 67 leads to,

$$\begin{array}{cc}\hfill y& =P_{\zeta }\left(s\right)u=P_{\zeta }\left(s\right)\left(-C_{\mathrm{i}}\left(s\right)y+u_{\mathrm{o}}\right)\hfill \\ \hfill & \Downarrow \hfill \\ \hfill \left(1+P_{\zeta }\left(s\right)C\left(s\right)\right)y& =P\left(s\right)u_{\mathrm{o}}\hfill \end{array}$$

which after some re-arrangement gives:

$$\left(1+\frac{K{K}_{\mathrm{c}}^{\mathrm{i}}}{1+\alpha T_{1}s}+2\zeta \frac{s}{{\omega }_0}+{\left(\frac{s}{{\omega }_0}\right)}^2\right)y=\frac{K}{1+T_{1}s}{u}_{\mathrm{o}}.$$

For small values of $\alpha$, $\alpha \to 0$:

$$\begin{array}{cc}\hfill \left(1+K{K}_{\mathrm{c}}^{\mathrm{i}}\right)+2\zeta \frac{s}{{\omega }_0}& +{\left(\frac{s}{{\omega }_0}\right)}^2=\left(1+K{K}_{\mathrm{c}}^{\mathrm{i}}\right)\times \hfill \\ \hfill & \cdots \times \left(1+2{\zeta }_{\mathrm{i}}\frac{s}{{\omega }_{\mathrm{i}}}+{\left(\frac{s}{{\omega }_{\mathrm{i}}}\right)}^2\right)\hfill \end{array}$$

where closed loop damping ${\zeta }_{\mathrm{i}}$ given by

$$\begin{array}{cc}\hfill {\zeta }_{\mathrm{i}}& =\zeta /\sqrt{1+K{K}_{\mathrm{c}}^{\mathrm{i}}}\hfill \\ \hfill & \Downarrow \hfill \\ \hfill K_{\mathrm{c}}^{\mathrm{i}}& =\frac{{\left(\zeta /{\zeta }_{\mathrm{i}}\right)}^2-1}{K}.\hfill \end{array}$$

(71)

A design procedure could thus be:

• Specify inner loop damping, ${\zeta }_{\mathrm{i}}\ge 1$ to provide (over-) damping.
• Compute inner gain $K_{\mathrm{c}}^{\mathrm{i}}$ from Eq. 71.
• Choose a “small” value for $\alpha$ to make the above design valid.

We choose ${\zeta }_{\mathrm{i}}=1$, and compute $K_{\mathrm{c}}^{\mathrm{i}}$ from Eq. 71. Next, closing the loop from $u_{\mathrm{o}}$ to y where we allow for $\alpha >0$ to have a realizable controller, the step response is as in Figure 21.

Based on closed inner loop with $\alpha =0.1$ as in the lower plot of Figure 21, we find model parameters in the model of Eq. 61 to be:

$$\begin{array}{cc}\hfill K& \approx 0.762·{10}^{-3}\hfill \\ \hfill \tau & \approx 0.5.\hfill \end{array}$$

Tuning an outer loop PI controller with Skogestad’s method with an integrator + delay approximation, this time with $T_{\ell }=2·\tau$ , the result is as in Figure 20, red curve. Although the result is considerable faster than what is achieved with the single-loop controller (same figure, blue curve), the result is surprisingly oscillatoric.

It is interesting to observe that with the same value for $T_{\ell }=2·\tau$, and changing $\alpha$ from $\alpha =0.1$ to $\alpha =5$, this gives a considerably better response, Figure 20, green curve.

5.6. Controller implementation with nonlinear model

The inner loop controller of Eq. 68 with $C_{\mathrm{i}}\left(s\right)$ as in Eq. 70 admits the following state space realization

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{z}_{\mathrm{i}}}{\mathrm{d}t}& =-\frac{1}{\alpha T_1}{z}_{\mathrm{i}}+y\hfill \\ \hfill u_{\mathrm{i}}& =\frac{K_{\mathrm{c}}^{\mathrm{i}}}{{\alpha }^2{T}_1}\left(1-\alpha \right)z_{\mathrm{i}}-\frac{K_{\mathrm{c}}^{\mathrm{i}}}{\alpha }y,\hfill \end{array}$$

where we must require $\alpha >0$. The PI controller for the outer loop as in Eq. 62 admits the following state space realization:

$$\begin{array}{cc}\hfill e& =r-y\hfill \\ \hfill \frac{\mathrm{d}{z}_{\mathrm{o}}}{\mathrm{d}t}& =\frac{1}{T_{\mathrm{i}}^{\mathrm{o}}}e\hfill \\ \hfill u_{\mathrm{o}}& =K_{\mathrm{p}}^{\mathrm{o}}\left(e+z_{\mathrm{o}}\right).\hfill \end{array}$$

With $y\leftarrow {\dot{V}}_{\mathrm{t}}$ as input to the controller together with a reference value $r\leftarrow {\dot{V}}_{\mathrm{t}}^{\mathrm{ref}}$, the controller produces output $u=u_{\mathrm{i}}+u_{\mathrm{o}}\to f_{\mathrm{p}}$ and is straightforward to implement in ModelingToolkit (or Modelica) and connect with the reservoir-to-separator model. Using the inner-outer model with $\alpha =5$, starting in steady state and injecting a step change in ${\dot{V}}_{\mathrm{t}}^{\mathrm{ref}}$ gives the response as in Figure 22.

The response in Figure 22 [red curve] is similar to the corresponding response [green curve] in Figure 20.

The response in Figure 22 is achieved with an ESP pump speed as in Figure 23.

In reality, it is not possible to “require” more production than the reservoir can produce. To handle this correctly, a more complex interaction with the reservoir than the productivity index model is required.

6. Conclusions

Sharma and Glemmestad [1], Sharma [3] provided a novel model of an ESP-lifted oil production system from reservoir to separator. In Lie [2], the model of Sharma and Glemmestad [1] was provided with a more structured model development, and with dimensionless equipment models. Also, some information that is difficult to find in the original publications (e.g., ESP model parameters) was included. A basic comparison of equation based modeling languages Modelica and ModelingToolkit for Julia was given.

Here, the work in Lie [2] is extended in several directions: (i) a couple of model typos in [2] are corrected, (ii) some more fluid details are given, (iii) more details of the ESP model is provided for reference (max/min flow rates, power consumption, efficiency), (iv) an extended discussion of pump control inputs is included, (v) a more thorough discussion is given of advantages with dimensionless models, (vi) a more detailed comparison of model implementation in Modelica vs. ModelingToolkit is provided, (vii) with a comparison of simulation results using OpenModelica vs. ModelingToolkit with Julia, (viii) more simulation results are included for reference, (ix) model linearization using ModelingToolkit with ControlSystems.jl for Julia is discussed, (x) controller tuning for the linearized model is illustrated using Skogestad’s method applied to a simple integrator+time delay approximation, as well as a double-loop design to dampen oscillations, (xi) the double-loop controller is applied to the nonlinear ModelingToolkit model to confirm the controller tuning. Model parameters, inputs, and initial states are provided in Appendix A.

A key difference between the model description here and the original one [1] is the focus on dimensionless models, which greatly reduces the chance of errors. The inclusion of the ESP power consumption model and flow rate constraints/best-efficiency-point for reference, makes it possible with realistic studies of constraints in controller design. Although not implemented here, an extension of more realistic pump control inputs by inclusion of the pump-motor aggregate moment of inertia is discussed. For some systems, such an extension may be of interest.

Both Modelica and ModelingToolkit are suitable and free languages for structured model development. Modelica is a mature language and a good choice, while ModelingToolkit is a more general language with an extensive eco-system. Model implementation is very similar for the two languages; ModelingToolkit does support more general model classes such as distributed models, stochastic models, etc., and has an extensive tool-set for control design, optimization, model fitting, etc. The numeric solutions are virtually identical for the two implementations; a minor discrepancy in friction pressure loss in one case is most likely due to a difference in handling implicit algebraic equations.

Model linearization is provided by both OpenModelica and ModelingToolkit for Julia; in the case of ModelingToolkit, the eco-system of Julia for further analysis and design can directly be used, while for OpenModelica, the linearization can be carried out via a script language API (e.g., from Julia, Python, or MATLAB) and analyzed in the host language.

For illustration purpose, linearization of the system from ESP pump speed to volumetric separator flow using ModelingToolkit/ControlSystems.jl is carried out in Julia. Then, several linear controllers are developed based on Skogestad’s method for an integrator+time delay model approximation, as well for a double-loop controller to dampen oscillations in the system caused by the compressible fluid in the manifold. A re-tuning of the double-loop controller with Skogestad’s method in the outer loop gives best performance. Implementation of the best controller in ModelingToolkit confirms that the linear controller works well with the nonlinear reservoir-to-separator system under the simulated experimental conditions.

The model development and discussion indicates a number of possible extensions such as (a) more realistic properties (density, viscosity), (b) allowing for distributed density along pipes19, (c) adding a more realistic system for water dilution in the manifold, (d) inclusion of valves in manifold–separator pipes, (e) integration with reservoir models, (f) use for advanced control design with constraints, (g) use for optimization, (h) integration with a reservoir model to study how to handle stiffness, and more Such extensions will give more insight into the industrial usefulness of the model.