*2.4. Gravity Term*

Figure 2 describes the evolution of the gravity term (see Equation (1)) along the emptying column *j*. *Lj*,*<sup>r</sup>* (*r* = 1, 2, ..., *h*) is the length of the pipe reach *r*, and *Lj* is the total pipe length *j* (*Lj* = ∑*<sup>r</sup>*=*<sup>h</sup> <sup>r</sup>*=1 *Lj*,*<sup>r</sup>*). Subscript *u* is used to identify the pipe reaches (1 to *h*) where the air–water interface is located.

**Figure 2.** Evolution of the gravity term of the emptying column *j*.

The gravity term of the emptying column *j* is computed by:

• When the air–water interface has not arrived the last reach (*r* = *h*):

$$\frac{\Delta z\_{\mathfrak{e},j}}{L\_{\mathfrak{e},j}} = \frac{\sum\_{r=\mathfrak{u}+1}^{r-h} L\_{j,r} \sin(\theta\_{j,r})}{L\_{\mathfrak{e},j}} + \left(1 - \frac{\sum\_{r=\mathfrak{u}+1}^{r-h} L\_{j,r}}{L\_{\mathfrak{e},j}}\right) \sin(\theta\_{j,\mathfrak{u}}) \,. \tag{9}$$

• When the air–water interface is located on the last reach (*r* = *h*):

$$\frac{\Delta z\_{\mathfrak{e},\mathfrak{j}}}{L\_{\mathfrak{e},\mathfrak{j}}} = \sin(\theta\_{\mathfrak{j},h}).\tag{10}$$

### **3. Application to Two Emptying Columns**

Figure 3 presents a case of two emptying columns. More complex systems can be treated in the same way based on the proposed model.

**Figure 3.** Two emptying columns in a pipeline.

The corresponding equations of the pipeline are: *Water* **2017**, *9*, 98

1. Mass oscillation equation applied to the emptying column 1

$$\frac{dv\_{\varepsilon,1}}{dt} = \frac{p\_1^\*-p\_{atm}^\*}{\rho\_{\rm w}L\_{\varepsilon,1}} + g\frac{\Delta z\_{\varepsilon,1}}{L\_{\varepsilon,1}} - f\_1\frac{v\_{\varepsilon,1}|v\_{\varepsilon,1}|}{2D} - \frac{gA\_1^2v\_{\varepsilon,1}|v\_{\varepsilon,1}|}{K\_1^2L\_{\varepsilon,1}}.\tag{11}$$

2. Emptying column 1 position

$$\frac{dL\_{\mathbf{c},1}}{dt} = -\upsilon\_{\mathbf{c},1} \quad \left(L\_{\mathbf{c},1} = L\_{\mathbf{c},1,0} - \int\_0^t \upsilon\_1 \mathbf{d}t\right). \tag{12}$$

3. Mass oscillation equation applied to the emptying column 2

$$\frac{dv\_{\varepsilon,2}}{dt} = \frac{p\_1^\*-p\_{\text{atm}}^\*}{\rho\_{\text{w}}L\_{\varepsilon,2}} + g\frac{\Delta z\_{\varepsilon,2}}{L\_{\varepsilon,2}} - f\_2\frac{v\_{\varepsilon,2}|v\_{\varepsilon,2}|}{2D} - \frac{gA\_2^2 v\_{\varepsilon,2}|v\_{\varepsilon,2}|}{K\_2^2 L\_{\varepsilon,2}}.\tag{13}$$

4. Emptying column 2 position

$$\frac{dL\_{\varepsilon,2}}{dt} = -v\_{\varepsilon,2} \quad \left(L\_{\varepsilon,2} = L\_{\varepsilon,2,0} - \int\_0^t v\_{\varepsilon,2} \mathbf{d}t\right). \tag{14}$$

5. Evolution of the air pocket

$$\frac{dp\_1^\*}{dt} = -k \frac{p\_1^\*(\upsilon\_{c,1}A\_1 + \upsilon\_{c,2}A\_2)}{A\_2(L\_2 - L\_{c,2}) + A\_1(L\_1 - L\_{c,1})} + \frac{p\_1^\* \rho\_{a,nc} \upsilon\_{a,nc,1} A\_{adm,1}}{A\_2(L\_2 - L\_{c,2}) + A\_1(L\_1 - L\_{c,1})} \frac{k}{\rho\_{a,1}}.\tag{15}$$

6. Continuity equation of the air pocket

$$\frac{d\rho\_{a,1}}{dt} = \frac{\rho\_{a,\rm nc}\upsilon\_{a,\rm nc,1}A\_{\rm adm,1} - (\upsilon\_{c,1}A\_1 + \upsilon\_{c,2}A\_2)\rho\_{a,1}}{A\_2(L\_2 - L\_{c,2}) + A\_1(L\_1 - L\_{c,1})}.\tag{16}$$

7. Air valve characterization

$$Q\_{a,nc,1} = C\_{adm,1} A\_{adm,1} \sqrt{T p\_{atm}^\* p\_{a,nc} \left[ \left( \frac{p\_1^\*}{p\_{atm}^\*} \right)^{1.4286} - \left( \frac{p\_1^\*}{p\_{atm}^\*} \right)^{1.714} \right]}.\tag{17}$$

This set of seven differential-algebraic equations (Equations (11) to (17)), together with the initial condition given by *ve*,<sup>1</sup>(0) = 0, *Le*,<sup>1</sup>(0) = *Le*,1,0, *ve*,<sup>2</sup>(0) = 0, *Le*,<sup>2</sup>(0) = *Le*,2,0, *<sup>p</sup>*<sup>∗</sup>1,0 = *p*<sup>∗</sup>*atm* = 101,325 Pa, *ρ<sup>a</sup>*,1,0 = *ρ<sup>a</sup>*,*nc* = 1.205 kg/m<sup>3</sup> and *Qa*,*nc*,<sup>1</sup>(0) = 0, allow for describing the whole problem. The Simulink Library in Matlab ( The MathWorks, Inc., Natick, MA, USA) was used in order to compute the seven unknown functions: *ve*,1, *Le*,1, *ve*,2, *Le*,2, *p*<sup>∗</sup>1, *ρ<sup>a</sup>*,<sup>1</sup> and *Qa*,*nc*,1.
