Electric Submersible Pump Lifted Oil Field: Basic Model for Control, and Comparison of Simulation Tools

Abstract

Optimal operation of petroleum production is important in a transition from energy systems based on fossil fuel to sustainable systems. One sub-process in petroleum production deals with transport from the (subsea) well-bore to a topside separator. Good control design for such operation requires a dynamic model of the petroleum flow from the well-bore to the separator. Here, such a dynamic model is considered for liquid production (oil/water) using an electric submersible pump (ESP) to aid in counteracting gravity and friction forces. Based on an existing model used for industrial control design, one goal is to report a complete dynamic model in a single paper. Emphasis is put on dimensionless equipment models for the simple change of units, and the model is developed from physical laws for easy replacement of sub-models, if needed. All the necessary information (equations, parameters) for model implementation is provided, and two candidate equation-based modeling languages are selected and compared: Modelica and ModelingToolkit [MTK] for Julia. The simulation results are virtually identical for the two languages and make sense from physics; however, there is a minor discrepancy in one plot—likely caused by slight differences in accuracy in handling initialization in the implicit algebraic equations. The implementation structures of the model in Modelica and MTK are similar. Modelica is a mature and excellent modeling tool, handles large-scale models, and has tools for producing C code and integration with other tools. MTK is still in rapid development, supports more model types than Modelica, and is integrated in an eco-system with excellent support for control design, optimization, model fitting, and more. To illustrate the suitability of using the developed model for control design, a simple PI controller is designed within the eco-system of MTK/Julia.

1. Introduction

This paper builds on material in Lie [1]. Petroleum products have been key energy carriers for more than a century. The current focus on climate (https://sdgs.un.org/goals; accessed on 13 December 2023) 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 the on-going research project “DigiWell” (DigiWell: see Funding). The vertical transport of petroleum from the oil well to the 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 the 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 the reservoir to the 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 simple, yet stringent, model development, while avoiding unit complexities in the variables.

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

Sharma’s model considers a case with four vertical pipes from the oil reservoirs to a single manifold, with two 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. The ESP model in Sharma and Glemmestad [2] is novel—the authors use a simple model of a booster pump, and use a valve model based on the ANSI/ISA S75.01 standard (http://integrated.cc/cse/ISA_750101_SPBd.pdf; accessed on 13 December 2023). The model in Sharma and Glemmestad [2] was re-structured and simplified in Lie [1], emphasizing dimensionless equipment models, thereby eliminating some of the complexity found in common industry models.

The dynamic model with ESP in Sharma and Glemmestad [2] is mainly relevant for the production of heavy oil. Several papers use this model in advanced industrial control studies, including Krishnamoorthy et al. [5], Santana et al. [6], Delou 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 the water content (Justiniano and Romero [8]).  Sharma and Glemmestad [2] assume an unrealistic linear viscosity dependence on the water fraction in pipe friction. Creation of an emulsion will likely modify the behavior of the ESP. Furthermore, real systems will contain some gas, which will also modify the behavior of the ESP (see, e.g., Zhu et al. [9,10,11]).

Although control-relevant models of ESP-aided oil production from bore-well heel to top-side separator are not new, it is difficult to find presentations with complete information provided in a concise way in a single location. This paper aims to fill that gap, while also demonstrating the advantage of using dimensionless equipment models. As a consequence, the paper contains all the necessary equations for a complete model, with the model parameters provided in an appendix. Although some models are over-simplified, the presentation should make it clear how more realistic models can replace the models used here. The simulation of realistic model topologies is greatly enhanced by using equation-based modeling languages, and a secondary aim is to demonstrate this. Simulation results are provided both to demonstrate the physical realism of the model, and also to allow for comparison with other implementations. Simulation languages embedded in an extensive eco-system will make control for optimal operation much simpler. A third aim is, thus, to demonstrate this possibility via a simple tuning of a controller based on the model at hand.

Section 2 gives an overview of the transport system from the oil reservoir via a manifold to a separator, and the 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 the liquid phase is considered, where the evaporation of liquids is assumed to be negligible. In the sequel, for an extensive variable x, $\dot{x}$ denotes the flow rate of x. More generally (i.e., for both extensive and intensive variables), $\frac{\mathrm{d}x}{\mathrm{d}t}$ denotes the time derivative of any variable x.

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—too large a diameter in pipes leads to unwanted laminar flow, while too small a diameter leads to high velocity and high friction loss. Water is commonly injected into the manifold to reduce the friction loss in the horizontal pipes. The added water is typically recycled from the separator, and is at close to the production temperature, thus making inclusion of the energy balance less important. 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 1, $p_{\mathrm{f}}^j$ is the reservoir formation pressure for well j, $p_{\mathrm{h}}^j$ is the bore-hole heel pressure, ${\dot{V}}_{\mathrm{h}}^j$ is the volumetric heel flow rate, $A_{\mathrm{v}}^{-,j}$ is the vertical cross-sectional area below the ESP, ${\ell }_{\mathrm{v}}^{-,j}$ is the length of the vertical pipe below the ESP, $p_{\mathrm{p}}^{\mathrm{i},j}$ is the inlet pressure to the (ESP) pump, $p_{\mathrm{p}}^{\mathrm{e},j}$ is the effluent pressure after the (ESP) pump, $A_{\mathrm{v}}^{+,j}$ is the vertical cross-sectional area above the ESP, ${\ell }_{\mathrm{v}}^{+,j}$ is the length of the vertical pipe above the ESP, ${\dot{V}}_{\mathrm{c}}^{\mathrm{i},j}$ is the volumetric flow rate at the inlet to the choke valve, $p_{\mathrm{c}}^{\mathrm{i},j}$ is the pressure at the inlet to the choke valve, and $p_{\mathrm{c}}^{\mathrm{e},j}$ is the pressure at the effluent from the choke valve. Next, $p_{\mathrm{m}}$ is the manifold pressure, while ${\dot{V}}_{\mathrm{w}}$ is the water volumetric flow rate added to the manifold to reduce the 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 the effluent pressure after the booster pump, ${\dot{V}}_{\mathrm{t}}^j$ is the volumetric flow rate in a transport pipe from the manifold to the separator, $p_{\mathrm{s}}^{\mathrm{i},j}$ is the 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 the vertical pipes are assumed to be connected to the same manifold pressure $p_{\mathrm{m}}$; hence, the 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 the transport pipes end up in the same separator: $p_{\mathrm{s}}^{\mathrm{e},j}=p_{\mathrm{s}}$ for all j. Here, a choke at the end of the transport pipes has been neglected.

2.2. Fluid Properties

The petroleum fluid properties are important. The density $\rho$ varies with the pressure p and the temperature T, $\rho \left(p,T\right)$. Neglecting the temperature dependence, and assuming constant isothermal compressibility ${\beta }_T$ (good thermodynamics books define the isothermal compressibility; see also https://en.wikipedia.org/wiki/Compressibility; accessed on 13 December 2023), the isothermal compressibility is the inverse of the bulk modulus.

$${\beta }_T=\frac{1}{\rho }{\left.\frac{\partial \rho }{\partial p}\right|}_T=\mathrm{constant}$$

leads to $\rho \left(p\right)$, 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 the water cut ${\chi }_{\mathrm{w}}$ as ${\chi }_{\mathrm{w}}\triangleq {\dot{V}}_{\mathrm{w}}/\dot{V}$, the volumetric flow rate of water divided by the total flow rate of the fluid, the total density $\rho$ can be expressed as

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

(2)

where ${\rho }_{\mathrm{w}}$ and ${\rho }_{\mathrm{o}}$ are the 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 the data in

Table A1. Parameters: petroleum liquid.

Parameter
${\beta }_T=\frac{1}{1.5\times {10}^9}\approx 6.67\times {10}^{-10}\mathrm{P}{\mathrm{a}}^{-1}$
$p_0=1\mathrm{bar}$
${\rho }_{\mathrm{o}}=900\mathrm{kg}/{\mathrm{m}}^3$
${\rho }_{\mathrm{w}}=1000\mathrm{kg}/{\mathrm{m}}^3$
${\chi }_{\mathrm{w}}=0.35$
${\rho }_0={\chi }_{\mathrm{w}}{\rho }_{\mathrm{w}}+\left(1-{\chi }_{\mathrm{w}}\right){\rho }_{\mathrm{o}}$
${\chi }_{\mathrm{w}}^{\mathrm{m}}=0.5$
${\rho }_0^{\mathrm{m}}={\chi }_{\mathrm{w}}^{\mathrm{m}}{\rho }_{\mathrm{w}}+\left(1-{\chi }_{\mathrm{w}}^{\mathrm{m}}\right){\rho }_{\mathrm{o}}$
${\nu }_{\mathrm{o}}=100\mathrm{cSt}=100\times {10}^{-6}{\mathrm{m}}^2/\mathrm{s}$
${\nu }_{\mathrm{w}}=1\mathrm{cSt}={10}^{-6}{\mathrm{m}}^2/\mathrm{s}$

of Appendix A, the density $\rho$ varies by ca. 10 kg/m³ with pressure variation in the range 25–225 bar (Figure 2).

Thus, we assume constant density in the pipes, but a pressure-dependent density is assumed in the manifold.

The authors of Sharma and Glemmestad [2] proposed a simple linear mixing rule for the kinematic viscosity $\nu$:   

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

(3)

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

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

The linear interpolation model of Equation (3) is used here to faciliate comparison with the results in Sharma and Glemmestad [2], even though it is not physically realistic [8].

2.3. Well-Bore Production

The total production from the reservoir (formation pressure $p_{\mathrm{f}}$) relates the 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 the 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 the scaling pressure with the same unit as $p_{\mathrm{f}},p_{\mathrm{h}}$.

2.4. Pump Models

2.4.1. Electric Submersible Pump

Pump models are often given as

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

(5)

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

The authors of Sharma and Glemmestad [2] provide values for the 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 [2], a comprehensive model for the pump head of a multi-stage ESP is developed. To ease change in the units, their model is rewritten here 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 Equation (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 }$ is the scaling flow rate, and $a_1,\dots ,a_3$ are dimensionless model parameters. (Here, $a_j$ is dimensionless, while in Sharma [3], the parameters $a_j$ have dimensions). This implies that the values of $a_j$ here are different from those of $a_j$ in Sharma and Glemmestad [2], Sharma [3]. In Sharma and Glemmestad [2], a head curve plot is provided; the result in Figure 3 based on the dimensionless model in Equation (9) is identical to their plot.

In addition, Sharma and Glemmestad [2] 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 pump (“brake horse power”, BHP, in the original publication), 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 Equation (10), ${\dot{W}}_{\mathrm{p},0}^{\mathrm{m}}$ is the 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 [2] include an efficiency curve plot; the result of Equation (12) based on a dimensionless model is identical to their plot.

2.4.2. Booster Pump

For the booster pump in the horizontal pipes, a simpler model is suggested in Sharma and Glemmestad [2], 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)

where $\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 units as $\Delta p_{\mathrm{bp},0}$, which is the pressure increase at the nominal pump frequency $f_{\mathrm{bp},0}$.

2.4.3. 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 the motor and the pump.

In Sharma [3], a model of the induction motor driving the ESP is developed. The experience of [3] is that the motor dynamics are much faster than those 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. These 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}}$$

with kinetic energy K

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

where $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}}$ is as in Equation (11).

In Sharma [3], the moment of inertia is given as approximately $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

The authors of Sharma and Glemmestad [2] base their valve models on the ANSI/ISA S75.01 standard (http://integrated.cc/cse/ISA_750101_SPBd.pdf; accessed on 13 December 2023). 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 the influent and effluent densities, respectively, $p_{\mathrm{i}},p_{\mathrm{e}}$ are the influent and effluent pressures, respectively, while ${\rho }^{\varsigma },p^{\varsigma }$ are the scaling density and pressure, respectively.

2.6. Friction Loss

The friction drop along the pipe can be given by the Darcy–Weisbach model (e.g., https://en.wikipedia.org/wiki/Darcy%E2%80%93Weisbach_equation; accessed on 13 December 2023) 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’s implicit expression (the Colebrook equation, sometimes known as the Colebrook–White equation). One explicit approximation to Colebrook’s expression is due to Swamee and Jain [12],

$$\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)

where $\mu$ is the dynamic viscosity, $\nu$ is the kinematic viscosity, and $ϵ$ is the “roughness height” of the pipe internal surface. The linear velocity v is related to the 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?

2.7.1. Example: ESP Pump Model

As a first example, consider the ESP model in Equation (9). In the original formulation in Sharma and Glemmestad [2] (slight change in notation),

$$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 the parameters ${\overline{a}}_j$ have rather complicated units and the equation is hard-coded with $\left[\dot{V}\right]=\mathrm{gal}/\min$, while $\left[h_{\mathrm{p}}\right]=\mathrm{ft}$ (for quantity x, $\left[x\right]$ is the unit of the quantity). In practice, (i) either the remainder of the model has to be posed using these units, or (ii) one has to operate with several copies of the variables, e.g., ${\dot{V}}_{\mathrm{gal}/\min }$ 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 the variables.

A much better solution is to write the model in dimensionless form. The simplest way to do this for the model in Equation (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{f}\mathrm{t}$ and ${\dot{V}}^{\varsigma }\equiv 1\mathrm{g}\mathrm{a}\mathrm{l}/\mathrm{m}\mathrm{i}\mathrm{n}$, then, $\left[{\tilde{a}}_j\right]\equiv \left[{\overline{a}}_j\right]$. 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}$ (the “native” unit is $\mathrm{gal}/\min$), this result is produced by choosing ${\dot{V}}^{\varsigma }=1\mathrm{gal}/\min =34.29\mathrm{bbl}/\mathrm{d}$, which can easily be found using, e.g., the WolframAlpha app (e.g., the Microsoft Store).

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 the parameters ${\overline{a}}_j$; for the parameters in Equation (9), the scaling parameters are in SI units, where ${\overline{a}}_j\to a_j$ and $a_j$ are given in

Table A2. Parameters: vertical pipe.

Parameter
${\ell }^{-}=100\mathrm{m}$
${\ell }^{+}=2000\mathrm{m}$
$d=0.1569\mathrm{m}$
$ϵ=0.0018\mathrm{inch}=45.7\mathsf{\mu }\mathrm{m}=45.7\times {10}^{-6}\mathrm{m}$
${\dot{V}}_{min,0}=14.43\times {10}^{-3}{\mathrm{m}}^3/\mathrm{s}$
${\dot{V}}_{\eta ,0}=19.83\times {10}^{-3}{\mathrm{m}}^3/\mathrm{s}$
${\dot{V}}_{max,0}=25.24\times {10}^{-3}{\mathrm{m}}^3/\mathrm{s}$
$h_{\mathrm{p},0}=1210.6\mathrm{m}$
$f_{\mathrm{p},0}=60\mathrm{Hz}$
${\dot{V}}^{\varsigma }=1{\mathrm{m}}^3/\mathrm{s}$
$a_1=-37.57$
$a_2=2.864\times {10}^3$
$a_3=-8.668\times {10}^4$
${\dot{W}}_{\mathrm{p},0}^{\mathrm{m}}=167.733\mathrm{k}\mathrm{W}$
$b_1=52.12$
$b_2=-768.7$
$b_3=38.544\times {10}^3$
$b_4=-1.534\times {10}^6$
${\dot{m}}_{\mathrm{v}}^{\mathrm{c}}=25.9\times {10}^3\mathrm{kg}/\mathrm{h}$
$f\left(u_{\mathrm{v}}\right)=\left\{\begin{array}{cc}0,\hfill & u_{\mathrm{v}}\le 0.05\hfill \\ \frac{11.1u_{\mathrm{v}}-0.556}{30},\hfill & 0.05<u_{\mathrm{v}}\le 0.5\hfill \\ \frac{50u_{\mathrm{v}}-20}{30},\hfill & 0.5<u_{\mathrm{v}}\le 1\hfill \end{array}\right.$
$p^{\varsigma }=1\mathrm{bar}$
${\rho }^{\varsigma }=1000\mathrm{kg}/{\mathrm{m}}^3$
${\dot{V}}_{\mathrm{pi}}^{\mathrm{c}}=7\times {10}^{-4}{\mathrm{m}}^3/\mathrm{s}$

. To find $a_j$, choose $h_{\mathrm{p},0}=1\mathrm{f}\mathrm{t}=0.3048\mathrm{m}$, ${\dot{V}}^{\varsigma }=1\mathrm{gal}/\min =6.309\times {10}^{-5}{\mathrm{m}}^3/\mathrm{s}$, and write Equation (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 Equation (9), while the 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 Equation (9).

2.7.2. Example: Control Valve

The ANSI/ISA S75.01 standard (http://integrated.cc/cse/ISA_750101_SPBd.pdf; accessed on 13 December 2023) 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)

where, 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, and $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 to the units rather complicated. The dimensionless formulation as in Equation (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) (see any good book on balance models, such as Bird et al. [13]). 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 the accumulated mass in the system, t is the time, $\dot{m}$ is the mass flow rate, and the 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 the linear momentum given as $\mathfrak{m}=mv$ with linear velocity v, $\dot{\mathfrak{m}}$ is the momentum flow rate given as $\dot{\mathfrak{m}}=\dot{m}v$, and F is the total force. With constant fluid density, ${\dot{\mathfrak{m}}}_{\mathrm{i}}={\dot{\mathfrak{m}}}_{\mathrm{e}}$, 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 the 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, Equation (23) reduces to Newton’s law. The momentum is given as $\mathfrak{m}=\dot{m}v$ with $\dot{m}=\rho {\dot{V}}_{\mathrm{v}}$, and v is related to ${\dot{V}}_{\mathrm{v}}$ by Equation (18). The total force is $F=F_{\mathrm{p}}+F_{\mathrm{b}}-F_{\mathrm{f}}-F_{\mathrm{g}}$, with

• Pressure forces at the 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 Equations (5) and (9).
• Friction loss,

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

  (26)

  with $\Delta p_{\mathrm{f}}$ given by Equations (15)–(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 the flow rate ${\dot{V}}_{\mathrm{v}}$ relates to the bottom hole pressure via the productivity index, Equation (4), and how the flow rate ${\dot{V}}_{\mathrm{v}}$ relates to the choke valve flow, Equation (14).

The most structured formulation would be to pose the momentum balance (here, Newton’s law) as the differential equation, and add all the necessary algebraic equations. However, the OpenModelica DAE solver struggles with such a formulation. The valve equation Equation (14) is implicit in the 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 (it was not tested whether ModelingToolkit can handle this implicit algebraic equation). Instead, if the differential variable is 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 the water cut ${\chi }_{\mathrm{w}}$, and (iii) all the 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 a constant manifold volume $V_{\mathrm{m}}$, and adding water at the 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 the 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)

The 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 the 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 Equation (44). Here, for simplicity, Equation (44) is still used.

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 the manifold to the separator, as well as the manifold water cut ${\chi }_{\mathrm{w}}^{\mathrm{m}}$, and (iii) the 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 the manifold to the separator is almost identical to the vertical pipe from the reservoir to the manifold. The essential differences are (i) no gravity pressure drop, (ii) a simpler booster pump model, (iii) neglect of the pressure drop from the pipe into the separator, and (iv) no need for a productivity 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 the manifold to the separator; the authors of Sharma and Glemmestad [2] use four 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 to 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 (two 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}}$, and the influent volumetric flows to the manifold should be the sum of the flows from the vertical pipes. 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, the connections should be performed 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 a general scripting language (Julia). Currently, MTK does not support a free 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 was solved using the free languages/tools OpenModelica [14,15] and ModelingToolkit [16] for Julia. To illustrate the similarity between the OpenModelica code and a current formulation using MTK, parts of the Modelica code and the MTK code for the reservoir heel–to–manifold are provided in Appendix B. The listings show that the Modelica code and the ModelingToolkit code have a high degree of similarity. A few things to note:

1.

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 the 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, e.g., Dt = Differential(t) or similar.

2.

In Modelica, the quantities need to be specified with a type (e.g., Real), and are prepended with a qualifier (e.g., constant, parameter) except for the variables. For Julia and ModelingToolkit, the data type is inferred, unless explicitly stated. In the code (Appendix B), quantities in MTK are grouped within begin...end blocks in macros (identifiers prepended by @, e.g., @parameters).

3.

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, the initial values for the unknowns after structural simplification (“states”) must be provided with numeric values (see variable p_c__i in the MTK code, Appendix B).

4.

In ModelingToolkit, the initial values of the 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.

5.

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 was posed as a DAE formulation with momentum as a differential variable. ModelingToolkit can use solvers from the large, high-quality DifferentialEquations.jl package [17]. With ModelingToolkit, more thought is currently required when choosing the 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 the step changes occur. On the other hand, the solutions from ModelingToolkit include interpolation functions, which yield smooth solutions with considerably fewer data points than for Modelica.

The 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 [18]. ModelingToolkit is integrated in the Julia eco-system, with support for linearization, plotting, control systems analysis, random variables, etc., and overall has more possibilities than 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 (on-going work on a JuliaSimCompiler.jl for a commercial extension of Julia will increase the possible system size). 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 two “states” and two 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 solvers (Julia’s DifferentialEquations.jl package can be accessed from Python and R). Neither MATLAB nor Python offer equation-based modeling languages with library/re-use support as do Modelica or ModelingToolkit; MathWorks does offer Simscape (commercial) with MATLAB integration for such use, however (https://se.mathworks.com/products/simscape.html; accessed on 13 December 2023).

5. Results

5.1. Reservoir Heel to Manifold

The parameters, initial conditions, and system inputs are given in Appendix A. Figure 4 shows the input variation for the reservoir heel–to–manifold ((R2M) case.

A step change in the formation pressure (red curve), as in Figure 4, is not very realistic; such changes are normally slow. The manifold pressure (blue curve) is not normally an input function, but rather is 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 rapidly. Nonetheless, the inputs in Figure 4 help provide useful information about time constants in the system.

Figure 5 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 4. 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 6 shows the response in the choke valve inlet pressure, $p_{\mathrm{c}}^{\mathrm{i}}$.

Observe that the choke valve inlet pressure for many types of step changes is continuous, but that it changes discontinuously upon a step change in manifold pressure at $t=$ 1.5 s. Again, this makes sense—when the manifold pressure drops (see Figure 4), the choke valve inlet pressure must also drop in proportion so that the flow through the valve changes continuously (see Equation (14)).

So far, ModelingToolkit/Julia and OpenModelica have provided virtually identical simulation results (Figure 5 and Figure 6). Figure 7 shows the pipe friction loss, $\Delta p_{\mathrm{f}}$.

It is interesting to contrast the magnitude of the friction loss in Figure 7 and how much smaller it is than the pressure boost in the ESP. Obviously, if a more realistic viscosity model than linear interpolation had been used, and in particular if emulsification occurs [8], the friction pressure drop might increase considerably 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 [17] and the solvers supported by OpenModelica. It appears that there is a minor discrepancy at $t=0$ for $\Delta p_{\mathrm{f}}$, which propagates throughout the solution. Because the codes and the 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 the initial time, and that $\Delta p_{\mathrm{f}}$ is rather sensitive to such an inaccuracy.

5.2. Reservoir Heel-to-Separator

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

Figure 8 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 4) has been replaced by a step change in the separator pressure (Figure 8).

For comments on the 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 9 shows the pressures in front of the choke valves for the vertical pipes, as well as the manifold pressure.

Figure 9 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 9. In general, for this extended system there is no visible difference between the OpenModelica and Julia solutions.

The resulting time constants and the overall behavior for the extended system (represented by Figure 9) are similar to those in Sharma [3]. Note that in Figure 9, some oscillatory behavior/overshoot is seen. This is to be expected due to the elasticity of the oil/water mixture with the given non-zero value of the isothermal compressibility ${\beta }_T$.

In order to understand the effect of parameter uncertainty in a model, Monte Carlo simulations are useful. Figure 10 shows the vertical flow rates from the reservoir to the manifold in the two pipes, as well as the flow from the manifold to the 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\times {10}^{-4},{10}^{-4}\right)$, and uncertain isothermal compressibility in the petroleum fluid, ${\beta }_T\sim {\mathcal{U}}_{\left[0.3/1.5\times {10}^9,3/1.5\times {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 the 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 11), 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 11, linearization is performed using the named_ss function in ControlSystemsMTK.jl, which has similar arguments as the function linearize in ModelingToolkit.jl.

Figure 11 suggests six poles/states in the system, and three transmission zeros. Obviously, the system has four 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 the 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 four 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 four 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 12 shows a Bode plot of the transfer function from the ESP rotational speed ($f_{\mathrm{p}}$) (it is assumed that the same speed is used for both ESPs in the vertical pipes) to the 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 11. The package ControlSystems.jl also has a function balance_statespace(), which, in combination with minimal realization, provides a transfer function with three poles or one time constant and one damped resonator,

$$P\left(s\right)\approx \frac{1.094\times {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 12).

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 also used for modeling and simulation.

5.4. Single-Loop Controller Tuning

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

A crude approximation of the system is read off Figure 13 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\times {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 [19,20,21] 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 the tuning parameters are $T_{\ell }$, the closed-loop time constant, and $\kappa$, the integral time modifier.

Figure 14 shows a unit reference step response for the closed-loop system [blue line]; here, the utility function feedback(P*C) was 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 such as the one described in Equation (60).

5.5. Double-Loop Controller Tuning

In order to better handle the oscillatoric behavior of the system, which has the 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 [22] 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 the 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 the controller Equation (68) into $y=P_{\zeta }\left(s\right)u$ with u, as in Equation (69), and $P_{\zeta }\left(s\right)$, as in Equation (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 the closed-loop damping ${\zeta }_{\mathrm{i}}$ is 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:

1.

Specify inner loop damping, ${\zeta }_{\mathrm{i}}\ge 1$ to provide (over-)damping.

2.

Compute inner gain $K_{\mathrm{c}}^{\mathrm{i}}$ from Equation (71).

3.

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 Equation (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 15.

Based on a closed inner loop with $\alpha =0.1$, as in the lower plot of Figure 15, we find the model parameters in the model of Equation (61) to be:

$$\begin{array}{cc}\hfill K& \approx 0.762\times {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 14 (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 (see Figure 14, green curve).

5.6. Controller Implementation with Non-Linear Model

The inner loop controller of Equation (68) with $C_{\mathrm{i}}\left(s\right)$ as in Equation (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 Equation (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 to connect with the reservoir-to-separator model. Using the inner–outer model with $\alpha =5$, starting in a steady state and injecting a step change in ${\dot{V}}_{\mathrm{t}}^{\mathrm{ref}}$ gives the response as shown in Figure 16.

The response in Figure 16 (red curve; nonlinear model) is similar to the corresponding response (green curve; linear approximation) in Figure 14.

The response in Figure 16 is achieved with an ESP pump speed as in Figure 17.

6. Conclusions

Sharma and Glemmestad [2], Sharma [3] provided a novel dynamic model of an ESP-lifted oil production system from the reservoir to the separator. Their pump models neglect the effect of the gas content and fluid viscosity variation, and they use a simplistic viscosity model in friction terms. Also, some of the information is spread over several publications, and the ESP model parameters are difficult to find.

Nevertheless, their model is used in the present work to facilitate comparison. Here, a structured development of the model in Sharma and Glemmestad [2] is provided from the fundamental balance laws with clearly stated simplifying assumptions. Their equipment models are rephrased in dimensionless form, and the merits of dimensionless forms are emphasized. A complete, rephrased mathematical model from the reservoir to the separator is provided, with all the necessary model parameters, initial states, and input functions. This achieves the goal of providing a complete model that can be implemented in a single presentation. Because of the structured model development, it is also relatively straightforward to replace the equipment models with alternative models.

The implementation of realistic and complex oil field topologies is greatly simplified by using modeling languages that facilitate use of model libraries and the re-use of models. The rephrased model of Sharma [3] is implemented in two such languages, Modelica and ModelingToolkit, for comparison. Modelica is a mature language for differential algebraic models, with excellent solvers that can handle large-scale models (in the range of ${10}^6$+ unknowns), and support for C code generation and integration via the FMI standard (FMI: functional mock-up interface, e.g., https://fmi-standard.org/; accessed on 13 December 2023) to other tools. ModelingToolkit [MTK] is newer, is in rapid development, handles a larger class of models (including distributed models), and is tightly integrated within the Julia scripting language, which has excellent support for control design, optimization, and more. It is demonstrated that the form of Modelica and MTK models can be very similar. The purpose of this work is not to implement a library unit of the model. However, the advantages of using modeling languages are clearly pointed out.

To illustrate the possibilities of modeling languages integrated with script languages, an example is given where the MTK model is linearized in Julia. Based on the linearized model, a simple PI controller from the ESP speed to the production flow rate is designed. This design is not meant for industrial use, but illustrates the simplicity of a control design workflow with modern modeling languages. The example also illustrates that standard PI tuning rules struggle with oscillatoric systems (here, due to the compressibility of the oil/water mixture), and that some additional refinement may improve the closed-loop system. However, one would anticipate that an even more complex controller, such as model predictive control [23], would further improve the performance.

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