*2.1. Case h* ≤ *H*<sup>t</sup>

When the water level is up to the tip of the access tunnel or below the tip of the access tunnel, *m*<sup>w</sup> is given by *m*<sup>w</sup> = *ρA*t where is the slant height for *h* as shown in Figure 1a. *m*w is further expressed as

$$m\_{\rm w} = \rho A\_{\rm t} h \frac{L}{H}.\tag{10}$$

The pressure force *F*<sup>p</sup> is formulated based on the pressure difference at the manifold and the air pressure with an expression

$$F\_{\rm P} = (p\_{\rm m} - p\_{\rm c})A\_{\rm t.} \tag{11}$$

The frictional force *F*<sup>f</sup> is expressed as

$$F\_{\rm f} = F\_{\rm D,W} + F\_{\rm D,a} \tag{12}$$

where *F*D,w is the frictional force formulated for water flow inside the surge tank based on Darcy's friction factor for water, *f*D,w. Similarly, *F*D,a is the frictional force formulated for air flow inside the surge tank based on Darcy's friction factor for air, *f*D,a. Both *f*D,w and *f*D,a are calculated as in [9]. The general expression for Darcy's friction factor *f*<sup>D</sup> is based on Reynolds' number *<sup>N</sup>*Re = *<sup>ρ</sup>*|*v*|*<sup>D</sup> <sup>μ</sup>* and expressed as

$$f\_D = \begin{cases} \frac{64}{N\_{\rm Re}} & N\_{\rm Re} < 2100\\ aN\_{\rm Re}^3 + bN\_{\rm Re}^2 + cN\_{\rm Re} + d & 2100 \le N\_{\rm Re} \le 2300\\ \frac{1}{\left(2\log\_{10}\left(\frac{\delta}{3\cdot\gamma D} + \frac{5\cdot\gamma}{N\_{\rm Re}^{10}}\right)\right)^2} & N\_{\rm Re} > 2300 \end{cases}$$

where *μ* is the dynamic viscosity of the fluid, *ε* is the pipe roughness height. For the region 2100 ≤ *N*Re ≤ 2300, *f*<sup>D</sup> is calculated from a cubic interpolation, with the coefficients *a*, *b*, *c*, and *d*, differentiable at the boundaries. The final expression for *F*<sup>f</sup> is calculated as in [9] given as

$$F\_{\rm f} = \frac{1}{2} \rho v \mid v \mid \left( A\_{\rm W,W} \frac{f\_{\rm D,W}}{4} + A\_{\rm W,a} \frac{f\_{\rm D,a}}{4} \right) \tag{13}$$

where | *v* | preserves the fluid frictional force against both directions of flow; flow induced from the access tunnel towards the air chamber, and vice-versa. *A*w,w is the wetted area due to water flow inside the surge tank given by

$$A\_{\rm W,W} = \pi D\_{\rm t} \ell \tag{14}$$

and *A*w,a is the wetted area due to the air during adiabatic compression and rarefaction inside the surge tank, and expressed as

$$A\_{\rm w,a} = \pi[D(L - L\_{\rm t}) + D\_{\rm t}(L\_{\rm t} - \ell)].\tag{15}$$

*2.2. Case h* > *H*<sup>t</sup>

When the water level inside the surge tank is above the access tunnel expression for *m*w is formulated by summing the mass of water inside the access tunnel and the mass of water inside the air chamber, and is expressed as

$$m\_{\rm W} = \rho \left[ A\_{\rm t} L\_{\rm t} + A(\ell - L\_{\rm t}) \right]. \tag{16}$$

For - > *L*<sup>t</sup> we consider Figure 2 for finding the total pressure force *F*<sup>p</sup> in the direction of the flow. The calculation of the fluid frictional force is given in Figure 3. From Figure 2, the pressure force *F*<sup>p</sup> is calculated based on the junction pressure *p*<sup>j</sup> between the junction of the access tunnel and the air chamber. *p*<sup>j</sup> is expressed as the sum of the air pressure *p*<sup>c</sup> and the hydrostatic pressure due to the difference in liquid-level *h* − *H*t. The junction pressure is then expressed as

$$p\_{\rm b} = p\_{\rm c} + \rho g (\ell - L\_{\rm t}) \frac{H}{L} \tag{17}$$

which relates in the final expression for *F*<sup>f</sup> as

$$F\_{\mathbb{P}} = (p\_{\mathrm{m}} - p\_{\mathrm{j}})A\_{\mathrm{t}} + (p\_{\mathrm{j}} - p\_{\mathrm{c}})A. \tag{18}$$

From Figure 2, the overall fluid frictional force *F*<sup>f</sup> is calculated with an expression given as

$$F\_{\rm f} = F\_{\rm D,w} + F\_{\rm \phi} + F\_{\rm D,a} \tag{19}$$

where *F*D,w + *F*D,a is given as

$$F\_{\rm D,w} + F\_{\rm D,a} = \frac{1}{2} \rho v \mid v \mid \left( A\_{\rm w,w} \frac{f\_{\rm D,w}}{4} + A\_{\rm w,a} \frac{f\_{\rm D,a}}{4} \right)$$

where *A*w,w = *π*[*D*t*L*<sup>t</sup> + *D*(- − *L*t)] and *A*w,a = *πD*(*L* − -); the calculations were similarly performed as in Equations (14) and (15).

**Figure 2.** Considering junction pressure *p*<sup>j</sup> for evaluating the overall pressure force *F*<sup>p</sup> in the direction of flow. *pj* is the pressure calculated based on the sum of air pressure *p*<sup>c</sup> and hydrostatic pressure due to liquid-level *h* − *H*t.

**Figure 3.** Expressions for fluid frictional force *F*<sup>f</sup> considering (**a**) the square expansion type fitting for the flow towards the chamber through the access tunnel and (**b**) the square reduction type fitting for the flow through the chamber to the access tunnel. In the figures, *φ*se and *φ*sr are the generalized friction factors for the square expansion and the square reduction type fittings, respectively, taken from [15].

In Equation (19), *Fφ* is the fluid frictional force due to water flow from the access tunnel towards the air chamber, and vice-versa. *Fφ* can be expressed in terms of the pressure drop (alternatively can be expressed in terms of the head loss). When the water is flowing *from the access tunnel towards the air chamber*, we consider the pressure drop due to the *square expansion* type of fitting as shown in Figure 3a, and when the water is flowing *from the air chamber towards the access tunnel*, we consider the pressure drop due to the *square reduction* type of fitting as shown in Figure 3b. Thus, *Fφ* is calculated based on the generalized friction factors *φ*se for the square expansion type of fitting and *φ*sr for the square reduction type of fitting. Additionally, for both types of flows as shown in Figure 3, we assume an average cross-sectional area

$$
\bar{A} = \frac{A + A\_{\rm t}}{2}.
$$

If Δ*p<sup>φ</sup>* is the pressure drop due to the fittings, there exists a relationship between Δ*pφ*, the average kinetic energy of the fluid per volume *K* = <sup>1</sup> <sup>2</sup> *ρv* | *v* | and the friction factor *φ* = {*φ*se, *φ*sr}. The relationship between Δ*pφ*, *K*, and *φ* is given by

$$
\Delta p\_{\phi} = \phi \mathcal{K}^{\prime\prime\prime}.
$$

The pressure drop Δ*p<sup>φ</sup>* is related to *F<sup>φ</sup>* through the average cross-sectional area *A*¯ and given as

$$F\_{\Phi} \approx \Delta p\_{\Phi} \bar{A}$$

which can be further expressed as

$$F\_{\Phi} \approx \frac{1}{2} \rho \upsilon \mid \upsilon \mid \bar{A} \phi , \ \ \phi = \{ \phi\_{\text{sc}} , \phi\_{\text{sr}} \} .$$

The final expression for overall fluid frictional force *F*<sup>f</sup> is then given as

$$F\_{\rm f} \approx \frac{1}{2} \rho v \mid v \mid \left( A\_{\rm W,W} \frac{f\_{\rm D,W}}{4} + A\_{\rm W,a} \frac{f\_{\rm D,a}}{4} + \bar{A} \phi \right) \quad \phi = \{\phi\_{\rm sc}, \phi\_{\rm sr}\}.\tag{20}$$

This completes the expressions for variables *m*, *F*<sup>p</sup> and *F*<sup>f</sup> for the two scenarios of the liquid level inside the surge tank, viz., *h* ≤ *H*<sup>t</sup> and *h* > *H*t. To further complete the information of variables in Equation (6), the expression for *FV*˙ <sup>g</sup> is calculated as

$$F\_{\mathfrak{F}}^{V} = mg\frac{H}{L},\tag{21}$$

as shown in the flow diagram of Figure 1a. Finally, the mechanistic model of the ACST needs an expression for the average velocity *v* expressed as

$$
v = \frac{\dot{V}}{\overline{A}}.\tag{22}$$

Equations (1)–(6), in addition to other associated algebraic relations from Equations (7)–(22), represent a semi-explicit DAEs formulation for the ACST, and can be modeled in a equation-based modeling language like Modelica. The developed mechanistic model of the ACST is implemented in OpenHPL as a feature extension, and the case study was carried out for Torpa HPP.

#### **3. Case Study**

Figure 4a shows the layout diagram of Torpa HPP. Similarly, Figure 4b shows the simulation model of Torpa HPP created in OpenHPL. In Figure 4b, the reservoir model, the intake tunnel model, the penstock model, and the discharge model are developed as in [21]. A detailed model of the penstock considering water compressibility and pipe elasticity can be formulated from [22]. However, we consider the penstock model as a simple pipe model. Similarly, the Francis turbine mechanistic model for the case study is modeled as in [23]. The mechanistic model for the tailrace is taken as an exact mirror replica of the reservoir model.

The dimensions of the ACST shown in Figure 4a are found based on the piezometric diagram for Torpa HPP from [10]. The model developed in Section 2 is based on a cylindrical access tunnel and a cylindrical air chamber. Thus, the hydraulic diameters for the access tunnel *D*<sup>t</sup> and the air chamber *D* are evaluated based on the volume of air inside the chamber using the operating conditions. Table 1 shows the parameters and the operating conditions of the ACST for Torpa HPP.

**Figure 4.** (**a**) Layout diagram for Torpa HPP. Nominal head, nominal discharge, and nominal power output are 445 m, 40 m3/s and 150 MW, respectively. The ACST has air volume of 13,000 m3, initially pressurized at 41 · <sup>10</sup><sup>5</sup> Pa. Similarly, both of the headrace and tailrace tunnels are 7 m in diameter. Torpa HPP consists of two turbine units each rated at 75 MW with rated discharge at 20 m3/s. Torpa HPP also consists of a tailrace surge tank not shown in the figure. (**b**) Simulation model of Torpa HPP implemented in OpenHPL from the head reservoir to the tail reservoir.

For the model created in Figure 4b, it is of interest to:


#### *3.1. Simulation Versus Real Measurements*

Figure 5 shows the simulated versus real measurement for Torpa HPP. As shown in Figure 4b, *u*v1 and *u*v2 are the turbine valve signals for the turbine unit-1 and the turbine unit-2, respectively, for controlling the volumetric discharge through the turbines. The input turbine valve signal for unit-1 is given by

$$u\_{\rm v1} = \begin{cases} 0.68 & 0 < t \le 500 \,\text{s} \\ \frac{0.68}{50} (t - 550) + 0.98 & 500 \,\text{s} < t \le 550 \,\text{s} \\ 0.98 & 550 \,\text{s} < t \le 1200 \,\text{s} \end{cases}$$

and the input turbine valve signal for unit-2 is given by,

$$u\_{\nabla 2} = \begin{cases} 0.55 & 0 < t \le 500 \,\mathrm{s} \\ \frac{0.55}{50}(t - 550) + 0.93 & 500 \,\mathrm{s} < t \le 550 \,\mathrm{s} \\ 0.93 & 550 \,\mathrm{s} < t \le 1200 \,\mathrm{s} . \end{cases}$$

For inputs *u*v1 and *u*v2, the mechanical power outputs for the turbine unit-1 (Figure 5c) and the turbine unit-2 (Figure 5d), the turbines inlet pressure *p*tr (Figure 5e), and the air pressure inside the surge tank *p*c (Figure 5f) are recorded for 1200 s with the measurement samples taken at each second. The air pressure *p*<sup>c</sup> is measured using the pressure sensor PARO scientific 8DP000-S with an error of less than 0.01% of full scale of 6 Mpa, the turbine inlet pressure *p*tr is measured using the pressure sensor PARO scientific DIQ 73K with an error of less than 0.04% of full scale of 20 Mpa, and the measurements for the mechanical power outputs are provided by the plant owner from Torpa HPP. The information about Torpa HPP and its experimental procedures are taken from [24]. Figure 5 shows that the simulation corresponds well with the real measurements in the case of power productions from the turbines (Figure 5c,d). In the case of the turbine inlet pressure *p*tr (Figure 5e) there is an steady-state error of 0.6 bar for 0 < *t* ≤ 500 s. We believe that the steady-state error in *p*tr for 0 < *t* ≤ 500 s can be eradicated by the inclusion of detailed geometrical dimensions for the headrace tunnel. In this paper, the headrace tunnel is considerd with a simple slanted pipe geometry as shown in Figure 4a. Similar steady-state error can be seen in the case of the height of water level inside the ACST *h* (Figure 5g) with negligible error of 0.05 m. In the case of air pressure inside the ACST *p*c, the simulation and the measurement data are in good agreement. The measurement sampling rate in the case of water level *h*, air pressure *p*c, and turbine power outputs are slower and oscillatory because the data are only recorded after a minimum change in the measured value, which may be the reason for the steady-state errors and phase difference between the simulation and measurements shown in Figure 5c,d,f,g. In addition, in Figure 5f,g for 800 s < *t* ≤ 1200 s, the simulated values have poorly damped oscillation while the measurement quickly reaches a steady value. The simulated and the experimental dynamics of the variables (*p*c and *h*) are not captured well because of the slower and oscillatory sampling rate of the sensors. The simulation and the real measurements are matched by manual tuning of pipe roughness height of the headrace tunnel (*ε* ≈ 0.4 mm), hydraulic diameter of the access tunnel *D*<sup>t</sup> ≈ 15 m, and hydraulic diameter of the air chamber *D* ≈ 24 m.


**Table 1.** Parameters and operating conditions of the ACST for Torpa HPP.
