**2. Mathematical Model**

Figure 1 shows a typical configuration of an irregular profile in a pipeline which consists of *n* entrapped air pockets, *d* air valves, *p* pipes, and *o* drain valves located at low points for draining the system, *D* being the internal pipe diameter ( *m*); *A* being the cross section area of the pipe (m2); *f* being the Darcy–Weisbach friction factor; and *g* being the gravity acceleration (m/s2). *Lj* (*j* = 1, 2, ..., *p*) is the total length of each pipe, *AVm* (*m* = 1, 2, ..., *d*) are the air valves, and *Ks* (*s* = 1, 2, ..., *o*) represents the flow factor for each drain valve. The emptying process starts when the drain valves are opened, thus air valves start to admit air into the pipeline. After that, the air pocket *i* will begin to expand generating subatmospheric pressure. At the same time, the drainage of the water column starts until the entire pipeline is completely emptied. Figure 1 gives the evolution at time *t*, the length of the emptying column *Le*,*<sup>j</sup>* (*j* = 1, 2, ..., *p*), the absolute pressure of the air pocket *p*∗ *i* (*i* = 1, 2, ..., *n*), and the water velocity *ve*,*j* (*j* = 1, 2, ..., *p*). The expansion of the air pocket *i* can be computed as *xi* = *Lj* − *Le*,*<sup>j</sup>* + *Lj*+<sup>1</sup> − *Le*,*j*+1.

**Figure 1.** Scheme of entrapped air pockets in a pipeline with irregular profile while water empties.

The one-dimensional (1D) proposed model has the following assumptions: (1) water column has been modeled by the rigid model; (2) the Darcy–Weisbach equation was considered to evaluate friction losses; (3) the thermodynamic behaviour of the air pocket is analyzed using a polytropic model; and (4) the air–water interface is perpendicular to the pipe direction. The proposed model can be used for pipelines with small diameters and with hydraulic slopes enough to prevent downstream air intrusion where open-channel flow does not occur [20,24].

Under these hypotheses, the problem is modeled by the following set of equations.

### *2.1. Equations for the Water Phase*

• Mass oscillation equation:

> The rigid model was used in order to compute the evolution of the water column [11,25], considering that the elasticity of the air is much higher than the elasticity of the pipe and the water [12]. Applying the rigid model to the emptying column *j* and considering that the drain valve *s* joins pipes *Lj* and *Lj*−1, which is a common point in a pipeline, then:

$$\frac{d\upsilon\_{\varepsilon,\mathbf{j}}}{dt} = \frac{p\_{\mathrm{i}}^{\*} - p\_{\mathrm{atm}}^{\*}}{\rho\_{\mathrm{w}}L\_{\varepsilon,\mathbf{j}}} + g\frac{\Delta z\_{\varepsilon,\mathbf{j}}}{L\_{\varepsilon,\mathbf{j}}} - f\_{\mathrm{j}}\frac{\upsilon\_{\varepsilon,\mathbf{j}}|\upsilon\_{\varepsilon,\mathbf{j}}|}{2D} - \frac{gA^{2}(\upsilon\_{\varepsilon,\mathbf{j}} + \upsilon\_{\varepsilon,\mathbf{j}-1})|\upsilon\_{\varepsilon,\mathbf{j}} + \upsilon\_{\varepsilon,\mathbf{j}-1}|}{L\_{\varepsilon,\mathbf{j}}K\_{\mathrm{s}}^{2}},\tag{1}$$

where <sup>Δ</sup>*ze*,*<sup>j</sup>* is the elevation difference (m), *ρw* is the water density (kg/m3), *p*<sup>∗</sup>*atm* is the atmospheric pressure (101,325 Pa) and *g* is the gravity acceleration (m/s2).

The expression *hm*,*<sup>s</sup>* = *Q*<sup>2</sup>*w*,*<sup>s</sup>*/*K*2*<sup>s</sup>* was used to estimate the local loss of the drain valve *s* in Equation (1). *Ks* is the flow factor and *Qw*,*<sup>s</sup>* is the water discharge.


The position of the air–water interface is considered perpendicular to the pipe direction [11,18,26]. The continuity equation for the moving air-water interface *j* is:

$$\frac{dL\_{\mathfrak{e},\mathfrak{j}}}{dt} = -\upsilon\_{\mathfrak{e},\mathfrak{j}} \quad \left(L\_{\mathfrak{e},\mathfrak{j}} = L\_{\mathfrak{e},\mathfrak{j},0} - \int\_{0}^{t} \upsilon\_{\mathfrak{e},\mathfrak{j}} \mathrm{d}t\right) ,\tag{2}$$

where subscript 0 refers to the initial condition.

### *2.2. Equations for Air Pockets*

• Continuity equation [3,17]:

$$\frac{dm\_{a,i}}{dt} = \rho\_{a,nc}v\_{a,nc,m}A\_{adm,m} \tag{3}$$

and by deriving:

$$\frac{dm\_{a,i}}{dt} = \frac{d(\rho\_{a,i}V\_{a,i})}{dt} = \frac{d\rho\_{a,i}}{dt}V\_{a,i} + \frac{dV\_{a,i}}{dt}\rho\_{a,i} \tag{4}$$

where *ma*,*<sup>i</sup>* is the air mass and *Va*,*<sup>i</sup>* is the air volume of the air pocket *i*.

Due to the air pocket *i* located between pipes *Lj* and *Lj*+1, then *Va*,*<sup>i</sup>* = (*Lj* − *Le*,*j*)*<sup>A</sup>* + (*Lj*+<sup>1</sup> − *Le*,*j*+<sup>1</sup>)*<sup>A</sup>* and *dVa*,*<sup>i</sup>*/*dt* = <sup>−</sup>(*dLe*,*j*/*dt*)*<sup>A</sup>* − (*dLe*,*j*+1/*dt*)*<sup>A</sup>* = *<sup>A</sup>*(*ve*,*<sup>j</sup>* + *ve*,*j*+<sup>1</sup>), thus:

$$\frac{d\rho\_{a,j}}{dt} = \frac{\rho\_{a,nc}v\_{a,nc,m}A\_{adm,m} - (v\_{c,j+1} + v\_{c,j})A\rho\_{a,i}}{A\left(L\_j - L\_{c,j} + L\_{j+1} - L\_{c,j+1}\right)},\tag{5}$$

where *ρ<sup>a</sup>*,*<sup>i</sup>* is the air density of the air pocket *i*, *ρ<sup>a</sup>*,*nc* is the air density in normal conditions (1.205 kg/m3), *Aadm*,*<sup>m</sup>* is the cross section area (m2) of the air valve *m* and *va*,*nc*,*<sup>m</sup>* is the air velocity in normal conditions admitted by the air valve *m*.

• Expansion equation for the air pocket *i*:

> the thermodynamic behaviour of the air pocket [27] is treated by using a polytropic model [22,28]

$$\frac{dp\_i^\*}{dt} = -k \frac{p\_i^\*}{V\_{a,i}} \frac{dV\_{a,i}}{dt} + \frac{p\_i^\*}{V\_{a,i}} \frac{k}{\rho\_{a,i}} \frac{dm\_{a,i}}{dt} \,. \tag{6}$$

Considering the equations presented before, then:

$$\frac{dp\_{\dot{i}}^{\*}}{dt} = \frac{kp\_{\dot{i}}^{\*}}{A(L\_{\dot{j}} - L\_{\mathfrak{e},\dot{j}} + L\_{\mathfrak{j}+1} - L\_{\mathfrak{e},\dot{j}+1})} \left(\frac{\rho\_{a,nc}v\_{a,nc,\dot{i}}A\_{adm,m}}{\rho\_{a,\dot{i}}} - A(v\_{\mathfrak{e},\dot{j}+1} + v\_{\mathfrak{e},\dot{j}})\right),\tag{7}$$

where *k* is the polytropic coefficient: if *k* = 1.0, then the process is isothermal, but, if *k* = 1.4, the process is adiabatic.

In Equations (5) and (7), if the air pocket *i* is only located in pipe *Lj*, then *ve*,*j*+<sup>1</sup> = 0 and *Lj*+<sup>1</sup> = *Le*,*j*+<sup>1</sup> = 0.

• Air valve characterization:

The formulation proposed by Wylie and Streeter [23,29] was used to represent the air admission into the system for air valves. Ideally, the air valve *m* should be working in subsonic flow (*p*<sup>∗</sup>*atm* > *p*∗*i*> 0.528*p*<sup>∗</sup>*atm*), then:

$$Q\_{a,nc,m} = C\_{adm,m} A\_{adm,m} \sqrt{7 p\_{atm}^\* p\_{a,nc} \left[ \left( \frac{p\_i^\*}{p\_{atm}^\*} \right)^{1.4286} - \left( \frac{p\_i^\*}{p\_{atm}^\*} \right)^{1.714} \right]},\tag{8}$$

where *Qa*,*nc*,*<sup>m</sup>* is the air discharge in normal conditions admitted by the air valve *m* and *Qa*,*nc*,*<sup>m</sup>* = *va*,*nc*,*<sup>m</sup> Aadm*,*m*.

When air valves are not located in high points of the system, then the position of them should be specified. When the air–water interface reaches an air valve, then it starts to admit air inside the system.

In summary, a set of 2*p* + 2*n* + *d* equations describes the whole system. Together with the corresponding boundary conditions, it can be solved for the 2*p* + 2*n* + *d* unknowns *ve*,*j*, *Le*,*j*, *p*∗*i* , *ρ<sup>a</sup>*,*i*, *va*,*nc*,*<sup>m</sup>*(*j* = 1, 2, ..., *p*; *i* = 1, 2, ..., *n*; *m* = 1, 2, ..., *d*).

### *2.3. Initial and Boundary Conditions*

When the pipeline is at rest, the initial conditions are described as follows: *ve*,*j*(0) = 0(*j* = 1, 2, .., *p*), *Le*,*j*(0) = *Le*,*j*,<sup>0</sup>(*j* = 1, 2, .., *p*), *p*∗*i* (0) = *<sup>p</sup>*<sup>∗</sup>*i*,<sup>0</sup>(*<sup>i</sup>* = 1, 2, .., *<sup>n</sup>*), *<sup>ρ</sup><sup>a</sup>*,*<sup>i</sup>*(0) = *<sup>ρ</sup><sup>a</sup>*,*i*,<sup>0</sup>(*<sup>i</sup>* = 1, 2, .., *<sup>n</sup>*), and *va*,*nc*,*<sup>m</sup>*(0) = 0(*i* = 1, 2, .., *d*).

The upstream boundary condition is given by *p*∗*i* = *<sup>p</sup>*<sup>∗</sup>*i*,0, the downstream is given by the flow factor *Ks* of the drain valve *s*, and the atmospheric pressure *p*<sup>∗</sup>*atm* due to the free discharge.
